[2879] | 1 | """Least squares interpolation. |
---|
| 2 | |
---|
| 3 | Implements a least-squares interpolation. |
---|
[3409] | 4 | Putting mesh data onto points. |
---|
[2879] | 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 | |
---|
[3319] | 13 | * The public interface to Interpolate |
---|
[2879] | 14 | __init__ |
---|
| 15 | interpolate |
---|
| 16 | interpolate_block |
---|
| 17 | |
---|
| 18 | """ |
---|
| 19 | |
---|
| 20 | import time |
---|
| 21 | import os |
---|
| 22 | from warnings import warn |
---|
[3411] | 23 | from math import sqrt |
---|
| 24 | from csv import writer, DictWriter |
---|
[2879] | 25 | |
---|
| 26 | from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \ |
---|
| 27 | ArrayType, allclose, take, NewAxis, arange |
---|
| 28 | |
---|
[3514] | 29 | from anuga.caching.caching import cache |
---|
[3560] | 30 | from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh |
---|
[3514] | 31 | from anuga.utilities.sparse import Sparse, Sparse_CSR |
---|
| 32 | from anuga.utilities.cg_solve import conjugate_gradient, VectorShapeError |
---|
| 33 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
| 34 | from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN |
---|
| 35 | from anuga.utilities.polygon import in_and_outside_polygon |
---|
[3768] | 36 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
| 37 | from anuga.geospatial_data.geospatial_data import ensure_absolute |
---|
[3514] | 38 | from anuga.fit_interpolate.search_functions import search_tree_of_vertices |
---|
| 39 | from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate |
---|
[3560] | 40 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
[2879] | 41 | |
---|
| 42 | |
---|
| 43 | class Interpolate (FitInterpolate): |
---|
| 44 | |
---|
| 45 | def __init__(self, |
---|
| 46 | vertex_coordinates, |
---|
| 47 | triangles, |
---|
| 48 | mesh_origin=None, |
---|
| 49 | verbose=False, |
---|
[4614] | 50 | max_vertices_per_cell=None): |
---|
[2879] | 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 | """ |
---|
[3689] | 78 | |
---|
| 79 | # FIXME (Ole): Need an input check |
---|
[2879] | 80 | |
---|
| 81 | # Initialise variabels |
---|
| 82 | self._A_can_be_reused = False |
---|
| 83 | self._point_coordinates = None |
---|
| 84 | |
---|
| 85 | FitInterpolate.__init__(self, |
---|
[4779] | 86 | vertex_coordinates=vertex_coordinates, |
---|
| 87 | triangles=triangles, |
---|
| 88 | mesh_origin=mesh_origin, |
---|
| 89 | verbose=verbose, |
---|
| 90 | max_vertices_per_cell=max_vertices_per_cell) |
---|
[2879] | 91 | |
---|
[5361] | 92 | def interpolate_polyline(self, |
---|
| 93 | f, |
---|
| 94 | vertex_coordinates, |
---|
[5372] | 95 | gauge_neighbour_id, |
---|
[5361] | 96 | point_coordinates=None, |
---|
| 97 | verbose=False): |
---|
[5403] | 98 | """Interpolate linearly between values f on nodes (vertex coordinates) of a polyline to midpoints of triangles |
---|
| 99 | of boundary. |
---|
| 100 | |
---|
| 101 | f is the data on the polyline nodes. |
---|
| 102 | |
---|
| 103 | The mesh values representing a smooth surface are |
---|
| 104 | assumed to be specified in f. |
---|
| 105 | |
---|
| 106 | Inputs: |
---|
| 107 | f: Vector or array of data at the polyline nodes. |
---|
| 108 | If f is an array, interpolation will be done for each column as |
---|
| 109 | per underlying matrix-matrix multiplication |
---|
| 110 | point_coordinates: Interpolate polyline data to these positions. |
---|
| 111 | List of coordinate pairs [x, y] of |
---|
| 112 | data points or an nx2 Numeric array or a Geospatial_data object |
---|
| 113 | |
---|
| 114 | Output: |
---|
| 115 | Interpolated values at inputted points (z). |
---|
| 116 | """ |
---|
| 117 | |
---|
| 118 | |
---|
| 119 | # FIXME: There is an option of passing a tolerance into this |
---|
| 120 | |
---|
[5361] | 121 | if isinstance(point_coordinates, Geospatial_data): |
---|
| 122 | point_coordinates = point_coordinates.get_data_points( \ |
---|
| 123 | absolute = True) |
---|
| 124 | |
---|
[5403] | 125 | from utilities.polygon import point_on_line |
---|
[5361] | 126 | from Numeric import ones |
---|
| 127 | z=ones(len(point_coordinates),Float) |
---|
| 128 | |
---|
| 129 | msg='point coordinates are not given (interpolate.py)' |
---|
| 130 | assert point_coordinates is not None, msg |
---|
| 131 | msg='function value must be specified at every interpolation node' |
---|
[5403] | 132 | assert f.shape[0]==vertex_coordinates.shape[0], msg |
---|
[5361] | 133 | msg='Must define function value at one or more nodes' |
---|
[5403] | 134 | assert f.shape[0]>0, msg |
---|
[5361] | 135 | |
---|
| 136 | n=f.shape[0] |
---|
| 137 | if n==1: |
---|
| 138 | z=f*z |
---|
[5403] | 139 | msg = 'Polyline contained only one point. I need more. ', str(f) |
---|
| 140 | raise Exception, msg |
---|
| 141 | |
---|
| 142 | # FIXME (John): add unit test for only One vertex point. Exception should be thrown. |
---|
| 143 | |
---|
| 144 | |
---|
[5361] | 145 | elif n>1: |
---|
| 146 | for i in range(len(point_coordinates)): |
---|
| 147 | found = False |
---|
[5372] | 148 | for j in range(n): |
---|
| 149 | if gauge_neighbour_id[j]>=0: |
---|
[5403] | 150 | if point_on_line(point_coordinates[i], |
---|
| 151 | [vertex_coordinates[j], vertex_coordinates[gauge_neighbour_id[j]]], |
---|
| 152 | rtol=1.0e-6): |
---|
[5372] | 153 | found=True |
---|
| 154 | x0=vertex_coordinates[j][0] |
---|
| 155 | y0=vertex_coordinates[j][1] |
---|
| 156 | x1=vertex_coordinates[gauge_neighbour_id[j]][0] |
---|
| 157 | y1=vertex_coordinates[gauge_neighbour_id[j]][1] |
---|
| 158 | x2=point_coordinates[i][0] |
---|
| 159 | y2=point_coordinates[i][1] |
---|
| 160 | |
---|
| 161 | segment_len=sqrt((x1-x0)**2+(y1-y0)**2) |
---|
| 162 | dist=sqrt((x2-x0)**2+(y2-y0)**2) |
---|
| 163 | z[i]=(f[gauge_neighbour_id[j]]-f[j])/segment_len*dist+f[j] |
---|
| 164 | #print 'element found on segment' |
---|
| 165 | break |
---|
| 166 | |
---|
[5361] | 167 | if not found: |
---|
| 168 | z[i]=0.0 |
---|
| 169 | #print 'point not on urs boundary' |
---|
| 170 | return z |
---|
| 171 | |
---|
[3689] | 172 | # FIXME: What is a good start_blocking_len value? |
---|
[3896] | 173 | def interpolate(self, |
---|
| 174 | f, |
---|
| 175 | point_coordinates=None, |
---|
| 176 | start_blocking_len=500000, |
---|
| 177 | verbose=False): |
---|
[2879] | 178 | """Interpolate mesh data f to determine values, z, at points. |
---|
| 179 | |
---|
| 180 | f is the data on the mesh vertices. |
---|
| 181 | |
---|
| 182 | The mesh values representing a smooth surface are |
---|
| 183 | assumed to be specified in f. |
---|
| 184 | |
---|
| 185 | Inputs: |
---|
| 186 | f: Vector or array of data at the mesh vertices. |
---|
| 187 | If f is an array, interpolation will be done for each column as |
---|
| 188 | per underlying matrix-matrix multiplication |
---|
| 189 | point_coordinates: Interpolate mesh data to these positions. |
---|
| 190 | List of coordinate pairs [x, y] of |
---|
| 191 | data points or an nx2 Numeric array or a Geospatial_data object |
---|
| 192 | |
---|
| 193 | If point_coordinates is absent, the points inputted last time |
---|
| 194 | this method was called are used, if possible. |
---|
| 195 | start_blocking_len: If the # of points is more or greater than this, |
---|
| 196 | start blocking |
---|
| 197 | |
---|
| 198 | Output: |
---|
| 199 | Interpolated values at inputted points (z). |
---|
| 200 | """ |
---|
[3689] | 201 | |
---|
[3768] | 202 | # FIXME (Ole): Why is the interpolation matrix rebuilt everytime the |
---|
| 203 | # method is called even if interpolation points are unchanged. |
---|
[5222] | 204 | |
---|
[5321] | 205 | #print "point_coordinates interpolate.interpolate", point_coordinates |
---|
| 206 | if verbose: print 'Build intepolation object' |
---|
[3689] | 207 | if isinstance(point_coordinates, Geospatial_data): |
---|
[2879] | 208 | point_coordinates = point_coordinates.get_data_points( \ |
---|
| 209 | absolute = True) |
---|
[3689] | 210 | |
---|
[2879] | 211 | # Can I interpolate, based on previous point_coordinates? |
---|
| 212 | if point_coordinates is None: |
---|
| 213 | if self._A_can_be_reused is True and \ |
---|
[3689] | 214 | len(self._point_coordinates) < start_blocking_len: |
---|
[2879] | 215 | z = self._get_point_data_z(f, |
---|
| 216 | verbose=verbose) |
---|
| 217 | elif self._point_coordinates is not None: |
---|
| 218 | # if verbose, give warning |
---|
| 219 | if verbose: |
---|
| 220 | print 'WARNING: Recalculating A matrix, due to blocking.' |
---|
| 221 | point_coordinates = self._point_coordinates |
---|
| 222 | else: |
---|
| 223 | #There are no good point_coordinates. import sys; sys.exit() |
---|
| 224 | msg = 'ERROR (interpolate.py): No point_coordinates inputted' |
---|
| 225 | raise Exception(msg) |
---|
[5321] | 226 | |
---|
| 227 | if point_coordinates is not None: |
---|
[2879] | 228 | self._point_coordinates = point_coordinates |
---|
| 229 | if len(point_coordinates) < start_blocking_len or \ |
---|
| 230 | start_blocking_len == 0: |
---|
| 231 | self._A_can_be_reused = True |
---|
| 232 | z = self.interpolate_block(f, point_coordinates, |
---|
| 233 | verbose=verbose) |
---|
| 234 | else: |
---|
[3689] | 235 | #print 'BLOCKING' |
---|
[2879] | 236 | #Handle blocking |
---|
| 237 | self._A_can_be_reused = False |
---|
[3689] | 238 | start = 0 |
---|
[2879] | 239 | # creating a dummy array to concatenate to. |
---|
| 240 | |
---|
| 241 | f = ensure_numeric(f, Float) |
---|
| 242 | #print "f.shape",f.shape |
---|
| 243 | if len(f.shape) > 1: |
---|
| 244 | z = zeros((0,f.shape[1])) |
---|
| 245 | else: |
---|
| 246 | z = zeros((0,)) |
---|
| 247 | |
---|
[3689] | 248 | for end in range(start_blocking_len, |
---|
| 249 | len(point_coordinates), |
---|
| 250 | start_blocking_len): |
---|
[3768] | 251 | |
---|
[2879] | 252 | t = self.interpolate_block(f, point_coordinates[start:end], |
---|
| 253 | verbose=verbose) |
---|
| 254 | #print "t", t |
---|
| 255 | #print "z", z |
---|
| 256 | z = concatenate((z,t)) |
---|
| 257 | start = end |
---|
[3768] | 258 | |
---|
[2879] | 259 | end = len(point_coordinates) |
---|
| 260 | t = self.interpolate_block(f, point_coordinates[start:end], |
---|
| 261 | verbose=verbose) |
---|
| 262 | z = concatenate((z,t)) |
---|
| 263 | return z |
---|
[5222] | 264 | |
---|
[2879] | 265 | |
---|
[5222] | 266 | def interpolate_block(self, f, point_coordinates, verbose=False): |
---|
[2879] | 267 | """ |
---|
| 268 | Call this if you want to control the blocking or make sure blocking |
---|
| 269 | doesn't occur. |
---|
| 270 | |
---|
| 271 | Return the point data, z. |
---|
| 272 | |
---|
| 273 | See interpolate for doc info. |
---|
| 274 | """ |
---|
| 275 | if isinstance(point_coordinates,Geospatial_data): |
---|
[5222] | 276 | point_coordinates = point_coordinates.get_data_points(\ |
---|
| 277 | absolute=True) |
---|
| 278 | |
---|
| 279 | # Convert lists to Numeric arrays if necessary |
---|
| 280 | point_coordinates = ensure_numeric(point_coordinates, Float) |
---|
| 281 | f = ensure_numeric(f, Float) |
---|
| 282 | |
---|
| 283 | self._A = self._build_interpolation_matrix_A(point_coordinates, |
---|
| 284 | verbose=verbose) |
---|
| 285 | |
---|
| 286 | |
---|
| 287 | # Check that input dimensions are compatible |
---|
| 288 | msg = 'Two colums must be specified in point coordinates. I got shape=%s'\ |
---|
| 289 | %(str(point_coordinates.shape)) |
---|
| 290 | assert point_coordinates.shape[1] == 2, msg |
---|
| 291 | |
---|
| 292 | msg = 'The number of rows in matrix A must be the same as the number of points supplied.' |
---|
| 293 | msg += ' I got %d points and %d matrix rows.'\ |
---|
| 294 | %(point_coordinates.shape[0], self._A.shape[0]) |
---|
| 295 | assert point_coordinates.shape[0] == self._A.shape[0], msg |
---|
| 296 | |
---|
| 297 | msg = 'The number of columns in matrix A must be the same as the number of mesh vertices.' |
---|
| 298 | msg += ' I got %d vertices and %d matrix columns.'\ |
---|
| 299 | %(f.shape[0], self._A.shape[1]) |
---|
| 300 | assert self._A.shape[1] == f.shape[0], msg |
---|
| 301 | |
---|
| 302 | # Compute Matrix vector product and return |
---|
[2879] | 303 | return self._get_point_data_z(f) |
---|
[5222] | 304 | |
---|
[2879] | 305 | |
---|
| 306 | def _get_point_data_z(self, f, verbose=False): |
---|
| 307 | """ |
---|
| 308 | Return the point data, z. |
---|
| 309 | |
---|
| 310 | Precondition, |
---|
| 311 | The _A matrix has been created |
---|
| 312 | """ |
---|
[5222] | 313 | |
---|
[2879] | 314 | z = self._A * f |
---|
| 315 | # Taking into account points outside the mesh. |
---|
| 316 | #print "self.outside_poly_indices", self.outside_poly_indices |
---|
| 317 | #print "self.inside_poly_indices", self.inside_poly_indices |
---|
| 318 | #print "z", z |
---|
| 319 | for i in self.outside_poly_indices: |
---|
[3452] | 320 | z[i] = NAN |
---|
[2879] | 321 | return z |
---|
| 322 | |
---|
| 323 | def _build_interpolation_matrix_A(self, |
---|
[3689] | 324 | point_coordinates, |
---|
[3896] | 325 | verbose=False): |
---|
[2879] | 326 | """Build n x m interpolation matrix, where |
---|
| 327 | n is the number of data points and |
---|
| 328 | m is the number of basis functions phi_k (one per vertex) |
---|
| 329 | |
---|
| 330 | This algorithm uses a quad tree data structure for fast binning |
---|
| 331 | of data points |
---|
| 332 | origin is a 3-tuple consisting of UTM zone, easting and northing. |
---|
| 333 | If specified coordinates are assumed to be relative to this origin. |
---|
| 334 | |
---|
| 335 | This one will override any data_origin that may be specified in |
---|
| 336 | instance interpolation |
---|
| 337 | |
---|
| 338 | Preconditions |
---|
| 339 | Point_coordindates and mesh vertices have the same origin. |
---|
| 340 | """ |
---|
| 341 | |
---|
[5321] | 342 | if verbose: print 'Building interpolation matrix' |
---|
[5222] | 343 | |
---|
| 344 | # Convert point_coordinates to Numeric arrays, in case it was a list. |
---|
[2879] | 345 | point_coordinates = ensure_numeric(point_coordinates, Float) |
---|
[4648] | 346 | |
---|
[5222] | 347 | |
---|
[2879] | 348 | if verbose: print 'Getting indices inside mesh boundary' |
---|
| 349 | self.inside_poly_indices, self.outside_poly_indices = \ |
---|
[3686] | 350 | in_and_outside_polygon(point_coordinates, |
---|
| 351 | self.mesh.get_boundary_polygon(), |
---|
| 352 | closed = True, verbose = verbose) |
---|
[4648] | 353 | |
---|
[2879] | 354 | #Build n x m interpolation matrix |
---|
[4589] | 355 | if verbose and len(self.outside_poly_indices) > 0: |
---|
[3018] | 356 | print '\n WARNING: Points outside mesh boundary. \n' |
---|
[3019] | 357 | # Since you can block, throw a warning, not an error. |
---|
| 358 | if verbose and 0 == len(self.inside_poly_indices): |
---|
| 359 | print '\n WARNING: No points within the mesh! \n' |
---|
| 360 | |
---|
[3945] | 361 | m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex) |
---|
| 362 | n = point_coordinates.shape[0] # Nbr of data points |
---|
[2879] | 363 | |
---|
| 364 | if verbose: print 'Number of datapoints: %d' %n |
---|
| 365 | if verbose: print 'Number of basis functions: %d' %m |
---|
| 366 | |
---|
| 367 | A = Sparse(n,m) |
---|
| 368 | |
---|
| 369 | n = len(self.inside_poly_indices) |
---|
| 370 | #Compute matrix elements for points inside the mesh |
---|
[4589] | 371 | if verbose: print 'Building interpolation matrix from %d points' %n |
---|
| 372 | for d, i in enumerate(self.inside_poly_indices): |
---|
| 373 | # For each data_coordinate point |
---|
| 374 | if verbose and d%((n+10)/10)==0: print 'Doing %d of %d' %(d, n) |
---|
[2879] | 375 | x = point_coordinates[i] |
---|
| 376 | element_found, sigma0, sigma1, sigma2, k = \ |
---|
| 377 | search_tree_of_vertices(self.root, self.mesh, x) |
---|
[4589] | 378 | |
---|
| 379 | # Update interpolation matrix A if necessary |
---|
[2879] | 380 | if element_found is True: |
---|
[4589] | 381 | # Assign values to matrix A |
---|
[2879] | 382 | |
---|
[4589] | 383 | j0 = self.mesh.triangles[k,0] # Global vertex id for sigma0 |
---|
| 384 | j1 = self.mesh.triangles[k,1] # Global vertex id for sigma1 |
---|
| 385 | j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2 |
---|
[2879] | 386 | |
---|
| 387 | sigmas = {j0:sigma0, j1:sigma1, j2:sigma2} |
---|
| 388 | js = [j0,j1,j2] |
---|
| 389 | |
---|
| 390 | for j in js: |
---|
| 391 | A[i,j] = sigmas[j] |
---|
| 392 | else: |
---|
| 393 | msg = 'Could not find triangle for point', x |
---|
| 394 | raise Exception(msg) |
---|
| 395 | return A |
---|
| 396 | |
---|
[4651] | 397 | def benchmark_interpolate(vertices, |
---|
| 398 | vertex_attributes, |
---|
| 399 | triangles, points, |
---|
| 400 | max_points_per_cell=None, |
---|
[4839] | 401 | start_blocking_len=500000, |
---|
| 402 | mesh_origin=None): |
---|
[4651] | 403 | """ |
---|
[4839] | 404 | points: Interpolate mesh data to these positions. |
---|
| 405 | List of coordinate pairs [x, y] of |
---|
| 406 | data points or an nx2 Numeric array or a Geospatial_data object |
---|
| 407 | |
---|
[4651] | 408 | No test for this yet. |
---|
| 409 | Note, this has no time the input data has no time dimension. Which is |
---|
| 410 | different from most of the data we interpolate, eg sww info. |
---|
| 411 | |
---|
| 412 | Output: |
---|
| 413 | Interpolated values at inputted points. |
---|
| 414 | """ |
---|
| 415 | interp = Interpolate(vertices, |
---|
| 416 | triangles, |
---|
[4839] | 417 | max_vertices_per_cell=max_points_per_cell, |
---|
| 418 | mesh_origin=mesh_origin) |
---|
[4651] | 419 | |
---|
| 420 | calc = interp.interpolate(vertex_attributes |
---|
| 421 | ,points |
---|
| 422 | ,start_blocking_len=start_blocking_len) |
---|
[4839] | 423 | #print "calc", calc |
---|
| 424 | |
---|
[3419] | 425 | def interpolate_sww2csv(sww_file, |
---|
| 426 | points, |
---|
[3411] | 427 | depth_file, |
---|
[3456] | 428 | velocity_x_file, |
---|
| 429 | velocity_y_file, |
---|
[5196] | 430 | stage_file=None, |
---|
[5496] | 431 | froude_file=None, |
---|
[3417] | 432 | #quantities = ['depth', 'velocity'], |
---|
[5447] | 433 | time_thinning=1, |
---|
[3417] | 434 | verbose=True, |
---|
| 435 | use_cache = True): |
---|
[3411] | 436 | """ |
---|
| 437 | Interpolate the quantities at a given set of locations, given |
---|
| 438 | an sww file. |
---|
[5454] | 439 | The results are written to csv files. |
---|
| 440 | |
---|
| 441 | sww_file is the input sww file. |
---|
| 442 | points is a list of the 'gauges' x,y location. |
---|
| 443 | depth_file is the name of the output depth file |
---|
| 444 | velocity_x_file is the name of the output x velocity file. |
---|
| 445 | velocity_y_file is the name of the output y velocity file. |
---|
| 446 | stage_file is the name of the output stage file. |
---|
| 447 | |
---|
| 448 | In the csv files columns represents the gauges and each row is a |
---|
| 449 | time slice. |
---|
[5447] | 450 | |
---|
[5454] | 451 | |
---|
[5447] | 452 | Time_thinning_number controls how many timesteps to use. Only |
---|
| 453 | timesteps with index%time_thinning_number == 0 will used, or |
---|
| 454 | in other words a value of 3, say, will cause the algorithm to |
---|
| 455 | use every third time step. |
---|
[3411] | 456 | |
---|
[4663] | 457 | In the future let points be a points file. |
---|
[5454] | 458 | And let the user choose the quantities. |
---|
[3417] | 459 | |
---|
| 460 | This is currently quite specific. |
---|
[5454] | 461 | If it is need to be more general, change things. |
---|
[3445] | 462 | |
---|
[3411] | 463 | """ |
---|
| 464 | quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
| 465 | #print "points",points |
---|
| 466 | points = ensure_absolute(points) |
---|
| 467 | point_count = len(points) |
---|
| 468 | callable_sww = file_function(sww_file, |
---|
| 469 | quantities=quantities, |
---|
| 470 | interpolation_points=points, |
---|
| 471 | verbose=verbose, |
---|
[5447] | 472 | time_thinning=time_thinning, |
---|
[3411] | 473 | use_cache=use_cache) |
---|
| 474 | |
---|
| 475 | depth_writer = writer(file(depth_file, "wb")) |
---|
[3456] | 476 | velocity_x_writer = writer(file(velocity_x_file, "wb")) |
---|
| 477 | velocity_y_writer = writer(file(velocity_y_file, "wb")) |
---|
[5196] | 478 | if stage_file is not None: |
---|
| 479 | stage_writer = writer(file(stage_file, "wb")) |
---|
[5496] | 480 | if froude_file is not None: |
---|
| 481 | froude_writer = writer(file(froude_file, "wb")) |
---|
[3411] | 482 | # Write heading |
---|
[3419] | 483 | heading = [str(x[0])+ ':' + str(x[1]) for x in points] |
---|
| 484 | heading.insert(0, "time") |
---|
| 485 | depth_writer.writerow(heading) |
---|
[3456] | 486 | velocity_x_writer.writerow(heading) |
---|
| 487 | velocity_y_writer.writerow(heading) |
---|
[5196] | 488 | if stage_file is not None: |
---|
[5496] | 489 | stage_writer.writerow(heading) |
---|
| 490 | if froude_file is not None: |
---|
| 491 | froude_writer.writerow(heading) |
---|
[3456] | 492 | |
---|
[3411] | 493 | for time in callable_sww.get_time(): |
---|
[3419] | 494 | depths = [time] |
---|
[3456] | 495 | velocity_xs = [time] |
---|
| 496 | velocity_ys = [time] |
---|
[5196] | 497 | if stage_file is not None: |
---|
| 498 | stages = [time] |
---|
[5496] | 499 | if froude_file is not None: |
---|
| 500 | froudes = [time] |
---|
[3411] | 501 | for point_i, point in enumerate(points): |
---|
| 502 | quantities = callable_sww(time,point_i) |
---|
| 503 | #print "quantities", quantities |
---|
| 504 | |
---|
| 505 | w = quantities[0] |
---|
| 506 | z = quantities[1] |
---|
[3456] | 507 | momentum_x = quantities[2] |
---|
| 508 | momentum_y = quantities[3] |
---|
[3411] | 509 | depth = w - z |
---|
| 510 | |
---|
[3456] | 511 | if w == NAN or z == NAN or momentum_x == NAN: |
---|
| 512 | velocity_x = NAN |
---|
[3411] | 513 | else: |
---|
| 514 | if depth > 1.e-30: # use epsilon |
---|
[3456] | 515 | velocity_x = momentum_x / depth #Absolute velocity |
---|
[3411] | 516 | else: |
---|
[3456] | 517 | velocity_x = 0 |
---|
| 518 | if w == NAN or z == NAN or momentum_y == NAN: |
---|
| 519 | velocity_y = NAN |
---|
| 520 | else: |
---|
| 521 | if depth > 1.e-30: # use epsilon |
---|
| 522 | velocity_y = momentum_y / depth #Absolute velocity |
---|
| 523 | else: |
---|
| 524 | velocity_y = 0 |
---|
[5496] | 525 | if depth < 1.e-30: # use epsilon |
---|
| 526 | froude = NAN |
---|
| 527 | else: |
---|
| 528 | froude = sqrt(velocity_x*velocity_x + velocity_y*velocity_y)/ \ |
---|
| 529 | sqrt(depth * 9.8066) # gravity m/s/s |
---|
[3411] | 530 | depths.append(depth) |
---|
[3456] | 531 | velocity_xs.append(velocity_x) |
---|
| 532 | velocity_ys.append(velocity_y) |
---|
[5196] | 533 | if stage_file is not None: |
---|
| 534 | stages.append(w) |
---|
[5496] | 535 | if froude_file is not None: |
---|
| 536 | froudes.append(froude) |
---|
[3411] | 537 | depth_writer.writerow(depths) |
---|
[3456] | 538 | velocity_x_writer.writerow(velocity_xs) |
---|
| 539 | velocity_y_writer.writerow(velocity_ys) |
---|
[5196] | 540 | if stage_file is not None: |
---|
[5496] | 541 | stage_writer.writerow(stages) |
---|
| 542 | if froude_file is not None: |
---|
| 543 | froude_writer.writerow(froudes) |
---|
[3411] | 544 | |
---|
[3768] | 545 | |
---|
[2879] | 546 | class Interpolation_function: |
---|
| 547 | """Interpolation_interface - creates callable object f(t, id) or f(t,x,y) |
---|
| 548 | which is interpolated from time series defined at vertices of |
---|
| 549 | triangular mesh (such as those stored in sww files) |
---|
| 550 | |
---|
| 551 | Let m be the number of vertices, n the number of triangles |
---|
[5288] | 552 | and p the number of timesteps. |
---|
| 553 | Also, let N be the number of interpolation points. |
---|
[2879] | 554 | |
---|
| 555 | Mandatory input |
---|
| 556 | time: px1 array of monotonously increasing times (Float) |
---|
| 557 | quantities: Dictionary of arrays or 1 array (Float) |
---|
| 558 | The arrays must either have dimensions pxm or mx1. |
---|
| 559 | The resulting function will be time dependent in |
---|
[5288] | 560 | the former case while it will be constant with |
---|
[2879] | 561 | respect to time in the latter case. |
---|
| 562 | |
---|
| 563 | Optional input: |
---|
[5288] | 564 | quantity_names: List of keys into the quantities dictionary for |
---|
| 565 | imposing a particular order on the output vector. |
---|
[2879] | 566 | vertex_coordinates: mx2 array of coordinates (Float) |
---|
| 567 | triangles: nx3 array of indices into vertex_coordinates (Int) |
---|
| 568 | interpolation_points: Nx2 array of coordinates to be interpolated to |
---|
| 569 | verbose: Level of reporting |
---|
| 570 | |
---|
| 571 | |
---|
| 572 | The quantities returned by the callable object are specified by |
---|
| 573 | the list quantities which must contain the names of the |
---|
| 574 | quantities to be returned and also reflect the order, e.g. for |
---|
| 575 | the shallow water wave equation, on would have |
---|
| 576 | quantities = ['stage', 'xmomentum', 'ymomentum'] |
---|
| 577 | |
---|
| 578 | The parameter interpolation_points decides at which points interpolated |
---|
| 579 | quantities are to be computed whenever object is called. |
---|
| 580 | If None, return average value |
---|
[3896] | 581 | |
---|
| 582 | FIXME (Ole): Need to allow vertex coordinates and interpolation points to be |
---|
| 583 | geospatial data objects |
---|
[3900] | 584 | |
---|
[5288] | 585 | Time assumed to be relative to starttime (FIXME (Ole): This comment should be removed) |
---|
| 586 | All coordinates assume origin of (0,0) - e.g. georeferencing must be taken care of |
---|
| 587 | outside this function |
---|
[2879] | 588 | """ |
---|
[3768] | 589 | |
---|
[2879] | 590 | |
---|
| 591 | def __init__(self, |
---|
| 592 | time, |
---|
| 593 | quantities, |
---|
[3850] | 594 | quantity_names=None, |
---|
| 595 | vertex_coordinates=None, |
---|
| 596 | triangles=None, |
---|
| 597 | interpolation_points=None, |
---|
[3900] | 598 | time_thinning=1, |
---|
[5372] | 599 | verbose=False, |
---|
| 600 | gauge_neighbour_id=None): |
---|
[2879] | 601 | """Initialise object and build spatial interpolation if required |
---|
[3900] | 602 | |
---|
| 603 | Time_thinning_number controls how many timesteps to use. Only timesteps with |
---|
| 604 | index%time_thinning_number == 0 will used, or in other words a value of 3, say, |
---|
| 605 | will cause the algorithm to use every third time step. |
---|
[2879] | 606 | """ |
---|
| 607 | |
---|
| 608 | from Numeric import array, zeros, Float, alltrue, concatenate,\ |
---|
| 609 | reshape, ArrayType |
---|
| 610 | |
---|
| 611 | |
---|
[3514] | 612 | from anuga.config import time_format |
---|
[2879] | 613 | import types |
---|
| 614 | |
---|
[4396] | 615 | |
---|
[3850] | 616 | # Check temporal info |
---|
[2879] | 617 | time = ensure_numeric(time) |
---|
| 618 | msg = 'Time must be a monotonuosly ' |
---|
| 619 | msg += 'increasing sequence %s' %time |
---|
| 620 | assert alltrue(time[1:] - time[:-1] >= 0 ), msg |
---|
| 621 | |
---|
| 622 | |
---|
[3850] | 623 | # Check if quantities is a single array only |
---|
[2879] | 624 | if type(quantities) != types.DictType: |
---|
| 625 | quantities = ensure_numeric(quantities) |
---|
| 626 | quantity_names = ['Attribute'] |
---|
| 627 | |
---|
[5566] | 628 | # Make it a dictionary |
---|
[2879] | 629 | quantities = {quantity_names[0]: quantities} |
---|
| 630 | |
---|
| 631 | |
---|
[3850] | 632 | # Use keys if no names are specified |
---|
[2879] | 633 | if quantity_names is None: |
---|
| 634 | quantity_names = quantities.keys() |
---|
| 635 | |
---|
| 636 | |
---|
[3850] | 637 | # Check spatial info |
---|
[2879] | 638 | if vertex_coordinates is None: |
---|
| 639 | self.spatial = False |
---|
[5288] | 640 | else: |
---|
| 641 | # FIXME (Ole): Try ensure_numeric here - |
---|
| 642 | #this function knows nothing about georefering. |
---|
[3941] | 643 | vertex_coordinates = ensure_absolute(vertex_coordinates) |
---|
[2879] | 644 | |
---|
[5361] | 645 | if triangles is not None: |
---|
| 646 | triangles = ensure_numeric(triangles) |
---|
| 647 | self.spatial = True |
---|
[2879] | 648 | |
---|
[3900] | 649 | # Thin timesteps if needed |
---|
| 650 | # Note array() is used to make the thinned arrays contiguous in memory |
---|
| 651 | self.time = array(time[::time_thinning]) |
---|
| 652 | for name in quantity_names: |
---|
| 653 | if len(quantities[name].shape) == 2: |
---|
| 654 | quantities[name] = array(quantities[name][::time_thinning,:]) |
---|
[2879] | 655 | |
---|
[3850] | 656 | # Save for use with statistics |
---|
[3050] | 657 | self.quantities_range = {} |
---|
| 658 | for name in quantity_names: |
---|
| 659 | q = quantities[name][:].flat |
---|
| 660 | self.quantities_range[name] = [min(q), max(q)] |
---|
| 661 | |
---|
[2879] | 662 | self.quantity_names = quantity_names |
---|
| 663 | self.vertex_coordinates = vertex_coordinates |
---|
| 664 | self.interpolation_points = interpolation_points |
---|
[3900] | 665 | |
---|
| 666 | |
---|
[2879] | 667 | self.index = 0 # Initial time index |
---|
| 668 | self.precomputed_values = {} |
---|
[3050] | 669 | |
---|
[4396] | 670 | |
---|
[3850] | 671 | # Precomputed spatial interpolation if requested |
---|
[2879] | 672 | if interpolation_points is not None: |
---|
[5361] | 673 | #no longer true. sts files have spatial = True but |
---|
| 674 | #if self.spatial is False: |
---|
| 675 | # raise 'Triangles and vertex_coordinates must be specified' |
---|
| 676 | # |
---|
[2879] | 677 | try: |
---|
[4833] | 678 | self.interpolation_points = interpolation_points = ensure_numeric(interpolation_points) |
---|
[2879] | 679 | except: |
---|
| 680 | msg = 'Interpolation points must be an N x 2 Numeric array '+\ |
---|
| 681 | 'or a list of points\n' |
---|
| 682 | msg += 'I got: %s.' %(str(self.interpolation_points)[:60] +\ |
---|
| 683 | '...') |
---|
| 684 | raise msg |
---|
| 685 | |
---|
[5361] | 686 | if triangles is not None and vertex_coordinates is not None: |
---|
| 687 | # Check that all interpolation points fall within |
---|
| 688 | # mesh boundary as defined by triangles and vertex_coordinates. |
---|
| 689 | from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh |
---|
| 690 | from anuga.utilities.polygon import outside_polygon |
---|
[2879] | 691 | |
---|
[5361] | 692 | # Create temporary mesh object from mesh info passed |
---|
| 693 | # into this function. |
---|
| 694 | mesh = Mesh(vertex_coordinates, triangles) |
---|
| 695 | mesh_boundary_polygon = mesh.get_boundary_polygon() |
---|
[4833] | 696 | |
---|
[4820] | 697 | |
---|
[5361] | 698 | indices = outside_polygon(interpolation_points, |
---|
| 699 | mesh_boundary_polygon) |
---|
[4820] | 700 | |
---|
[5361] | 701 | # Record result |
---|
| 702 | #self.mesh_boundary_polygon = mesh_boundary_polygon |
---|
| 703 | self.indices_outside_mesh = indices |
---|
[4820] | 704 | |
---|
[5361] | 705 | # Report |
---|
| 706 | if len(indices) > 0: |
---|
| 707 | msg = 'Interpolation points in Interpolation function fall ' |
---|
| 708 | msg += 'outside specified mesh. ' |
---|
| 709 | msg += 'Offending points:\n' |
---|
| 710 | out_interp_pts = [] |
---|
| 711 | for i in indices: |
---|
| 712 | msg += '%d: %s\n' %(i, interpolation_points[i]) |
---|
| 713 | out_interp_pts.append(ensure_numeric(interpolation_points[i])) |
---|
[4820] | 714 | |
---|
[5361] | 715 | if verbose is True: |
---|
| 716 | import sys |
---|
| 717 | if sys.platform == 'win32': |
---|
| 718 | from anuga.utilities.polygon import plot_polygons |
---|
| 719 | #out_interp_pts = take(interpolation_points,[indices]) |
---|
| 720 | title = 'Interpolation points fall outside specified mesh' |
---|
| 721 | plot_polygons([mesh_boundary_polygon,interpolation_points,out_interp_pts], |
---|
| 722 | ['line','point','outside'],figname='points_boundary_out',label=title,verbose=verbose) |
---|
[4851] | 723 | |
---|
[5361] | 724 | # Joaquim Luis suggested this as an Exception, so |
---|
| 725 | # that the user can now what the problem is rather than |
---|
| 726 | # looking for NaN's. However, NANs are handy as they can |
---|
| 727 | # be ignored leaving good points for continued processing. |
---|
| 728 | if verbose: |
---|
| 729 | print msg |
---|
| 730 | #raise Exception(msg) |
---|
| 731 | elif triangles is None and vertex_coordinates is not None:#jj |
---|
| 732 | #Dealing with sts file |
---|
| 733 | pass |
---|
| 734 | else: |
---|
| 735 | msg = 'Sww file function requires both triangles and vertex_coordinates. sts file file function requires the later.' |
---|
| 736 | raise Exception(msg) |
---|
[4820] | 737 | |
---|
[4833] | 738 | # Plot boundary and interpolation points |
---|
[4845] | 739 | if verbose is True: |
---|
[4869] | 740 | import sys |
---|
| 741 | if sys.platform == 'win32': |
---|
| 742 | from anuga.utilities.polygon import plot_polygons |
---|
| 743 | title = 'Interpolation function: Polygon and interpolation points' |
---|
| 744 | plot_polygons([mesh_boundary_polygon,interpolation_points], |
---|
| 745 | ['line','point'],figname='points_boundary',label=title,verbose=verbose) |
---|
[4820] | 746 | |
---|
[2879] | 747 | m = len(self.interpolation_points) |
---|
[2884] | 748 | p = len(self.time) |
---|
[4396] | 749 | |
---|
[2879] | 750 | for name in quantity_names: |
---|
| 751 | self.precomputed_values[name] = zeros((p, m), Float) |
---|
| 752 | |
---|
[3850] | 753 | # Build interpolator |
---|
[3900] | 754 | if verbose: |
---|
[5361] | 755 | if triangles is not None and vertex_coordinates is not None: |
---|
| 756 | msg = 'Building interpolation matrix from source mesh ' |
---|
| 757 | msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0], |
---|
| 758 | triangles.shape[0]) |
---|
| 759 | elif triangles is None and vertex_coordinates is not None: |
---|
| 760 | msg = 'Building interpolation matrix from source points' |
---|
| 761 | |
---|
[3900] | 762 | print msg |
---|
| 763 | |
---|
| 764 | |
---|
[2879] | 765 | interpol = Interpolate(vertex_coordinates, |
---|
[3689] | 766 | triangles, |
---|
[3850] | 767 | verbose=verbose) |
---|
[2879] | 768 | |
---|
[3900] | 769 | if verbose: |
---|
| 770 | print 'Interpolating (%d interpolation points, %d timesteps).'\ |
---|
| 771 | %(self.interpolation_points.shape[0], self.time.shape[0]), |
---|
| 772 | |
---|
| 773 | if time_thinning > 1: |
---|
| 774 | print 'Timesteps were thinned by a factor of %d' %time_thinning |
---|
| 775 | else: |
---|
| 776 | print |
---|
| 777 | |
---|
[2884] | 778 | for i, t in enumerate(self.time): |
---|
[3768] | 779 | # Interpolate quantities at this timestep |
---|
[2879] | 780 | if verbose and i%((p+10)/10)==0: |
---|
[5222] | 781 | print ' time step %d of %d' %(i, p) |
---|
[2879] | 782 | |
---|
| 783 | for name in quantity_names: |
---|
| 784 | if len(quantities[name].shape) == 2: |
---|
[3768] | 785 | Q = quantities[name][i,:] # Quantities at timestep i |
---|
[2879] | 786 | else: |
---|
[3768] | 787 | Q = quantities[name][:] # No time dependency |
---|
[5222] | 788 | |
---|
| 789 | if verbose and i%((p+10)/10)==0: |
---|
| 790 | print ' quantity %s, size=%d' %(name, len(Q)) |
---|
| 791 | |
---|
[5361] | 792 | # Interpolate |
---|
| 793 | if triangles is not None and vertex_coordinates is not None: |
---|
| 794 | result = interpol.interpolate(Q, |
---|
| 795 | point_coordinates=\ |
---|
| 796 | self.interpolation_points, |
---|
| 797 | verbose=False) # Don't clutter |
---|
| 798 | elif triangles is None and vertex_coordinates is not None: |
---|
[5372] | 799 | result=interpol.interpolate_polyline(Q,vertex_coordinates,gauge_neighbour_id,point_coordinates=self.interpolation_points) |
---|
[4588] | 800 | |
---|
| 801 | #assert len(result), len(interpolation_points) |
---|
[2879] | 802 | self.precomputed_values[name][i, :] = result |
---|
[3768] | 803 | |
---|
[2879] | 804 | |
---|
[3768] | 805 | # Report |
---|
[2879] | 806 | if verbose: |
---|
| 807 | print self.statistics() |
---|
| 808 | #self.print_statistics() |
---|
| 809 | |
---|
| 810 | else: |
---|
[3768] | 811 | # Store quantitites as is |
---|
| 812 | for name in quantity_names: |
---|
[2879] | 813 | self.precomputed_values[name] = quantities[name] |
---|
| 814 | |
---|
| 815 | def __repr__(self): |
---|
[3768] | 816 | # return 'Interpolation function (spatio-temporal)' |
---|
[2879] | 817 | return self.statistics() |
---|
[3768] | 818 | |
---|
| 819 | def __call__(self, t, point_id=None, x=None, y=None): |
---|
[4370] | 820 | """Evaluate f(t) or f(t, point_id) |
---|
| 821 | |
---|
[2879] | 822 | Inputs: |
---|
| 823 | t: time - Model time. Must lie within existing timesteps |
---|
| 824 | point_id: index of one of the preprocessed points. |
---|
[4370] | 825 | |
---|
| 826 | |
---|
| 827 | If spatial info is present and all of point_id |
---|
[2879] | 828 | are None an exception is raised |
---|
| 829 | |
---|
[4370] | 830 | If no spatial info is present, point_id arguments are ignored |
---|
[2879] | 831 | making f a function of time only. |
---|
| 832 | |
---|
[4370] | 833 | |
---|
[4548] | 834 | FIXME: f(t, x, y) x, y could overrided location, point_id ignored |
---|
[2879] | 835 | FIXME: point_id could also be a slice |
---|
| 836 | FIXME: What if x and y are vectors? |
---|
| 837 | FIXME: What about f(x,y) without t? |
---|
| 838 | """ |
---|
| 839 | |
---|
| 840 | from math import pi, cos, sin, sqrt |
---|
| 841 | from Numeric import zeros, Float |
---|
[3560] | 842 | from anuga.abstract_2d_finite_volumes.util import mean |
---|
[2879] | 843 | |
---|
| 844 | if self.spatial is True: |
---|
| 845 | if point_id is None: |
---|
| 846 | if x is None or y is None: |
---|
| 847 | msg = 'Either point_id or x and y must be specified' |
---|
| 848 | raise Exception(msg) |
---|
| 849 | else: |
---|
| 850 | if self.interpolation_points is None: |
---|
| 851 | msg = 'Interpolation_function must be instantiated ' +\ |
---|
| 852 | 'with a list of interpolation points before parameter ' +\ |
---|
| 853 | 'point_id can be used' |
---|
| 854 | raise Exception(msg) |
---|
| 855 | |
---|
[2884] | 856 | msg = 'Time interval [%.16f:%.16f]' %(self.time[0], self.time[-1]) |
---|
[2879] | 857 | msg += ' does not match model time: %.16f\n' %t |
---|
[2884] | 858 | if t < self.time[0]: raise Exception(msg) |
---|
| 859 | if t > self.time[-1]: raise Exception(msg) |
---|
[2879] | 860 | |
---|
| 861 | oldindex = self.index #Time index |
---|
| 862 | |
---|
[4820] | 863 | # Find current time slot |
---|
[2884] | 864 | while t > self.time[self.index]: self.index += 1 |
---|
| 865 | while t < self.time[self.index]: self.index -= 1 |
---|
[2879] | 866 | |
---|
[2884] | 867 | if t == self.time[self.index]: |
---|
[4820] | 868 | # Protect against case where t == T[-1] (last time) |
---|
| 869 | # - also works in general when t == T[i] |
---|
[2879] | 870 | ratio = 0 |
---|
| 871 | else: |
---|
[4820] | 872 | # t is now between index and index+1 |
---|
[2884] | 873 | ratio = (t - self.time[self.index])/\ |
---|
| 874 | (self.time[self.index+1] - self.time[self.index]) |
---|
[2879] | 875 | |
---|
[4820] | 876 | # Compute interpolated values |
---|
[2879] | 877 | q = zeros(len(self.quantity_names), Float) |
---|
[4820] | 878 | # print "self.precomputed_values", self.precomputed_values |
---|
[2879] | 879 | for i, name in enumerate(self.quantity_names): |
---|
| 880 | Q = self.precomputed_values[name] |
---|
| 881 | |
---|
| 882 | if self.spatial is False: |
---|
[4820] | 883 | # If there is no spatial info |
---|
[2879] | 884 | assert len(Q.shape) == 1 |
---|
| 885 | |
---|
| 886 | Q0 = Q[self.index] |
---|
| 887 | if ratio > 0: Q1 = Q[self.index+1] |
---|
| 888 | |
---|
| 889 | else: |
---|
| 890 | if x is not None and y is not None: |
---|
[4820] | 891 | # Interpolate to x, y |
---|
[2879] | 892 | |
---|
| 893 | raise 'x,y interpolation not yet implemented' |
---|
| 894 | else: |
---|
[4820] | 895 | # Use precomputed point |
---|
[2879] | 896 | Q0 = Q[self.index, point_id] |
---|
| 897 | if ratio > 0: |
---|
| 898 | Q1 = Q[self.index+1, point_id] |
---|
[4588] | 899 | |
---|
[4820] | 900 | # Linear temporal interpolation |
---|
[2879] | 901 | if ratio > 0: |
---|
[3452] | 902 | if Q0 == NAN and Q1 == NAN: |
---|
[5482] | 903 | q[i] = Q0 |
---|
[2879] | 904 | else: |
---|
| 905 | q[i] = Q0 + ratio*(Q1 - Q0) |
---|
| 906 | else: |
---|
| 907 | q[i] = Q0 |
---|
| 908 | |
---|
| 909 | |
---|
[4820] | 910 | # Return vector of interpolated values |
---|
| 911 | # if len(q) == 1: |
---|
| 912 | # return q[0] |
---|
| 913 | # else: |
---|
| 914 | # return q |
---|
[2879] | 915 | |
---|
| 916 | |
---|
[4820] | 917 | # Return vector of interpolated values |
---|
| 918 | # FIXME: |
---|
[2879] | 919 | if self.spatial is True: |
---|
| 920 | return q |
---|
| 921 | else: |
---|
[4820] | 922 | # Replicate q according to x and y |
---|
| 923 | # This is e.g used for Wind_stress |
---|
[2879] | 924 | if x is None or y is None: |
---|
| 925 | return q |
---|
| 926 | else: |
---|
| 927 | try: |
---|
| 928 | N = len(x) |
---|
| 929 | except: |
---|
| 930 | return q |
---|
| 931 | else: |
---|
| 932 | from Numeric import ones, Float |
---|
[4820] | 933 | # x is a vector - Create one constant column for each value |
---|
[2879] | 934 | N = len(x) |
---|
| 935 | assert len(y) == N, 'x and y must have same length' |
---|
| 936 | res = [] |
---|
| 937 | for col in q: |
---|
| 938 | res.append(col*ones(N, Float)) |
---|
| 939 | |
---|
| 940 | return res |
---|
| 941 | |
---|
| 942 | |
---|
[2884] | 943 | def get_time(self): |
---|
| 944 | """Return model time as a vector of timesteps |
---|
| 945 | """ |
---|
| 946 | return self.time |
---|
| 947 | |
---|
[3768] | 948 | |
---|
[2879] | 949 | def statistics(self): |
---|
| 950 | """Output statistics about interpolation_function |
---|
| 951 | """ |
---|
| 952 | |
---|
| 953 | vertex_coordinates = self.vertex_coordinates |
---|
| 954 | interpolation_points = self.interpolation_points |
---|
| 955 | quantity_names = self.quantity_names |
---|
[3055] | 956 | #quantities = self.quantities |
---|
[2879] | 957 | precomputed_values = self.precomputed_values |
---|
| 958 | |
---|
| 959 | x = vertex_coordinates[:,0] |
---|
| 960 | y = vertex_coordinates[:,1] |
---|
| 961 | |
---|
| 962 | str = '------------------------------------------------\n' |
---|
| 963 | str += 'Interpolation_function (spatio-temporal) statistics:\n' |
---|
| 964 | str += ' Extent:\n' |
---|
| 965 | str += ' x in [%f, %f], len(x) == %d\n'\ |
---|
| 966 | %(min(x), max(x), len(x)) |
---|
| 967 | str += ' y in [%f, %f], len(y) == %d\n'\ |
---|
| 968 | %(min(y), max(y), len(y)) |
---|
| 969 | str += ' t in [%f, %f], len(t) == %d\n'\ |
---|
[2884] | 970 | %(min(self.time), max(self.time), len(self.time)) |
---|
[2879] | 971 | str += ' Quantities:\n' |
---|
| 972 | for name in quantity_names: |
---|
[3050] | 973 | minq, maxq = self.quantities_range[name] |
---|
| 974 | str += ' %s in [%f, %f]\n' %(name, minq, maxq) |
---|
| 975 | #q = quantities[name][:].flat |
---|
| 976 | #str += ' %s in [%f, %f]\n' %(name, min(q), max(q)) |
---|
[2879] | 977 | |
---|
| 978 | if interpolation_points is not None: |
---|
| 979 | str += ' Interpolation points (xi, eta):'\ |
---|
| 980 | ' number of points == %d\n' %interpolation_points.shape[0] |
---|
| 981 | str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]), |
---|
| 982 | max(interpolation_points[:,0])) |
---|
| 983 | str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]), |
---|
| 984 | max(interpolation_points[:,1])) |
---|
| 985 | str += ' Interpolated quantities (over all timesteps):\n' |
---|
| 986 | |
---|
| 987 | for name in quantity_names: |
---|
| 988 | q = precomputed_values[name][:].flat |
---|
| 989 | str += ' %s at interpolation points in [%f, %f]\n'\ |
---|
| 990 | %(name, min(q), max(q)) |
---|
| 991 | str += '------------------------------------------------\n' |
---|
| 992 | |
---|
| 993 | return str |
---|
| 994 | |
---|
[3768] | 995 | |
---|
[2879] | 996 | def interpolate_sww(sww_file, time, interpolation_points, |
---|
[3896] | 997 | quantity_names=None, verbose=False): |
---|
[2879] | 998 | """ |
---|
| 999 | obsolete. |
---|
| 1000 | use file_function in utils |
---|
| 1001 | """ |
---|
| 1002 | #open sww file |
---|
| 1003 | x, y, volumes, time, quantities = read_sww(sww_file) |
---|
| 1004 | print "x",x |
---|
| 1005 | print "y",y |
---|
| 1006 | |
---|
| 1007 | print "time", time |
---|
| 1008 | print "quantities", quantities |
---|
| 1009 | |
---|
| 1010 | #Add the x and y together |
---|
| 1011 | vertex_coordinates = concatenate((x[:,NewAxis], y[:,NewAxis]),axis=1) |
---|
| 1012 | |
---|
| 1013 | #Will return the quantity values at the specified times and locations |
---|
[3768] | 1014 | interp = Interpolation_interface(time, |
---|
| 1015 | quantities, |
---|
| 1016 | quantity_names=quantity_names, |
---|
| 1017 | vertex_coordinates=vertex_coordinates, |
---|
| 1018 | triangles=volumes, |
---|
| 1019 | interpolation_points=interpolation_points, |
---|
| 1020 | verbose=verbose) |
---|
[2879] | 1021 | |
---|
| 1022 | |
---|
| 1023 | def read_sww(file_name): |
---|
| 1024 | """ |
---|
| 1025 | obsolete - Nothing should be calling this |
---|
| 1026 | |
---|
| 1027 | Read in an sww file. |
---|
| 1028 | |
---|
| 1029 | Input; |
---|
| 1030 | file_name - the sww file |
---|
| 1031 | |
---|
| 1032 | Output; |
---|
| 1033 | x - Vector of x values |
---|
| 1034 | y - Vector of y values |
---|
| 1035 | z - Vector of bed elevation |
---|
| 1036 | volumes - Array. Each row has 3 values, representing |
---|
| 1037 | the vertices that define the volume |
---|
| 1038 | time - Vector of the times where there is stage information |
---|
| 1039 | stage - array with respect to time and vertices (x,y) |
---|
| 1040 | """ |
---|
[4548] | 1041 | |
---|
| 1042 | msg = 'Function read_sww in interpolat.py is obsolete' |
---|
| 1043 | raise Exception, msg |
---|
[2879] | 1044 | |
---|
| 1045 | #FIXME Have this reader as part of data_manager? |
---|
| 1046 | |
---|
| 1047 | from Scientific.IO.NetCDF import NetCDFFile |
---|
| 1048 | import tempfile |
---|
| 1049 | import sys |
---|
| 1050 | import os |
---|
| 1051 | |
---|
| 1052 | #Check contents |
---|
| 1053 | #Get NetCDF |
---|
| 1054 | |
---|
| 1055 | # see if the file is there. Throw a QUIET IO error if it isn't |
---|
| 1056 | # I don't think I could implement the above |
---|
| 1057 | |
---|
| 1058 | #throws prints to screen if file not present |
---|
| 1059 | junk = tempfile.mktemp(".txt") |
---|
| 1060 | fd = open(junk,'w') |
---|
| 1061 | stdout = sys.stdout |
---|
| 1062 | sys.stdout = fd |
---|
| 1063 | fid = NetCDFFile(file_name, 'r') |
---|
| 1064 | sys.stdout = stdout |
---|
| 1065 | fd.close() |
---|
| 1066 | #clean up |
---|
| 1067 | os.remove(junk) |
---|
| 1068 | |
---|
| 1069 | # Get the variables |
---|
| 1070 | x = fid.variables['x'][:] |
---|
| 1071 | y = fid.variables['y'][:] |
---|
| 1072 | volumes = fid.variables['volumes'][:] |
---|
| 1073 | time = fid.variables['time'][:] |
---|
| 1074 | |
---|
| 1075 | keys = fid.variables.keys() |
---|
| 1076 | keys.remove("x") |
---|
| 1077 | keys.remove("y") |
---|
| 1078 | keys.remove("volumes") |
---|
| 1079 | keys.remove("time") |
---|
| 1080 | #Turn NetCDF objects into Numeric arrays |
---|
| 1081 | quantities = {} |
---|
| 1082 | for name in keys: |
---|
| 1083 | quantities[name] = fid.variables[name][:] |
---|
| 1084 | |
---|
| 1085 | fid.close() |
---|
| 1086 | return x, y, volumes, time, quantities |
---|
| 1087 | |
---|
[3768] | 1088 | |
---|
[2879] | 1089 | #------------------------------------------------------------- |
---|
| 1090 | if __name__ == "__main__": |
---|
[3411] | 1091 | names = ["x","y"] |
---|
| 1092 | someiterable = [[1,2],[3,4]] |
---|
| 1093 | csvwriter = writer(file("some.csv", "wb")) |
---|
| 1094 | csvwriter.writerow(names) |
---|
| 1095 | for row in someiterable: |
---|
| 1096 | csvwriter.writerow(row) |
---|