| 1 | """Least squares interpolation. |
|---|
| 2 | |
|---|
| 3 | Implements a least-squares interpolation. |
|---|
| 4 | Putting mesh data onto points. |
|---|
| 5 | |
|---|
| 6 | Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou |
|---|
| 7 | Geoscience Australia, 2004. |
|---|
| 8 | |
|---|
| 9 | DESIGN ISSUES |
|---|
| 10 | * what variables should be global? |
|---|
| 11 | - if there are no global vars functions can be moved around alot easier |
|---|
| 12 | |
|---|
| 13 | * The public interface to Interpolate |
|---|
| 14 | __init__ |
|---|
| 15 | interpolate |
|---|
| 16 | interpolate_block |
|---|
| 17 | |
|---|
| 18 | """ |
|---|
| 19 | |
|---|
| 20 | import time |
|---|
| 21 | import os |
|---|
| 22 | from warnings import warn |
|---|
| 23 | from math import sqrt |
|---|
| 24 | from csv import writer, DictWriter |
|---|
| 25 | |
|---|
| 26 | from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \ |
|---|
| 27 | ArrayType, allclose, take, NewAxis, arange |
|---|
| 28 | |
|---|
| 29 | from caching.caching import cache |
|---|
| 30 | from pyvolution.neighbour_mesh import Mesh |
|---|
| 31 | from utilities.sparse import Sparse, Sparse_CSR |
|---|
| 32 | from utilities.cg_solve import conjugate_gradient, VectorShapeError |
|---|
| 33 | from coordinate_transforms.geo_reference import Geo_reference |
|---|
| 34 | from pyvolution.quad import build_quadtree |
|---|
| 35 | from utilities.numerical_tools import ensure_numeric, mean, NAN |
|---|
| 36 | from utilities.polygon import in_and_outside_polygon |
|---|
| 37 | from geospatial_data.geospatial_data import Geospatial_data, ensure_absolute |
|---|
| 38 | from fit_interpolate.search_functions import search_tree_of_vertices |
|---|
| 39 | from fit_interpolate.general_fit_interpolate import FitInterpolate |
|---|
| 40 | from pyvolution.util import file_function |
|---|
| 41 | |
|---|
| 42 | |
|---|
| 43 | class Interpolate (FitInterpolate): |
|---|
| 44 | |
|---|
| 45 | def __init__(self, |
|---|
| 46 | vertex_coordinates, |
|---|
| 47 | triangles, |
|---|
| 48 | mesh_origin=None, |
|---|
| 49 | verbose=False, |
|---|
| 50 | max_vertices_per_cell=30): |
|---|
| 51 | |
|---|
| 52 | |
|---|
| 53 | """ Build interpolation matrix mapping from |
|---|
| 54 | function values at vertices to function values at data points |
|---|
| 55 | |
|---|
| 56 | Inputs: |
|---|
| 57 | |
|---|
| 58 | vertex_coordinates: List of coordinate pairs [xi, eta] of |
|---|
| 59 | points constituting a mesh (or an m x 2 Numeric array or |
|---|
| 60 | a geospatial object) |
|---|
| 61 | Points may appear multiple times |
|---|
| 62 | (e.g. if vertices have discontinuities) |
|---|
| 63 | |
|---|
| 64 | triangles: List of 3-tuples (or a Numeric array) of |
|---|
| 65 | integers representing indices of all vertices in the mesh. |
|---|
| 66 | |
|---|
| 67 | mesh_origin: A geo_reference object or 3-tuples consisting of |
|---|
| 68 | UTM zone, easting and northing. |
|---|
| 69 | If specified vertex coordinates are assumed to be |
|---|
| 70 | relative to their respective origins. |
|---|
| 71 | |
|---|
| 72 | max_vertices_per_cell: Number of vertices in a quad tree cell |
|---|
| 73 | at which the cell is split into 4. |
|---|
| 74 | |
|---|
| 75 | Note: Don't supply a vertex coords as a geospatial object and |
|---|
| 76 | a mesh origin, since geospatial has its own mesh origin. |
|---|
| 77 | """ |
|---|
| 78 | |
|---|
| 79 | # Initialise variabels |
|---|
| 80 | self._A_can_be_reused = False |
|---|
| 81 | self._point_coordinates = None |
|---|
| 82 | |
|---|
| 83 | FitInterpolate.__init__(self, |
|---|
| 84 | vertex_coordinates, |
|---|
| 85 | triangles, |
|---|
| 86 | mesh_origin, |
|---|
| 87 | verbose, |
|---|
| 88 | max_vertices_per_cell) |
|---|
| 89 | |
|---|
| 90 | # FIXME: What is a good start_blocking_len value? |
|---|
| 91 | def interpolate(self, f, point_coordinates = None, |
|---|
| 92 | start_blocking_len = 500000, verbose=False): |
|---|
| 93 | """Interpolate mesh data f to determine values, z, at points. |
|---|
| 94 | |
|---|
| 95 | f is the data on the mesh vertices. |
|---|
| 96 | |
|---|
| 97 | The mesh values representing a smooth surface are |
|---|
| 98 | assumed to be specified in f. |
|---|
| 99 | |
|---|
| 100 | Inputs: |
|---|
| 101 | f: Vector or array of data at the mesh vertices. |
|---|
| 102 | If f is an array, interpolation will be done for each column as |
|---|
| 103 | per underlying matrix-matrix multiplication |
|---|
| 104 | point_coordinates: Interpolate mesh data to these positions. |
|---|
| 105 | List of coordinate pairs [x, y] of |
|---|
| 106 | data points or an nx2 Numeric array or a Geospatial_data object |
|---|
| 107 | |
|---|
| 108 | If point_coordinates is absent, the points inputted last time |
|---|
| 109 | this method was called are used, if possible. |
|---|
| 110 | start_blocking_len: If the # of points is more or greater than this, |
|---|
| 111 | start blocking |
|---|
| 112 | |
|---|
| 113 | Output: |
|---|
| 114 | Interpolated values at inputted points (z). |
|---|
| 115 | """ |
|---|
| 116 | #print "point_coordinates interpolate.interpolate",point_coordinates |
|---|
| 117 | if isinstance(point_coordinates,Geospatial_data): |
|---|
| 118 | point_coordinates = point_coordinates.get_data_points( \ |
|---|
| 119 | absolute = True) |
|---|
| 120 | |
|---|
| 121 | # Can I interpolate, based on previous point_coordinates? |
|---|
| 122 | if point_coordinates is None: |
|---|
| 123 | if self._A_can_be_reused is True and \ |
|---|
| 124 | len(self._point_coordinates) < start_blocking_len: |
|---|
| 125 | z = self._get_point_data_z(f, |
|---|
| 126 | verbose=verbose) |
|---|
| 127 | elif self._point_coordinates is not None: |
|---|
| 128 | # if verbose, give warning |
|---|
| 129 | if verbose: |
|---|
| 130 | print 'WARNING: Recalculating A matrix, due to blocking.' |
|---|
| 131 | point_coordinates = self._point_coordinates |
|---|
| 132 | else: |
|---|
| 133 | #There are no good point_coordinates. import sys; sys.exit() |
|---|
| 134 | msg = 'ERROR (interpolate.py): No point_coordinates inputted' |
|---|
| 135 | raise Exception(msg) |
|---|
| 136 | |
|---|
| 137 | |
|---|
| 138 | if point_coordinates is not None: |
|---|
| 139 | self._point_coordinates = point_coordinates |
|---|
| 140 | if len(point_coordinates) < start_blocking_len or \ |
|---|
| 141 | start_blocking_len == 0: |
|---|
| 142 | self._A_can_be_reused = True |
|---|
| 143 | z = self.interpolate_block(f, point_coordinates, |
|---|
| 144 | verbose=verbose) |
|---|
| 145 | else: |
|---|
| 146 | #Handle blocking |
|---|
| 147 | self._A_can_be_reused = False |
|---|
| 148 | start=0 |
|---|
| 149 | # creating a dummy array to concatenate to. |
|---|
| 150 | |
|---|
| 151 | f = ensure_numeric(f, Float) |
|---|
| 152 | #print "f.shape",f.shape |
|---|
| 153 | if len(f.shape) > 1: |
|---|
| 154 | z = zeros((0,f.shape[1])) |
|---|
| 155 | else: |
|---|
| 156 | z = zeros((0,)) |
|---|
| 157 | |
|---|
| 158 | for end in range(start_blocking_len |
|---|
| 159 | ,len(point_coordinates) |
|---|
| 160 | ,start_blocking_len): |
|---|
| 161 | t = self.interpolate_block(f, point_coordinates[start:end], |
|---|
| 162 | verbose=verbose) |
|---|
| 163 | #print "t", t |
|---|
| 164 | #print "z", z |
|---|
| 165 | z = concatenate((z,t)) |
|---|
| 166 | start = end |
|---|
| 167 | end = len(point_coordinates) |
|---|
| 168 | t = self.interpolate_block(f, point_coordinates[start:end], |
|---|
| 169 | verbose=verbose) |
|---|
| 170 | z = concatenate((z,t)) |
|---|
| 171 | return z |
|---|
| 172 | |
|---|
| 173 | def interpolate_block(self, f, point_coordinates = None, verbose=False): |
|---|
| 174 | """ |
|---|
| 175 | Call this if you want to control the blocking or make sure blocking |
|---|
| 176 | doesn't occur. |
|---|
| 177 | |
|---|
| 178 | Return the point data, z. |
|---|
| 179 | |
|---|
| 180 | See interpolate for doc info. |
|---|
| 181 | """ |
|---|
| 182 | if isinstance(point_coordinates,Geospatial_data): |
|---|
| 183 | point_coordinates = point_coordinates.get_data_points( \ |
|---|
| 184 | absolute = True) |
|---|
| 185 | if point_coordinates is not None: |
|---|
| 186 | self._A =self._build_interpolation_matrix_A(point_coordinates, |
|---|
| 187 | verbose=verbose) |
|---|
| 188 | return self._get_point_data_z(f) |
|---|
| 189 | |
|---|
| 190 | def _get_point_data_z(self, f, verbose=False): |
|---|
| 191 | """ |
|---|
| 192 | Return the point data, z. |
|---|
| 193 | |
|---|
| 194 | Precondition, |
|---|
| 195 | The _A matrix has been created |
|---|
| 196 | """ |
|---|
| 197 | z = self._A * f |
|---|
| 198 | # Taking into account points outside the mesh. |
|---|
| 199 | #print "self.outside_poly_indices", self.outside_poly_indices |
|---|
| 200 | #print "self.inside_poly_indices", self.inside_poly_indices |
|---|
| 201 | #print "z", z |
|---|
| 202 | for i in self.outside_poly_indices: |
|---|
| 203 | z[i] = NAN |
|---|
| 204 | return z |
|---|
| 205 | |
|---|
| 206 | |
|---|
| 207 | def _build_interpolation_matrix_A(self, |
|---|
| 208 | point_coordinates, |
|---|
| 209 | verbose = False): |
|---|
| 210 | """Build n x m interpolation matrix, where |
|---|
| 211 | n is the number of data points and |
|---|
| 212 | m is the number of basis functions phi_k (one per vertex) |
|---|
| 213 | |
|---|
| 214 | This algorithm uses a quad tree data structure for fast binning |
|---|
| 215 | of data points |
|---|
| 216 | origin is a 3-tuple consisting of UTM zone, easting and northing. |
|---|
| 217 | If specified coordinates are assumed to be relative to this origin. |
|---|
| 218 | |
|---|
| 219 | This one will override any data_origin that may be specified in |
|---|
| 220 | instance interpolation |
|---|
| 221 | |
|---|
| 222 | Preconditions |
|---|
| 223 | Point_coordindates and mesh vertices have the same origin. |
|---|
| 224 | """ |
|---|
| 225 | |
|---|
| 226 | |
|---|
| 227 | |
|---|
| 228 | #Convert point_coordinates to Numeric arrays, in case it was a list. |
|---|
| 229 | point_coordinates = ensure_numeric(point_coordinates, Float) |
|---|
| 230 | |
|---|
| 231 | if verbose: print 'Getting indices inside mesh boundary' |
|---|
| 232 | #print "self.mesh.get_boundary_polygon()",self.mesh.get_boundary_polygon() |
|---|
| 233 | self.inside_poly_indices, self.outside_poly_indices = \ |
|---|
| 234 | in_and_outside_polygon(point_coordinates, |
|---|
| 235 | self.mesh.get_boundary_polygon(), |
|---|
| 236 | closed = True, verbose = verbose) |
|---|
| 237 | #print "self.inside_poly_indices",self.inside_poly_indices |
|---|
| 238 | #print "self.outside_poly_indices",self.outside_poly_indices |
|---|
| 239 | #Build n x m interpolation matrix |
|---|
| 240 | if verbose and len(self.outside_poly_indices) >0: |
|---|
| 241 | print '\n WARNING: Points outside mesh boundary. \n' |
|---|
| 242 | # Since you can block, throw a warning, not an error. |
|---|
| 243 | if verbose and 0 == len(self.inside_poly_indices): |
|---|
| 244 | print '\n WARNING: No points within the mesh! \n' |
|---|
| 245 | |
|---|
| 246 | m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex) |
|---|
| 247 | n = point_coordinates.shape[0] #Nbr of data points |
|---|
| 248 | |
|---|
| 249 | if verbose: print 'Number of datapoints: %d' %n |
|---|
| 250 | if verbose: print 'Number of basis functions: %d' %m |
|---|
| 251 | |
|---|
| 252 | A = Sparse(n,m) |
|---|
| 253 | |
|---|
| 254 | n = len(self.inside_poly_indices) |
|---|
| 255 | #Compute matrix elements for points inside the mesh |
|---|
| 256 | for i in self.inside_poly_indices: |
|---|
| 257 | #For each data_coordinate point |
|---|
| 258 | if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n) |
|---|
| 259 | x = point_coordinates[i] |
|---|
| 260 | element_found, sigma0, sigma1, sigma2, k = \ |
|---|
| 261 | search_tree_of_vertices(self.root, self.mesh, x) |
|---|
| 262 | #Update interpolation matrix A if necessary |
|---|
| 263 | if element_found is True: |
|---|
| 264 | #Assign values to matrix A |
|---|
| 265 | |
|---|
| 266 | j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0 |
|---|
| 267 | j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1 |
|---|
| 268 | j2 = self.mesh.triangles[k,2] #Global vertex id for sigma2 |
|---|
| 269 | |
|---|
| 270 | sigmas = {j0:sigma0, j1:sigma1, j2:sigma2} |
|---|
| 271 | js = [j0,j1,j2] |
|---|
| 272 | |
|---|
| 273 | for j in js: |
|---|
| 274 | A[i,j] = sigmas[j] |
|---|
| 275 | else: |
|---|
| 276 | msg = 'Could not find triangle for point', x |
|---|
| 277 | raise Exception(msg) |
|---|
| 278 | return A |
|---|
| 279 | |
|---|
| 280 | def interpolate_sww2csv(sww_file, |
|---|
| 281 | points, |
|---|
| 282 | depth_file, |
|---|
| 283 | velocity_x_file, |
|---|
| 284 | velocity_y_file, |
|---|
| 285 | #quantities = ['depth', 'velocity'], |
|---|
| 286 | verbose=True, |
|---|
| 287 | use_cache = True): |
|---|
| 288 | """ |
|---|
| 289 | Interpolate the quantities at a given set of locations, given |
|---|
| 290 | an sww file. |
|---|
| 291 | The results are written to a csv file. |
|---|
| 292 | |
|---|
| 293 | In the future let points be a csv or xya file. |
|---|
| 294 | And the user choose the quantities. |
|---|
| 295 | |
|---|
| 296 | This is currently quite specific. |
|---|
| 297 | If it need to be more general, chagne things. |
|---|
| 298 | |
|---|
| 299 | This is really returning speed, not velocity. |
|---|
| 300 | """ |
|---|
| 301 | quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
|---|
| 302 | #print "points",points |
|---|
| 303 | points = ensure_absolute(points) |
|---|
| 304 | point_count = len(points) |
|---|
| 305 | callable_sww = file_function(sww_file, |
|---|
| 306 | quantities=quantities, |
|---|
| 307 | interpolation_points=points, |
|---|
| 308 | verbose=verbose, |
|---|
| 309 | use_cache=use_cache) |
|---|
| 310 | |
|---|
| 311 | depth_writer = writer(file(depth_file, "wb")) |
|---|
| 312 | velocity_x_writer = writer(file(velocity_x_file, "wb")) |
|---|
| 313 | velocity_y_writer = writer(file(velocity_y_file, "wb")) |
|---|
| 314 | # Write heading |
|---|
| 315 | heading = [str(x[0])+ ':' + str(x[1]) for x in points] |
|---|
| 316 | heading.insert(0, "time") |
|---|
| 317 | depth_writer.writerow(heading) |
|---|
| 318 | velocity_x_writer.writerow(heading) |
|---|
| 319 | velocity_y_writer.writerow(heading) |
|---|
| 320 | |
|---|
| 321 | for time in callable_sww.get_time(): |
|---|
| 322 | depths = [time] |
|---|
| 323 | velocity_xs = [time] |
|---|
| 324 | velocity_ys = [time] |
|---|
| 325 | for point_i, point in enumerate(points): |
|---|
| 326 | quantities = callable_sww(time,point_i) |
|---|
| 327 | #print "quantities", quantities |
|---|
| 328 | |
|---|
| 329 | w = quantities[0] |
|---|
| 330 | z = quantities[1] |
|---|
| 331 | momentum_x = quantities[2] |
|---|
| 332 | momentum_y = quantities[3] |
|---|
| 333 | depth = w - z |
|---|
| 334 | |
|---|
| 335 | if w == NAN or z == NAN or momentum_x == NAN: |
|---|
| 336 | velocity_x = NAN |
|---|
| 337 | else: |
|---|
| 338 | if depth > 1.e-30: # use epsilon |
|---|
| 339 | velocity_x = momentum_x / depth #Absolute velocity |
|---|
| 340 | else: |
|---|
| 341 | velocity_x = 0 |
|---|
| 342 | if w == NAN or z == NAN or momentum_y == NAN: |
|---|
| 343 | velocity_y = NAN |
|---|
| 344 | else: |
|---|
| 345 | if depth > 1.e-30: # use epsilon |
|---|
| 346 | velocity_y = momentum_y / depth #Absolute velocity |
|---|
| 347 | else: |
|---|
| 348 | velocity_y = 0 |
|---|
| 349 | depths.append(depth) |
|---|
| 350 | velocity_xs.append(velocity_x) |
|---|
| 351 | velocity_ys.append(velocity_y) |
|---|
| 352 | depth_writer.writerow(depths) |
|---|
| 353 | velocity_x_writer.writerow(velocity_xs) |
|---|
| 354 | velocity_y_writer.writerow(velocity_ys) |
|---|
| 355 | |
|---|
| 356 | class Interpolation_function: |
|---|
| 357 | """Interpolation_interface - creates callable object f(t, id) or f(t,x,y) |
|---|
| 358 | which is interpolated from time series defined at vertices of |
|---|
| 359 | triangular mesh (such as those stored in sww files) |
|---|
| 360 | |
|---|
| 361 | Let m be the number of vertices, n the number of triangles |
|---|
| 362 | and p the number of timesteps. |
|---|
| 363 | |
|---|
| 364 | Mandatory input |
|---|
| 365 | time: px1 array of monotonously increasing times (Float) |
|---|
| 366 | quantities: Dictionary of arrays or 1 array (Float) |
|---|
| 367 | The arrays must either have dimensions pxm or mx1. |
|---|
| 368 | The resulting function will be time dependent in |
|---|
| 369 | the former case while it will be constan with |
|---|
| 370 | respect to time in the latter case. |
|---|
| 371 | |
|---|
| 372 | Optional input: |
|---|
| 373 | quantity_names: List of keys into the quantities dictionary |
|---|
| 374 | vertex_coordinates: mx2 array of coordinates (Float) |
|---|
| 375 | triangles: nx3 array of indices into vertex_coordinates (Int) |
|---|
| 376 | interpolation_points: Nx2 array of coordinates to be interpolated to |
|---|
| 377 | verbose: Level of reporting |
|---|
| 378 | |
|---|
| 379 | |
|---|
| 380 | The quantities returned by the callable object are specified by |
|---|
| 381 | the list quantities which must contain the names of the |
|---|
| 382 | quantities to be returned and also reflect the order, e.g. for |
|---|
| 383 | the shallow water wave equation, on would have |
|---|
| 384 | quantities = ['stage', 'xmomentum', 'ymomentum'] |
|---|
| 385 | |
|---|
| 386 | The parameter interpolation_points decides at which points interpolated |
|---|
| 387 | quantities are to be computed whenever object is called. |
|---|
| 388 | If None, return average value |
|---|
| 389 | """ |
|---|
| 390 | |
|---|
| 391 | |
|---|
| 392 | |
|---|
| 393 | def __init__(self, |
|---|
| 394 | time, |
|---|
| 395 | quantities, |
|---|
| 396 | quantity_names = None, |
|---|
| 397 | vertex_coordinates = None, |
|---|
| 398 | triangles = None, |
|---|
| 399 | interpolation_points = None, |
|---|
| 400 | verbose = False): |
|---|
| 401 | """Initialise object and build spatial interpolation if required |
|---|
| 402 | """ |
|---|
| 403 | |
|---|
| 404 | from Numeric import array, zeros, Float, alltrue, concatenate,\ |
|---|
| 405 | reshape, ArrayType |
|---|
| 406 | |
|---|
| 407 | |
|---|
| 408 | #from util import mean, ensure_numeric |
|---|
| 409 | from config import time_format |
|---|
| 410 | import types |
|---|
| 411 | |
|---|
| 412 | |
|---|
| 413 | #Check temporal info |
|---|
| 414 | time = ensure_numeric(time) |
|---|
| 415 | msg = 'Time must be a monotonuosly ' |
|---|
| 416 | msg += 'increasing sequence %s' %time |
|---|
| 417 | assert alltrue(time[1:] - time[:-1] >= 0 ), msg |
|---|
| 418 | |
|---|
| 419 | |
|---|
| 420 | #Check if quantities is a single array only |
|---|
| 421 | if type(quantities) != types.DictType: |
|---|
| 422 | quantities = ensure_numeric(quantities) |
|---|
| 423 | quantity_names = ['Attribute'] |
|---|
| 424 | |
|---|
| 425 | #Make it a dictionary |
|---|
| 426 | quantities = {quantity_names[0]: quantities} |
|---|
| 427 | |
|---|
| 428 | |
|---|
| 429 | #Use keys if no names are specified |
|---|
| 430 | if quantity_names is None: |
|---|
| 431 | quantity_names = quantities.keys() |
|---|
| 432 | |
|---|
| 433 | |
|---|
| 434 | #Check spatial info |
|---|
| 435 | if vertex_coordinates is None: |
|---|
| 436 | self.spatial = False |
|---|
| 437 | else: |
|---|
| 438 | vertex_coordinates = ensure_numeric(vertex_coordinates) |
|---|
| 439 | |
|---|
| 440 | assert triangles is not None, 'Triangles array must be specified' |
|---|
| 441 | triangles = ensure_numeric(triangles) |
|---|
| 442 | self.spatial = True |
|---|
| 443 | |
|---|
| 444 | |
|---|
| 445 | #Save for use with statistics |
|---|
| 446 | |
|---|
| 447 | self.quantities_range = {} |
|---|
| 448 | for name in quantity_names: |
|---|
| 449 | q = quantities[name][:].flat |
|---|
| 450 | self.quantities_range[name] = [min(q), max(q)] |
|---|
| 451 | |
|---|
| 452 | self.quantity_names = quantity_names |
|---|
| 453 | #self.quantities = quantities #Takes too much memory |
|---|
| 454 | self.vertex_coordinates = vertex_coordinates |
|---|
| 455 | self.interpolation_points = interpolation_points |
|---|
| 456 | self.time = time[:] # Time assumed to be relative to starttime |
|---|
| 457 | self.index = 0 # Initial time index |
|---|
| 458 | self.precomputed_values = {} |
|---|
| 459 | |
|---|
| 460 | |
|---|
| 461 | |
|---|
| 462 | |
|---|
| 463 | |
|---|
| 464 | |
|---|
| 465 | #Precomputed spatial interpolation if requested |
|---|
| 466 | if interpolation_points is not None: |
|---|
| 467 | if self.spatial is False: |
|---|
| 468 | raise 'Triangles and vertex_coordinates must be specified' |
|---|
| 469 | |
|---|
| 470 | try: |
|---|
| 471 | self.interpolation_points = ensure_numeric(interpolation_points) |
|---|
| 472 | except: |
|---|
| 473 | msg = 'Interpolation points must be an N x 2 Numeric array '+\ |
|---|
| 474 | 'or a list of points\n' |
|---|
| 475 | msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\ |
|---|
| 476 | '...') |
|---|
| 477 | raise msg |
|---|
| 478 | |
|---|
| 479 | |
|---|
| 480 | m = len(self.interpolation_points) |
|---|
| 481 | p = len(self.time) |
|---|
| 482 | |
|---|
| 483 | for name in quantity_names: |
|---|
| 484 | self.precomputed_values[name] = zeros((p, m), Float) |
|---|
| 485 | |
|---|
| 486 | #Build interpolator |
|---|
| 487 | interpol = Interpolate(vertex_coordinates, |
|---|
| 488 | triangles, |
|---|
| 489 | #point_coordinates = \ |
|---|
| 490 | #self.interpolation_points, |
|---|
| 491 | #alpha = 0, |
|---|
| 492 | verbose = verbose) |
|---|
| 493 | |
|---|
| 494 | if verbose: print 'Interpolate' |
|---|
| 495 | for i, t in enumerate(self.time): |
|---|
| 496 | #Interpolate quantities at this timestep |
|---|
| 497 | if verbose and i%((p+10)/10)==0: |
|---|
| 498 | print ' time step %d of %d' %(i, p) |
|---|
| 499 | |
|---|
| 500 | for name in quantity_names: |
|---|
| 501 | if len(quantities[name].shape) == 2: |
|---|
| 502 | result = interpol.interpolate(quantities[name][i,:], |
|---|
| 503 | point_coordinates = \ |
|---|
| 504 | self.interpolation_points) |
|---|
| 505 | else: |
|---|
| 506 | #Assume no time dependency |
|---|
| 507 | result = interpol.interpolate(quantities[name][:], |
|---|
| 508 | point_coordinates = \ |
|---|
| 509 | self.interpolation_points) |
|---|
| 510 | |
|---|
| 511 | self.precomputed_values[name][i, :] = result |
|---|
| 512 | |
|---|
| 513 | #Report |
|---|
| 514 | if verbose: |
|---|
| 515 | print self.statistics() |
|---|
| 516 | #self.print_statistics() |
|---|
| 517 | |
|---|
| 518 | else: |
|---|
| 519 | #Store quantitites as is |
|---|
| 520 | for name in quantity_names: |
|---|
| 521 | self.precomputed_values[name] = quantities[name] |
|---|
| 522 | |
|---|
| 523 | |
|---|
| 524 | def __repr__(self): |
|---|
| 525 | #return 'Interpolation function (spatio-temporal)' |
|---|
| 526 | return self.statistics() |
|---|
| 527 | |
|---|
| 528 | |
|---|
| 529 | def __call__(self, t, point_id = None, x = None, y = None): |
|---|
| 530 | """Evaluate f(t), f(t, point_id) or f(t, x, y) |
|---|
| 531 | |
|---|
| 532 | Inputs: |
|---|
| 533 | t: time - Model time. Must lie within existing timesteps |
|---|
| 534 | point_id: index of one of the preprocessed points. |
|---|
| 535 | x, y: Overrides location, point_id ignored |
|---|
| 536 | |
|---|
| 537 | If spatial info is present and all of x,y,point_id |
|---|
| 538 | are None an exception is raised |
|---|
| 539 | |
|---|
| 540 | If no spatial info is present, point_id and x,y arguments are ignored |
|---|
| 541 | making f a function of time only. |
|---|
| 542 | |
|---|
| 543 | |
|---|
| 544 | FIXME: point_id could also be a slice |
|---|
| 545 | FIXME: What if x and y are vectors? |
|---|
| 546 | FIXME: What about f(x,y) without t? |
|---|
| 547 | """ |
|---|
| 548 | |
|---|
| 549 | from math import pi, cos, sin, sqrt |
|---|
| 550 | from Numeric import zeros, Float |
|---|
| 551 | from util import mean |
|---|
| 552 | |
|---|
| 553 | if self.spatial is True: |
|---|
| 554 | if point_id is None: |
|---|
| 555 | if x is None or y is None: |
|---|
| 556 | msg = 'Either point_id or x and y must be specified' |
|---|
| 557 | raise Exception(msg) |
|---|
| 558 | else: |
|---|
| 559 | if self.interpolation_points is None: |
|---|
| 560 | msg = 'Interpolation_function must be instantiated ' +\ |
|---|
| 561 | 'with a list of interpolation points before parameter ' +\ |
|---|
| 562 | 'point_id can be used' |
|---|
| 563 | raise Exception(msg) |
|---|
| 564 | |
|---|
| 565 | msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1]) |
|---|
| 566 | msg += ' does not match model time: %.16f\n' %t |
|---|
| 567 | if t < self.time[0]: raise Exception(msg) |
|---|
| 568 | if t > self.time[-1]: raise Exception(msg) |
|---|
| 569 | |
|---|
| 570 | oldindex = self.index #Time index |
|---|
| 571 | |
|---|
| 572 | #Find current time slot |
|---|
| 573 | while t > self.time[self.index]: self.index += 1 |
|---|
| 574 | while t < self.time[self.index]: self.index -= 1 |
|---|
| 575 | |
|---|
| 576 | if t == self.time[self.index]: |
|---|
| 577 | #Protect against case where t == T[-1] (last time) |
|---|
| 578 | # - also works in general when t == T[i] |
|---|
| 579 | ratio = 0 |
|---|
| 580 | else: |
|---|
| 581 | #t is now between index and index+1 |
|---|
| 582 | ratio = (t - self.time[self.index])/\ |
|---|
| 583 | (self.time[self.index+1] - self.time[self.index]) |
|---|
| 584 | |
|---|
| 585 | #Compute interpolated values |
|---|
| 586 | q = zeros(len(self.quantity_names), Float) |
|---|
| 587 | #print "self.precomputed_values", self.precomputed_values |
|---|
| 588 | for i, name in enumerate(self.quantity_names): |
|---|
| 589 | Q = self.precomputed_values[name] |
|---|
| 590 | |
|---|
| 591 | if self.spatial is False: |
|---|
| 592 | #If there is no spatial info |
|---|
| 593 | assert len(Q.shape) == 1 |
|---|
| 594 | |
|---|
| 595 | Q0 = Q[self.index] |
|---|
| 596 | if ratio > 0: Q1 = Q[self.index+1] |
|---|
| 597 | |
|---|
| 598 | else: |
|---|
| 599 | if x is not None and y is not None: |
|---|
| 600 | #Interpolate to x, y |
|---|
| 601 | |
|---|
| 602 | raise 'x,y interpolation not yet implemented' |
|---|
| 603 | else: |
|---|
| 604 | #Use precomputed point |
|---|
| 605 | Q0 = Q[self.index, point_id] |
|---|
| 606 | if ratio > 0: |
|---|
| 607 | Q1 = Q[self.index+1, point_id] |
|---|
| 608 | |
|---|
| 609 | #Linear temporal interpolation |
|---|
| 610 | if ratio > 0: |
|---|
| 611 | if Q0 == NAN and Q1 == NAN: |
|---|
| 612 | q[i] = Q0 |
|---|
| 613 | else: |
|---|
| 614 | q[i] = Q0 + ratio*(Q1 - Q0) |
|---|
| 615 | else: |
|---|
| 616 | q[i] = Q0 |
|---|
| 617 | |
|---|
| 618 | |
|---|
| 619 | #Return vector of interpolated values |
|---|
| 620 | #if len(q) == 1: |
|---|
| 621 | # return q[0] |
|---|
| 622 | #else: |
|---|
| 623 | # return q |
|---|
| 624 | |
|---|
| 625 | |
|---|
| 626 | #Return vector of interpolated values |
|---|
| 627 | #FIXME: |
|---|
| 628 | if self.spatial is True: |
|---|
| 629 | return q |
|---|
| 630 | else: |
|---|
| 631 | #Replicate q according to x and y |
|---|
| 632 | #This is e.g used for Wind_stress |
|---|
| 633 | if x is None or y is None: |
|---|
| 634 | return q |
|---|
| 635 | else: |
|---|
| 636 | try: |
|---|
| 637 | N = len(x) |
|---|
| 638 | except: |
|---|
| 639 | return q |
|---|
| 640 | else: |
|---|
| 641 | from Numeric import ones, Float |
|---|
| 642 | #x is a vector - Create one constant column for each value |
|---|
| 643 | N = len(x) |
|---|
| 644 | assert len(y) == N, 'x and y must have same length' |
|---|
| 645 | res = [] |
|---|
| 646 | for col in q: |
|---|
| 647 | res.append(col*ones(N, Float)) |
|---|
| 648 | |
|---|
| 649 | return res |
|---|
| 650 | |
|---|
| 651 | |
|---|
| 652 | def get_time(self): |
|---|
| 653 | """Return model time as a vector of timesteps |
|---|
| 654 | """ |
|---|
| 655 | return self.time |
|---|
| 656 | |
|---|
| 657 | def statistics(self): |
|---|
| 658 | """Output statistics about interpolation_function |
|---|
| 659 | """ |
|---|
| 660 | |
|---|
| 661 | vertex_coordinates = self.vertex_coordinates |
|---|
| 662 | interpolation_points = self.interpolation_points |
|---|
| 663 | quantity_names = self.quantity_names |
|---|
| 664 | #quantities = self.quantities |
|---|
| 665 | precomputed_values = self.precomputed_values |
|---|
| 666 | |
|---|
| 667 | x = vertex_coordinates[:,0] |
|---|
| 668 | y = vertex_coordinates[:,1] |
|---|
| 669 | |
|---|
| 670 | str = '------------------------------------------------\n' |
|---|
| 671 | str += 'Interpolation_function (spatio-temporal) statistics:\n' |
|---|
| 672 | str += ' Extent:\n' |
|---|
| 673 | str += ' x in [%f, %f], len(x) == %d\n'\ |
|---|
| 674 | %(min(x), max(x), len(x)) |
|---|
| 675 | str += ' y in [%f, %f], len(y) == %d\n'\ |
|---|
| 676 | %(min(y), max(y), len(y)) |
|---|
| 677 | str += ' t in [%f, %f], len(t) == %d\n'\ |
|---|
| 678 | %(min(self.time), max(self.time), len(self.time)) |
|---|
| 679 | str += ' Quantities:\n' |
|---|
| 680 | for name in quantity_names: |
|---|
| 681 | minq, maxq = self.quantities_range[name] |
|---|
| 682 | str += ' %s in [%f, %f]\n' %(name, minq, maxq) |
|---|
| 683 | #q = quantities[name][:].flat |
|---|
| 684 | #str += ' %s in [%f, %f]\n' %(name, min(q), max(q)) |
|---|
| 685 | |
|---|
| 686 | if interpolation_points is not None: |
|---|
| 687 | str += ' Interpolation points (xi, eta):'\ |
|---|
| 688 | ' number of points == %d\n' %interpolation_points.shape[0] |
|---|
| 689 | str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]), |
|---|
| 690 | max(interpolation_points[:,0])) |
|---|
| 691 | str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]), |
|---|
| 692 | max(interpolation_points[:,1])) |
|---|
| 693 | str += ' Interpolated quantities (over all timesteps):\n' |
|---|
| 694 | |
|---|
| 695 | for name in quantity_names: |
|---|
| 696 | q = precomputed_values[name][:].flat |
|---|
| 697 | str += ' %s at interpolation points in [%f, %f]\n'\ |
|---|
| 698 | %(name, min(q), max(q)) |
|---|
| 699 | str += '------------------------------------------------\n' |
|---|
| 700 | |
|---|
| 701 | return str |
|---|
| 702 | |
|---|
| 703 | def interpolate_sww(sww_file, time, interpolation_points, |
|---|
| 704 | quantity_names = None, verbose = False): |
|---|
| 705 | """ |
|---|
| 706 | obsolete. |
|---|
| 707 | use file_function in utils |
|---|
| 708 | """ |
|---|
| 709 | #open sww file |
|---|
| 710 | x, y, volumes, time, quantities = read_sww(sww_file) |
|---|
| 711 | print "x",x |
|---|
| 712 | print "y",y |
|---|
| 713 | |
|---|
| 714 | print "time", time |
|---|
| 715 | print "quantities", quantities |
|---|
| 716 | |
|---|
| 717 | #Add the x and y together |
|---|
| 718 | vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1) |
|---|
| 719 | |
|---|
| 720 | #Will return the quantity values at the specified times and locations |
|---|
| 721 | interp = Interpolation_interface( |
|---|
| 722 | time, |
|---|
| 723 | quantities, |
|---|
| 724 | quantity_names = quantity_names, |
|---|
| 725 | vertex_coordinates = vertex_coordinates, |
|---|
| 726 | triangles = volumes, |
|---|
| 727 | interpolation_points = interpolation_points, |
|---|
| 728 | verbose = verbose) |
|---|
| 729 | |
|---|
| 730 | |
|---|
| 731 | def read_sww(file_name): |
|---|
| 732 | """ |
|---|
| 733 | obsolete - Nothing should be calling this |
|---|
| 734 | |
|---|
| 735 | Read in an sww file. |
|---|
| 736 | |
|---|
| 737 | Input; |
|---|
| 738 | file_name - the sww file |
|---|
| 739 | |
|---|
| 740 | Output; |
|---|
| 741 | x - Vector of x values |
|---|
| 742 | y - Vector of y values |
|---|
| 743 | z - Vector of bed elevation |
|---|
| 744 | volumes - Array. Each row has 3 values, representing |
|---|
| 745 | the vertices that define the volume |
|---|
| 746 | time - Vector of the times where there is stage information |
|---|
| 747 | stage - array with respect to time and vertices (x,y) |
|---|
| 748 | """ |
|---|
| 749 | |
|---|
| 750 | #FIXME Have this reader as part of data_manager? |
|---|
| 751 | |
|---|
| 752 | from Scientific.IO.NetCDF import NetCDFFile |
|---|
| 753 | import tempfile |
|---|
| 754 | import sys |
|---|
| 755 | import os |
|---|
| 756 | |
|---|
| 757 | #Check contents |
|---|
| 758 | #Get NetCDF |
|---|
| 759 | |
|---|
| 760 | # see if the file is there. Throw a QUIET IO error if it isn't |
|---|
| 761 | # I don't think I could implement the above |
|---|
| 762 | |
|---|
| 763 | #throws prints to screen if file not present |
|---|
| 764 | junk = tempfile.mktemp(".txt") |
|---|
| 765 | fd = open(junk,'w') |
|---|
| 766 | stdout = sys.stdout |
|---|
| 767 | sys.stdout = fd |
|---|
| 768 | fid = NetCDFFile(file_name, 'r') |
|---|
| 769 | sys.stdout = stdout |
|---|
| 770 | fd.close() |
|---|
| 771 | #clean up |
|---|
| 772 | os.remove(junk) |
|---|
| 773 | |
|---|
| 774 | # Get the variables |
|---|
| 775 | x = fid.variables['x'][:] |
|---|
| 776 | y = fid.variables['y'][:] |
|---|
| 777 | volumes = fid.variables['volumes'][:] |
|---|
| 778 | time = fid.variables['time'][:] |
|---|
| 779 | |
|---|
| 780 | keys = fid.variables.keys() |
|---|
| 781 | keys.remove("x") |
|---|
| 782 | keys.remove("y") |
|---|
| 783 | keys.remove("volumes") |
|---|
| 784 | keys.remove("time") |
|---|
| 785 | #Turn NetCDF objects into Numeric arrays |
|---|
| 786 | quantities = {} |
|---|
| 787 | for name in keys: |
|---|
| 788 | quantities[name] = fid.variables[name][:] |
|---|
| 789 | |
|---|
| 790 | fid.close() |
|---|
| 791 | return x, y, volumes, time, quantities |
|---|
| 792 | |
|---|
| 793 | #------------------------------------------------------------- |
|---|
| 794 | if __name__ == "__main__": |
|---|
| 795 | names = ["x","y"] |
|---|
| 796 | someiterable = [[1,2],[3,4]] |
|---|
| 797 | csvwriter = writer(file("some.csv", "wb")) |
|---|
| 798 | csvwriter.writerow(names) |
|---|
| 799 | for row in someiterable: |
|---|
| 800 | csvwriter.writerow(row) |
|---|
