[6073] | 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 | import sys |
---|
| 23 | from warnings import warn |
---|
| 24 | from math import sqrt |
---|
| 25 | from csv import writer, DictWriter |
---|
| 26 | |
---|
| 27 | from anuga.caching.caching import cache |
---|
| 28 | from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh |
---|
| 29 | from anuga.utilities.sparse import Sparse, Sparse_CSR |
---|
| 30 | from anuga.utilities.cg_solve import conjugate_gradient, VectorShapeError |
---|
| 31 | from anuga.coordinate_transforms.geo_reference import Geo_reference |
---|
| 32 | from anuga.utilities.numerical_tools import ensure_numeric, mean, NAN |
---|
| 33 | from anuga.utilities.polygon import in_and_outside_polygon |
---|
| 34 | from anuga.geospatial_data.geospatial_data import Geospatial_data |
---|
| 35 | from anuga.geospatial_data.geospatial_data import ensure_absolute |
---|
| 36 | from anuga.fit_interpolate.search_functions import search_tree_of_vertices |
---|
| 37 | from anuga.fit_interpolate.general_fit_interpolate import FitInterpolate |
---|
| 38 | from anuga.abstract_2d_finite_volumes.util import file_function |
---|
[6086] | 39 | from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a |
---|
[6189] | 40 | from utilities.polygon import interpolate_polyline |
---|
[7317] | 41 | import anuga.utilities.log as log |
---|
[6073] | 42 | |
---|
[7276] | 43 | import numpy as num |
---|
[6152] | 44 | |
---|
| 45 | |
---|
[6073] | 46 | # Interpolation specific exceptions |
---|
| 47 | |
---|
| 48 | class Modeltime_too_late(Exception): pass |
---|
| 49 | class Modeltime_too_early(Exception): pass |
---|
| 50 | |
---|
| 51 | |
---|
| 52 | ## |
---|
| 53 | # @brief Interpolate vertex_values to interpolation points. |
---|
| 54 | # @param vertex_coordinates List of coordinate pairs making a mesh. |
---|
| 55 | # @param triangles Iterable of 3-tuples representing indices of mesh vertices. |
---|
| 56 | # @param vertex_values Array of data at mesh vertices. |
---|
| 57 | # @param interpolation_points Points to interpolate to. |
---|
| 58 | # @param mesh_origin A geo_ref object or 3-tuples of UTMzone, easting, northing. |
---|
| 59 | # @param max_vertices_per_cell Max number of vertices before splitting cell. |
---|
| 60 | # @param start_blocking_len Block if # of points greater than this. |
---|
| 61 | # @param use_cache If True, cache. |
---|
| 62 | # @param verbose True if this function is to be verbose. |
---|
| 63 | def interpolate(vertex_coordinates, |
---|
| 64 | triangles, |
---|
| 65 | vertex_values, |
---|
| 66 | interpolation_points, |
---|
| 67 | mesh_origin=None, |
---|
| 68 | max_vertices_per_cell=None, |
---|
| 69 | start_blocking_len=500000, |
---|
[7276] | 70 | use_cache=False, |
---|
[7673] | 71 | verbose=False, |
---|
| 72 | output_centroids=False): |
---|
[6073] | 73 | """Interpolate vertex_values to interpolation points. |
---|
[7276] | 74 | |
---|
[6073] | 75 | Inputs (mandatory): |
---|
| 76 | |
---|
[7276] | 77 | |
---|
[6073] | 78 | vertex_coordinates: List of coordinate pairs [xi, eta] of |
---|
[7276] | 79 | points constituting a mesh |
---|
| 80 | (or an m x 2 numeric array or |
---|
[6073] | 81 | a geospatial object) |
---|
| 82 | Points may appear multiple times |
---|
| 83 | (e.g. if vertices have discontinuities) |
---|
| 84 | |
---|
[7276] | 85 | triangles: List of 3-tuples (or a numeric array) of |
---|
| 86 | integers representing indices of all vertices |
---|
[6073] | 87 | in the mesh. |
---|
[7276] | 88 | |
---|
[6073] | 89 | vertex_values: Vector or array of data at the mesh vertices. |
---|
| 90 | If array, interpolation will be done for each column as |
---|
| 91 | per underlying matrix-matrix multiplication |
---|
[7276] | 92 | |
---|
[6073] | 93 | interpolation_points: Interpolate mesh data to these positions. |
---|
| 94 | List of coordinate pairs [x, y] of |
---|
[7276] | 95 | data points or an nx2 numeric array or a |
---|
[6073] | 96 | Geospatial_data object |
---|
[7276] | 97 | |
---|
[6073] | 98 | Inputs (optional) |
---|
[7276] | 99 | |
---|
[6073] | 100 | mesh_origin: A geo_reference object or 3-tuples consisting of |
---|
| 101 | UTM zone, easting and northing. |
---|
| 102 | If specified vertex coordinates are assumed to be |
---|
| 103 | relative to their respective origins. |
---|
| 104 | |
---|
| 105 | max_vertices_per_cell: Number of vertices in a quad tree cell |
---|
| 106 | at which the cell is split into 4. |
---|
| 107 | |
---|
| 108 | Note: Don't supply a vertex coords as a geospatial |
---|
| 109 | object and a mesh origin, since geospatial has its |
---|
| 110 | own mesh origin. |
---|
[7276] | 111 | |
---|
[6073] | 112 | start_blocking_len: If the # of points is more or greater than this, |
---|
[7276] | 113 | start blocking |
---|
| 114 | |
---|
[6073] | 115 | use_cache: True or False |
---|
| 116 | |
---|
[7276] | 117 | |
---|
[6073] | 118 | Output: |
---|
[7276] | 119 | |
---|
[6073] | 120 | Interpolated values at specified point_coordinates |
---|
| 121 | |
---|
[7276] | 122 | Note: This function is a simple shortcut for case where |
---|
[6073] | 123 | interpolation matrix is unnecessary |
---|
[7276] | 124 | Note: This function does not take blocking into account, |
---|
[6073] | 125 | but allows caching. |
---|
[7276] | 126 | |
---|
[6073] | 127 | """ |
---|
| 128 | |
---|
| 129 | # FIXME(Ole): Probably obsolete since I is precomputed and |
---|
| 130 | # interpolate_block caches |
---|
[7276] | 131 | |
---|
[6073] | 132 | from anuga.caching import cache |
---|
| 133 | |
---|
| 134 | # Create interpolation object with matrix |
---|
[7276] | 135 | args = (ensure_numeric(vertex_coordinates, num.float), |
---|
[6073] | 136 | ensure_numeric(triangles)) |
---|
| 137 | kwargs = {'mesh_origin': mesh_origin, |
---|
| 138 | 'max_vertices_per_cell': max_vertices_per_cell, |
---|
| 139 | 'verbose': verbose} |
---|
[7276] | 140 | |
---|
[6073] | 141 | if use_cache is True: |
---|
| 142 | if sys.platform != 'win32': |
---|
| 143 | I = cache(Interpolate, args, kwargs, verbose=verbose) |
---|
| 144 | else: |
---|
| 145 | # Messy wrapping of Interpolate to deal with win32 error |
---|
| 146 | def wrap_Interpolate(args,kwargs): |
---|
| 147 | I = apply(Interpolate, args, kwargs) |
---|
| 148 | return I |
---|
| 149 | I = cache(wrap_Interpolate, (args, kwargs), {}, verbose=verbose) |
---|
| 150 | else: |
---|
| 151 | I = apply(Interpolate, args, kwargs) |
---|
[7276] | 152 | |
---|
[6073] | 153 | # Call interpolate method with interpolation points |
---|
| 154 | result = I.interpolate_block(vertex_values, interpolation_points, |
---|
| 155 | use_cache=use_cache, |
---|
[7673] | 156 | verbose=verbose, |
---|
| 157 | output_centroids=output_centroids) |
---|
[7276] | 158 | |
---|
[6073] | 159 | return result |
---|
| 160 | |
---|
| 161 | |
---|
| 162 | ## |
---|
[7276] | 163 | # @brief |
---|
[6073] | 164 | class Interpolate (FitInterpolate): |
---|
| 165 | |
---|
| 166 | ## |
---|
| 167 | # @brief Build interpolation matrix. |
---|
| 168 | # @param vertex_coordinates List of pairs [xi, eta] of points making a mesh. |
---|
| 169 | # @param triangles List of 3-tuples of indices of all vertices in the mesh. |
---|
| 170 | # @param mesh_origin A geo_ref object (UTM zone, easting and northing). |
---|
| 171 | # @param verbose If True, this function is to be verbose. |
---|
| 172 | # @param max_vertices_per_cell Split quadtree cell if vertices >= this. |
---|
| 173 | def __init__(self, |
---|
| 174 | vertex_coordinates, |
---|
| 175 | triangles, |
---|
| 176 | mesh_origin=None, |
---|
| 177 | verbose=False, |
---|
| 178 | max_vertices_per_cell=None): |
---|
| 179 | |
---|
| 180 | """ Build interpolation matrix mapping from |
---|
| 181 | function values at vertices to function values at data points |
---|
| 182 | |
---|
| 183 | Inputs: |
---|
| 184 | vertex_coordinates: List of coordinate pairs [xi, eta] of |
---|
[7276] | 185 | points constituting a mesh (or an m x 2 numeric array or |
---|
[6073] | 186 | a geospatial object) |
---|
| 187 | Points may appear multiple times |
---|
| 188 | (e.g. if vertices have discontinuities) |
---|
| 189 | |
---|
[7276] | 190 | triangles: List of 3-tuples (or a numeric array) of |
---|
[6073] | 191 | integers representing indices of all vertices in the mesh. |
---|
| 192 | |
---|
| 193 | mesh_origin: A geo_reference object or 3-tuples consisting of |
---|
| 194 | UTM zone, easting and northing. |
---|
| 195 | If specified vertex coordinates are assumed to be |
---|
| 196 | relative to their respective origins. |
---|
| 197 | |
---|
| 198 | max_vertices_per_cell: Number of vertices in a quad tree cell |
---|
| 199 | at which the cell is split into 4. |
---|
| 200 | |
---|
| 201 | Note: Don't supply a vertex coords as a geospatial object and |
---|
| 202 | a mesh origin, since geospatial has its own mesh origin. |
---|
| 203 | """ |
---|
| 204 | |
---|
| 205 | # FIXME (Ole): Need an input check |
---|
[7276] | 206 | |
---|
[6073] | 207 | FitInterpolate.__init__(self, |
---|
| 208 | vertex_coordinates=vertex_coordinates, |
---|
| 209 | triangles=triangles, |
---|
| 210 | mesh_origin=mesh_origin, |
---|
| 211 | verbose=verbose, |
---|
| 212 | max_vertices_per_cell=max_vertices_per_cell) |
---|
[7276] | 213 | |
---|
[6073] | 214 | # Initialise variables |
---|
| 215 | self._A_can_be_reused = False # FIXME (Ole): Probably obsolete |
---|
| 216 | self._point_coordinates = None # FIXME (Ole): Probably obsolete |
---|
| 217 | self.interpolation_matrices = {} # Store precomputed matrices |
---|
| 218 | |
---|
| 219 | |
---|
| 220 | ## |
---|
| 221 | # @brief Interpolate mesh data f to determine values, z, at points. |
---|
| 222 | # @param f Data on the mesh vertices. |
---|
| 223 | # @param point_coordinates Interpolate mesh data to these positions. |
---|
| 224 | # @param start_blocking_len Block if # points >= this. |
---|
| 225 | # @param verbose True if this function is to be verbose. |
---|
| 226 | # FIXME: What is a good start_blocking_len value? |
---|
| 227 | def interpolate(self, |
---|
| 228 | f, |
---|
| 229 | point_coordinates=None, |
---|
| 230 | start_blocking_len=500000, |
---|
[7673] | 231 | verbose=False, |
---|
| 232 | output_centroids=False): |
---|
[6073] | 233 | """Interpolate mesh data f to determine values, z, at points. |
---|
| 234 | |
---|
| 235 | f is the data on the mesh vertices. |
---|
| 236 | |
---|
| 237 | The mesh values representing a smooth surface are |
---|
| 238 | assumed to be specified in f. |
---|
| 239 | |
---|
| 240 | Inputs: |
---|
| 241 | f: Vector or array of data at the mesh vertices. |
---|
| 242 | If f is an array, interpolation will be done for each column as |
---|
| 243 | per underlying matrix-matrix multiplication |
---|
| 244 | |
---|
| 245 | point_coordinates: Interpolate mesh data to these positions. |
---|
| 246 | List of coordinate pairs [x, y] of |
---|
[7276] | 247 | data points or an nx2 numeric array or a Geospatial_data object |
---|
| 248 | |
---|
| 249 | If point_coordinates is absent, the points inputted last time |
---|
[6073] | 250 | this method was called are used, if possible. |
---|
| 251 | |
---|
| 252 | start_blocking_len: If the # of points is more or greater than this, |
---|
[7276] | 253 | start blocking |
---|
[6073] | 254 | |
---|
[7276] | 255 | Output: |
---|
| 256 | Interpolated values at inputted points (z). |
---|
[6073] | 257 | """ |
---|
| 258 | |
---|
| 259 | # FIXME (Ole): Why is the interpolation matrix rebuilt everytime the |
---|
| 260 | # method is called even if interpolation points are unchanged. |
---|
| 261 | # This really should use some kind of caching in cases where |
---|
| 262 | # interpolation points are reused. |
---|
| 263 | # |
---|
| 264 | # This has now been addressed through an attempt in interpolate_block |
---|
| 265 | |
---|
[7317] | 266 | if verbose: log.critical('Build intepolation object') |
---|
[6073] | 267 | if isinstance(point_coordinates, Geospatial_data): |
---|
| 268 | point_coordinates = point_coordinates.get_data_points(absolute=True) |
---|
| 269 | |
---|
| 270 | # Can I interpolate, based on previous point_coordinates? |
---|
| 271 | if point_coordinates is None: |
---|
| 272 | if self._A_can_be_reused is True \ |
---|
| 273 | and len(self._point_coordinates) < start_blocking_len: |
---|
| 274 | z = self._get_point_data_z(f, verbose=verbose) |
---|
| 275 | elif self._point_coordinates is not None: |
---|
| 276 | # if verbose, give warning |
---|
| 277 | if verbose: |
---|
[7317] | 278 | log.critical('WARNING: Recalculating A matrix, ' |
---|
| 279 | 'due to blocking.') |
---|
[6073] | 280 | point_coordinates = self._point_coordinates |
---|
| 281 | else: |
---|
[6184] | 282 | # There are no good point_coordinates. import sys; sys.exit() |
---|
[6073] | 283 | msg = 'ERROR (interpolate.py): No point_coordinates inputted' |
---|
| 284 | raise Exception(msg) |
---|
[7276] | 285 | |
---|
[6073] | 286 | if point_coordinates is not None: |
---|
| 287 | self._point_coordinates = point_coordinates |
---|
| 288 | if len(point_coordinates) < start_blocking_len \ |
---|
| 289 | or start_blocking_len == 0: |
---|
| 290 | self._A_can_be_reused = True |
---|
| 291 | z = self.interpolate_block(f, point_coordinates, |
---|
[7673] | 292 | verbose=verbose, output_centroids=output_centroids) |
---|
[6073] | 293 | else: |
---|
[6184] | 294 | # Handle blocking |
---|
[6073] | 295 | self._A_can_be_reused = False |
---|
| 296 | start = 0 |
---|
| 297 | # creating a dummy array to concatenate to. |
---|
[7276] | 298 | |
---|
| 299 | f = ensure_numeric(f, num.float) |
---|
[6073] | 300 | if len(f.shape) > 1: |
---|
[7276] | 301 | z = num.zeros((0, f.shape[1]), num.int) #array default# |
---|
[6073] | 302 | else: |
---|
[7276] | 303 | z = num.zeros((0,), num.int) #array default# |
---|
| 304 | |
---|
[6073] | 305 | for end in range(start_blocking_len, |
---|
| 306 | len(point_coordinates), |
---|
| 307 | start_blocking_len): |
---|
| 308 | t = self.interpolate_block(f, point_coordinates[start:end], |
---|
[7673] | 309 | verbose=verbose, output_centroids=output_centroids) |
---|
[7276] | 310 | z = num.concatenate((z, t), axis=0) #??default# |
---|
[6073] | 311 | start = end |
---|
| 312 | |
---|
| 313 | end = len(point_coordinates) |
---|
| 314 | t = self.interpolate_block(f, point_coordinates[start:end], |
---|
[7673] | 315 | verbose=verbose, output_centroids=output_centroids) |
---|
[7276] | 316 | z = num.concatenate((z, t), axis=0) #??default# |
---|
[6073] | 317 | return z |
---|
| 318 | |
---|
[7276] | 319 | |
---|
[6073] | 320 | ## |
---|
[7673] | 321 | # @brief Interpolate a block of vertices |
---|
| 322 | # @param f Array of arbitrary data to be interpolated |
---|
| 323 | # @param point_coordinates List of vertices to intersect with the mesh |
---|
| 324 | # @param use_cache True if caching should be used to accelerate the calculations |
---|
[6073] | 325 | # @param verbose True if this function is verbose. |
---|
| 326 | # @return ?? |
---|
[7276] | 327 | def interpolate_block(self, f, point_coordinates, |
---|
[7673] | 328 | use_cache=False, verbose=False, output_centroids=False): |
---|
[6073] | 329 | """ |
---|
| 330 | Call this if you want to control the blocking or make sure blocking |
---|
| 331 | doesn't occur. |
---|
| 332 | |
---|
| 333 | Return the point data, z. |
---|
[7276] | 334 | |
---|
[6073] | 335 | See interpolate for doc info. |
---|
| 336 | """ |
---|
[7276] | 337 | |
---|
[6073] | 338 | # FIXME (Ole): I reckon we should change the interface so that |
---|
[7276] | 339 | # the user can specify the interpolation matrix instead of the |
---|
[6073] | 340 | # interpolation points to save time. |
---|
| 341 | |
---|
| 342 | if isinstance(point_coordinates, Geospatial_data): |
---|
| 343 | point_coordinates = point_coordinates.get_data_points(absolute=True) |
---|
| 344 | |
---|
[7276] | 345 | # Convert lists to numeric arrays if necessary |
---|
| 346 | point_coordinates = ensure_numeric(point_coordinates, num.float) |
---|
| 347 | f = ensure_numeric(f, num.float) |
---|
[6073] | 348 | |
---|
| 349 | from anuga.caching import myhash |
---|
[7276] | 350 | import sys |
---|
| 351 | |
---|
[6073] | 352 | if use_cache is True: |
---|
| 353 | if sys.platform != 'win32': |
---|
| 354 | # FIXME (Ole): (Why doesn't this work on windoze?) |
---|
| 355 | # Still absolutely fails on Win 24 Oct 2008 |
---|
[7276] | 356 | |
---|
[6073] | 357 | X = cache(self._build_interpolation_matrix_A, |
---|
[7673] | 358 | args=(point_coordinates, output_centroids), |
---|
[7276] | 359 | kwargs={'verbose': verbose}, |
---|
| 360 | verbose=verbose) |
---|
[6073] | 361 | else: |
---|
| 362 | # FIXME |
---|
| 363 | # Hash point_coordinates to memory location, reuse if possible |
---|
| 364 | # This will work on Linux as well if we want to use it there. |
---|
| 365 | key = myhash(point_coordinates) |
---|
[7276] | 366 | |
---|
| 367 | reuse_A = False |
---|
| 368 | |
---|
[6073] | 369 | if self.interpolation_matrices.has_key(key): |
---|
[7276] | 370 | X, stored_points = self.interpolation_matrices[key] |
---|
[6152] | 371 | if num.alltrue(stored_points == point_coordinates): |
---|
[7276] | 372 | reuse_A = True # Reuse interpolation matrix |
---|
| 373 | |
---|
[6073] | 374 | if reuse_A is False: |
---|
| 375 | X = self._build_interpolation_matrix_A(point_coordinates, |
---|
[7673] | 376 | output_centroids, |
---|
[6073] | 377 | verbose=verbose) |
---|
| 378 | self.interpolation_matrices[key] = (X, point_coordinates) |
---|
| 379 | else: |
---|
[7673] | 380 | X = self._build_interpolation_matrix_A(point_coordinates, output_centroids, |
---|
[7276] | 381 | verbose=verbose) |
---|
[6073] | 382 | |
---|
[7276] | 383 | # Unpack result |
---|
[7675] | 384 | self._A, self.inside_poly_indices, self.outside_poly_indices, self.centroids = X |
---|
[6073] | 385 | |
---|
| 386 | # Check that input dimensions are compatible |
---|
| 387 | msg = 'Two columns must be specified in point coordinates. ' \ |
---|
| 388 | 'I got shape=%s' % (str(point_coordinates.shape)) |
---|
| 389 | assert point_coordinates.shape[1] == 2, msg |
---|
| 390 | |
---|
| 391 | msg = 'The number of rows in matrix A must be the same as the ' |
---|
| 392 | msg += 'number of points supplied.' |
---|
| 393 | msg += ' I got %d points and %d matrix rows.' \ |
---|
| 394 | % (point_coordinates.shape[0], self._A.shape[0]) |
---|
[7276] | 395 | assert point_coordinates.shape[0] == self._A.shape[0], msg |
---|
[6073] | 396 | |
---|
| 397 | msg = 'The number of columns in matrix A must be the same as the ' |
---|
| 398 | msg += 'number of mesh vertices.' |
---|
| 399 | msg += ' I got %d vertices and %d matrix columns.' \ |
---|
[7276] | 400 | % (f.shape[0], self._A.shape[1]) |
---|
[6073] | 401 | assert self._A.shape[1] == f.shape[0], msg |
---|
| 402 | |
---|
| 403 | # Compute Matrix vector product and return |
---|
| 404 | return self._get_point_data_z(f) |
---|
| 405 | |
---|
| 406 | |
---|
| 407 | ## |
---|
[7673] | 408 | # @brief Get interpolated data at given points. |
---|
| 409 | # Applies a transform to all points to calculate the |
---|
| 410 | # interpolated values. Points outside the mesh are returned as NaN. |
---|
| 411 | # @note self._A matrix must be valid |
---|
| 412 | # @param f Array of arbitrary data |
---|
[6073] | 413 | # @param verbose True if this function is to be verbose. |
---|
[7673] | 414 | # @return f transformed by interpolation matrix (f') |
---|
[6073] | 415 | def _get_point_data_z(self, f, verbose=False): |
---|
| 416 | """ |
---|
| 417 | Return the point data, z. |
---|
[7276] | 418 | |
---|
[6073] | 419 | Precondition: The _A matrix has been created |
---|
| 420 | """ |
---|
| 421 | |
---|
| 422 | z = self._A * f |
---|
| 423 | |
---|
| 424 | # Taking into account points outside the mesh. |
---|
[7276] | 425 | for i in self.outside_poly_indices: |
---|
[6073] | 426 | z[i] = NAN |
---|
| 427 | return z |
---|
| 428 | |
---|
| 429 | |
---|
| 430 | ## |
---|
| 431 | # @brief Build NxM interpolation matrix. |
---|
[7673] | 432 | # @param point_coordinates Points to sample at |
---|
| 433 | # @param output_centroids set to True to always sample from the centre |
---|
| 434 | # of the intersected triangle, instead of the intersection |
---|
| 435 | # point. |
---|
[6073] | 436 | # @param verbose True if this function is to be verbose. |
---|
[7675] | 437 | # @return Interpolation matrix A, plus lists of the points inside and outside the mesh |
---|
| 438 | # and the list of centroids, if requested. |
---|
[6073] | 439 | def _build_interpolation_matrix_A(self, |
---|
| 440 | point_coordinates, |
---|
[7673] | 441 | output_centroids=False, |
---|
[6073] | 442 | verbose=False): |
---|
| 443 | """Build n x m interpolation matrix, where |
---|
| 444 | n is the number of data points and |
---|
| 445 | m is the number of basis functions phi_k (one per vertex) |
---|
| 446 | |
---|
| 447 | This algorithm uses a quad tree data structure for fast binning |
---|
| 448 | of data points |
---|
| 449 | origin is a 3-tuple consisting of UTM zone, easting and northing. |
---|
| 450 | If specified coordinates are assumed to be relative to this origin. |
---|
| 451 | |
---|
| 452 | This one will override any data_origin that may be specified in |
---|
| 453 | instance interpolation |
---|
| 454 | |
---|
| 455 | Preconditions: |
---|
| 456 | Point_coordindates and mesh vertices have the same origin. |
---|
| 457 | """ |
---|
| 458 | |
---|
[7317] | 459 | if verbose: log.critical('Building interpolation matrix') |
---|
[6073] | 460 | |
---|
[7276] | 461 | # Convert point_coordinates to numeric arrays, in case it was a list. |
---|
| 462 | point_coordinates = ensure_numeric(point_coordinates, num.float) |
---|
[6073] | 463 | |
---|
[7317] | 464 | if verbose: log.critical('Getting indices inside mesh boundary') |
---|
[6073] | 465 | |
---|
| 466 | inside_poly_indices, outside_poly_indices = \ |
---|
| 467 | in_and_outside_polygon(point_coordinates, |
---|
| 468 | self.mesh.get_boundary_polygon(), |
---|
| 469 | closed=True, verbose=verbose) |
---|
[7276] | 470 | |
---|
[6073] | 471 | # Build n x m interpolation matrix |
---|
| 472 | if verbose and len(outside_poly_indices) > 0: |
---|
[7317] | 473 | log.critical('WARNING: Points outside mesh boundary.') |
---|
[7276] | 474 | |
---|
[6073] | 475 | # Since you can block, throw a warning, not an error. |
---|
| 476 | if verbose and 0 == len(inside_poly_indices): |
---|
[7317] | 477 | log.critical('WARNING: No points within the mesh!') |
---|
[7276] | 478 | |
---|
[6073] | 479 | m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex) |
---|
| 480 | n = point_coordinates.shape[0] # Nbr of data points |
---|
| 481 | |
---|
[7317] | 482 | if verbose: log.critical('Number of datapoints: %d' % n) |
---|
| 483 | if verbose: log.critical('Number of basis functions: %d' % m) |
---|
[6073] | 484 | |
---|
| 485 | A = Sparse(n,m) |
---|
| 486 | |
---|
| 487 | n = len(inside_poly_indices) |
---|
[7276] | 488 | |
---|
[7675] | 489 | centroids = [] |
---|
| 490 | |
---|
[6073] | 491 | # Compute matrix elements for points inside the mesh |
---|
[7317] | 492 | if verbose: log.critical('Building interpolation matrix from %d points' |
---|
| 493 | % n) |
---|
[6073] | 494 | |
---|
| 495 | for d, i in enumerate(inside_poly_indices): |
---|
| 496 | # For each data_coordinate point |
---|
[7317] | 497 | if verbose and d%((n+10)/10)==0: log.critical('Doing %d of %d' |
---|
| 498 | %(d, n)) |
---|
[6073] | 499 | |
---|
| 500 | x = point_coordinates[i] |
---|
| 501 | element_found, sigma0, sigma1, sigma2, k = \ |
---|
[6545] | 502 | search_tree_of_vertices(self.root, x) |
---|
[6073] | 503 | |
---|
| 504 | # Update interpolation matrix A if necessary |
---|
| 505 | if element_found is True: |
---|
| 506 | # Assign values to matrix A |
---|
| 507 | j0 = self.mesh.triangles[k,0] # Global vertex id for sigma0 |
---|
| 508 | j1 = self.mesh.triangles[k,1] # Global vertex id for sigma1 |
---|
| 509 | j2 = self.mesh.triangles[k,2] # Global vertex id for sigma2 |
---|
[7673] | 510 | js = [j0, j1, j2] |
---|
| 511 | |
---|
| 512 | if output_centroids is False: |
---|
| 513 | # Weight each vertex according to its distance from x |
---|
| 514 | sigmas = {j0:sigma0, j1:sigma1, j2:sigma2} |
---|
| 515 | for j in js: |
---|
| 516 | A[i, j] = sigmas[j] |
---|
| 517 | else: |
---|
| 518 | # If centroids are needed, weight all 3 vertices equally |
---|
| 519 | for j in js: |
---|
[7675] | 520 | A[i, j] = 1.0/3.0 |
---|
| 521 | centroids.append(self.mesh.centroid_coordinates[k]) |
---|
[6073] | 522 | else: |
---|
[7276] | 523 | msg = 'Could not find triangle for point', x |
---|
[6073] | 524 | raise Exception(msg) |
---|
[7675] | 525 | return A, inside_poly_indices, outside_poly_indices, centroids |
---|
[6073] | 526 | |
---|
| 527 | |
---|
[7276] | 528 | |
---|
| 529 | |
---|
| 530 | |
---|
| 531 | |
---|
[6073] | 532 | ## |
---|
| 533 | # @brief ?? |
---|
| 534 | # @param vertices ?? |
---|
| 535 | # @param vertex_attributes ?? |
---|
| 536 | # @param triangles ?? |
---|
| 537 | # @param points ?? |
---|
| 538 | # @param max_points_per_cell ?? |
---|
| 539 | # @param start_blocking_len ?? |
---|
| 540 | # @param mesh_origin ?? |
---|
| 541 | def benchmark_interpolate(vertices, |
---|
| 542 | vertex_attributes, |
---|
| 543 | triangles, points, |
---|
| 544 | max_points_per_cell=None, |
---|
| 545 | start_blocking_len=500000, |
---|
| 546 | mesh_origin=None): |
---|
| 547 | """ |
---|
| 548 | points: Interpolate mesh data to these positions. |
---|
| 549 | List of coordinate pairs [x, y] of |
---|
[7276] | 550 | data points or an nx2 numeric array or a Geospatial_data object |
---|
| 551 | |
---|
[6073] | 552 | No test for this yet. |
---|
| 553 | Note, this has no time the input data has no time dimension. Which is |
---|
| 554 | different from most of the data we interpolate, eg sww info. |
---|
[7276] | 555 | |
---|
[6073] | 556 | Output: |
---|
| 557 | Interpolated values at inputted points. |
---|
| 558 | """ |
---|
| 559 | |
---|
| 560 | interp = Interpolate(vertices, |
---|
[7276] | 561 | triangles, |
---|
[6073] | 562 | max_vertices_per_cell=max_points_per_cell, |
---|
| 563 | mesh_origin=mesh_origin) |
---|
[7276] | 564 | |
---|
[6073] | 565 | calc = interp.interpolate(vertex_attributes, |
---|
| 566 | points, |
---|
| 567 | start_blocking_len=start_blocking_len) |
---|
| 568 | |
---|
| 569 | |
---|
| 570 | ## |
---|
| 571 | # @brief Interpolate quantities at given locations (from .SWW file). |
---|
| 572 | # @param sww_file Input .SWW file. |
---|
| 573 | # @param points A list of the 'gauges' x,y location. |
---|
| 574 | # @param depth_file The name of the output depth file. |
---|
| 575 | # @param velocity_x_file Name of the output x velocity file. |
---|
| 576 | # @param velocity_y_file Name of the output y velocity file. |
---|
| 577 | # @param stage_file Name of the output stage file. |
---|
[7276] | 578 | # @param froude_file |
---|
[6073] | 579 | # @param time_thinning Time thinning step to use. |
---|
| 580 | # @param verbose True if this function is to be verbose. |
---|
| 581 | # @param use_cache True if we are caching. |
---|
| 582 | def interpolate_sww2csv(sww_file, |
---|
| 583 | points, |
---|
| 584 | depth_file, |
---|
| 585 | velocity_x_file, |
---|
| 586 | velocity_y_file, |
---|
| 587 | stage_file=None, |
---|
| 588 | froude_file=None, |
---|
| 589 | time_thinning=1, |
---|
| 590 | verbose=True, |
---|
| 591 | use_cache = True): |
---|
| 592 | """ |
---|
| 593 | Interpolate the quantities at a given set of locations, given |
---|
| 594 | an sww file. |
---|
| 595 | The results are written to csv files. |
---|
| 596 | |
---|
| 597 | sww_file is the input sww file. |
---|
| 598 | points is a list of the 'gauges' x,y location. |
---|
| 599 | depth_file is the name of the output depth file |
---|
| 600 | velocity_x_file is the name of the output x velocity file. |
---|
| 601 | velocity_y_file is the name of the output y velocity file. |
---|
| 602 | stage_file is the name of the output stage file. |
---|
| 603 | |
---|
| 604 | In the csv files columns represents the gauges and each row is a |
---|
| 605 | time slice. |
---|
| 606 | |
---|
| 607 | Time_thinning_number controls how many timesteps to use. Only |
---|
| 608 | timesteps with index%time_thinning_number == 0 will used, or |
---|
| 609 | in other words a value of 3, say, will cause the algorithm to |
---|
| 610 | use every third time step. |
---|
| 611 | |
---|
| 612 | In the future let points be a points file. |
---|
| 613 | And let the user choose the quantities. |
---|
| 614 | |
---|
| 615 | This is currently quite specific. |
---|
| 616 | If it is need to be more general, change things. |
---|
| 617 | """ |
---|
| 618 | |
---|
| 619 | quantities = ['stage', 'elevation', 'xmomentum', 'ymomentum'] |
---|
| 620 | points = ensure_absolute(points) |
---|
| 621 | point_count = len(points) |
---|
| 622 | callable_sww = file_function(sww_file, |
---|
| 623 | quantities=quantities, |
---|
| 624 | interpolation_points=points, |
---|
| 625 | verbose=verbose, |
---|
| 626 | time_thinning=time_thinning, |
---|
| 627 | use_cache=use_cache) |
---|
[7276] | 628 | |
---|
[6073] | 629 | depth_writer = writer(file(depth_file, "wb")) |
---|
| 630 | velocity_x_writer = writer(file(velocity_x_file, "wb")) |
---|
| 631 | velocity_y_writer = writer(file(velocity_y_file, "wb")) |
---|
| 632 | if stage_file is not None: |
---|
| 633 | stage_writer = writer(file(stage_file, "wb")) |
---|
| 634 | if froude_file is not None: |
---|
| 635 | froude_writer = writer(file(froude_file, "wb")) |
---|
| 636 | |
---|
| 637 | # Write heading |
---|
| 638 | heading = [str(x[0])+ ':' + str(x[1]) for x in points] |
---|
| 639 | heading.insert(0, "time") |
---|
| 640 | depth_writer.writerow(heading) |
---|
| 641 | velocity_x_writer.writerow(heading) |
---|
| 642 | velocity_y_writer.writerow(heading) |
---|
| 643 | if stage_file is not None: |
---|
[7276] | 644 | stage_writer.writerow(heading) |
---|
[6073] | 645 | if froude_file is not None: |
---|
[7276] | 646 | froude_writer.writerow(heading) |
---|
| 647 | |
---|
[6073] | 648 | for time in callable_sww.get_time(): |
---|
| 649 | depths = [time] |
---|
| 650 | velocity_xs = [time] |
---|
| 651 | velocity_ys = [time] |
---|
[7276] | 652 | if stage_file is not None: |
---|
| 653 | stages = [time] |
---|
| 654 | if froude_file is not None: |
---|
| 655 | froudes = [time] |
---|
[6073] | 656 | for point_i, point in enumerate(points): |
---|
| 657 | quantities = callable_sww(time,point_i) |
---|
[7276] | 658 | |
---|
[6073] | 659 | w = quantities[0] |
---|
| 660 | z = quantities[1] |
---|
| 661 | momentum_x = quantities[2] |
---|
| 662 | momentum_y = quantities[3] |
---|
[7276] | 663 | depth = w - z |
---|
| 664 | |
---|
[6073] | 665 | if w == NAN or z == NAN or momentum_x == NAN: |
---|
| 666 | velocity_x = NAN |
---|
| 667 | else: |
---|
| 668 | if depth > 1.e-30: # use epsilon |
---|
| 669 | velocity_x = momentum_x / depth #Absolute velocity |
---|
| 670 | else: |
---|
| 671 | velocity_x = 0 |
---|
| 672 | |
---|
| 673 | if w == NAN or z == NAN or momentum_y == NAN: |
---|
| 674 | velocity_y = NAN |
---|
| 675 | else: |
---|
| 676 | if depth > 1.e-30: # use epsilon |
---|
| 677 | velocity_y = momentum_y / depth #Absolute velocity |
---|
| 678 | else: |
---|
| 679 | velocity_y = 0 |
---|
| 680 | |
---|
| 681 | if depth < 1.e-30: # use epsilon |
---|
| 682 | froude = NAN |
---|
| 683 | else: |
---|
| 684 | froude = sqrt(velocity_x*velocity_x + velocity_y*velocity_y)/ \ |
---|
| 685 | sqrt(depth * 9.8066) # gravity m/s/s |
---|
| 686 | |
---|
| 687 | depths.append(depth) |
---|
| 688 | velocity_xs.append(velocity_x) |
---|
| 689 | velocity_ys.append(velocity_y) |
---|
| 690 | |
---|
| 691 | if stage_file is not None: |
---|
| 692 | stages.append(w) |
---|
| 693 | if froude_file is not None: |
---|
| 694 | froudes.append(froude) |
---|
| 695 | |
---|
| 696 | depth_writer.writerow(depths) |
---|
| 697 | velocity_x_writer.writerow(velocity_xs) |
---|
| 698 | velocity_y_writer.writerow(velocity_ys) |
---|
| 699 | |
---|
| 700 | if stage_file is not None: |
---|
[7276] | 701 | stage_writer.writerow(stages) |
---|
[6073] | 702 | if froude_file is not None: |
---|
[7276] | 703 | froude_writer.writerow(froudes) |
---|
[6073] | 704 | |
---|
| 705 | |
---|
| 706 | ## |
---|
[7276] | 707 | # @brief |
---|
[6073] | 708 | class Interpolation_function: |
---|
| 709 | """Interpolation_interface - creates callable object f(t, id) or f(t,x,y) |
---|
| 710 | which is interpolated from time series defined at vertices of |
---|
| 711 | triangular mesh (such as those stored in sww files) |
---|
| 712 | |
---|
| 713 | Let m be the number of vertices, n the number of triangles |
---|
| 714 | and p the number of timesteps. |
---|
| 715 | Also, let N be the number of interpolation points. |
---|
| 716 | |
---|
| 717 | Mandatory input |
---|
[7276] | 718 | time: px1 array of monotonously increasing times (float) |
---|
| 719 | quantities: Dictionary of arrays or 1 array (float) |
---|
[6073] | 720 | The arrays must either have dimensions pxm or mx1. |
---|
| 721 | The resulting function will be time dependent in |
---|
| 722 | the former case while it will be constant with |
---|
| 723 | respect to time in the latter case. |
---|
[7276] | 724 | |
---|
[6073] | 725 | Optional input: |
---|
| 726 | quantity_names: List of keys into the quantities dictionary for |
---|
| 727 | imposing a particular order on the output vector. |
---|
[7276] | 728 | vertex_coordinates: mx2 array of coordinates (float) |
---|
| 729 | triangles: nx3 array of indices into vertex_coordinates (int) |
---|
| 730 | interpolation_points: Nx2 array of coordinates to be interpolated to |
---|
[6073] | 731 | verbose: Level of reporting |
---|
| 732 | |
---|
| 733 | The quantities returned by the callable object are specified by |
---|
| 734 | the list quantities which must contain the names of the |
---|
| 735 | quantities to be returned and also reflect the order, e.g. for |
---|
| 736 | the shallow water wave equation, on would have |
---|
| 737 | quantities = ['stage', 'xmomentum', 'ymomentum'] |
---|
| 738 | |
---|
| 739 | The parameter interpolation_points decides at which points interpolated |
---|
| 740 | quantities are to be computed whenever object is called. |
---|
| 741 | If None, return average value |
---|
| 742 | |
---|
| 743 | FIXME (Ole): Need to allow vertex coordinates and interpolation points to |
---|
| 744 | be geospatial data objects |
---|
| 745 | |
---|
| 746 | (FIXME (Ole): This comment should be removed) |
---|
| 747 | Time assumed to be relative to starttime |
---|
| 748 | All coordinates assume origin of (0,0) - e.g. georeferencing must be |
---|
| 749 | taken care of outside this function |
---|
| 750 | """ |
---|
| 751 | |
---|
| 752 | ## |
---|
| 753 | # @brief ?? |
---|
| 754 | # @param time ?? |
---|
| 755 | # @param quantities ?? |
---|
| 756 | # @param quantity_names ?? |
---|
| 757 | # @param vertex_coordinates ?? |
---|
| 758 | # @param triangles ?? |
---|
| 759 | # @param interpolation_points ?? |
---|
| 760 | # @param time_thinning ?? |
---|
| 761 | # @param verbose ?? |
---|
| 762 | # @param gauge_neighbour_id ?? |
---|
| 763 | def __init__(self, |
---|
| 764 | time, |
---|
| 765 | quantities, |
---|
[7276] | 766 | quantity_names=None, |
---|
[6073] | 767 | vertex_coordinates=None, |
---|
| 768 | triangles=None, |
---|
| 769 | interpolation_points=None, |
---|
| 770 | time_thinning=1, |
---|
| 771 | verbose=False, |
---|
[7673] | 772 | gauge_neighbour_id=None, |
---|
| 773 | output_centroids=False): |
---|
[6073] | 774 | """Initialise object and build spatial interpolation if required |
---|
| 775 | |
---|
| 776 | Time_thinning_number controls how many timesteps to use. Only timesteps |
---|
| 777 | with index%time_thinning_number == 0 will used, or in other words a |
---|
| 778 | value of 3, say, will cause the algorithm to use every third time step. |
---|
| 779 | """ |
---|
| 780 | |
---|
| 781 | from anuga.config import time_format |
---|
| 782 | import types |
---|
| 783 | |
---|
[6194] | 784 | if verbose is True: |
---|
[7317] | 785 | log.critical('Interpolation_function: input checks') |
---|
[6194] | 786 | |
---|
[7276] | 787 | # Check temporal info |
---|
[6194] | 788 | time = ensure_numeric(time) |
---|
[7326] | 789 | |
---|
[6194] | 790 | if not num.alltrue(time[1:] - time[:-1] >= 0): |
---|
| 791 | # This message is time consuming to form due to the conversion of |
---|
| 792 | msg = 'Time must be a monotonuosly increasing sequence %s' % time |
---|
| 793 | raise Exception, msg |
---|
[6073] | 794 | |
---|
| 795 | # Check if quantities is a single array only |
---|
| 796 | if type(quantities) != types.DictType: |
---|
| 797 | quantities = ensure_numeric(quantities) |
---|
| 798 | quantity_names = ['Attribute'] |
---|
| 799 | |
---|
| 800 | # Make it a dictionary |
---|
| 801 | quantities = {quantity_names[0]: quantities} |
---|
| 802 | |
---|
| 803 | # Use keys if no names are specified |
---|
| 804 | if quantity_names is None: |
---|
| 805 | quantity_names = quantities.keys() |
---|
| 806 | |
---|
| 807 | # Check spatial info |
---|
| 808 | if vertex_coordinates is None: |
---|
| 809 | self.spatial = False |
---|
| 810 | else: |
---|
| 811 | # FIXME (Ole): Try ensure_numeric here - |
---|
| 812 | # this function knows nothing about georefering. |
---|
| 813 | vertex_coordinates = ensure_absolute(vertex_coordinates) |
---|
| 814 | |
---|
| 815 | if triangles is not None: |
---|
| 816 | triangles = ensure_numeric(triangles) |
---|
[7276] | 817 | self.spatial = True |
---|
[6073] | 818 | |
---|
[6194] | 819 | if verbose is True: |
---|
[7317] | 820 | log.critical('Interpolation_function: thinning by %d' |
---|
| 821 | % time_thinning) |
---|
[6194] | 822 | |
---|
[7276] | 823 | |
---|
[6073] | 824 | # Thin timesteps if needed |
---|
| 825 | # Note array() is used to make the thinned arrays contiguous in memory |
---|
[7276] | 826 | self.time = num.array(time[::time_thinning]) |
---|
[6073] | 827 | for name in quantity_names: |
---|
| 828 | if len(quantities[name].shape) == 2: |
---|
[6152] | 829 | quantities[name] = num.array(quantities[name][::time_thinning,:]) |
---|
[6194] | 830 | |
---|
| 831 | if verbose is True: |
---|
[7317] | 832 | log.critical('Interpolation_function: precomputing') |
---|
[7276] | 833 | |
---|
[6073] | 834 | # Save for use with statistics |
---|
| 835 | self.quantities_range = {} |
---|
| 836 | for name in quantity_names: |
---|
[7276] | 837 | q = quantities[name][:].flatten() |
---|
[6073] | 838 | self.quantities_range[name] = [min(q), max(q)] |
---|
[7276] | 839 | |
---|
| 840 | self.quantity_names = quantity_names |
---|
| 841 | self.vertex_coordinates = vertex_coordinates |
---|
[6073] | 842 | self.interpolation_points = interpolation_points |
---|
| 843 | |
---|
| 844 | self.index = 0 # Initial time index |
---|
| 845 | self.precomputed_values = {} |
---|
[7675] | 846 | self.centroids = [] |
---|
[6073] | 847 | |
---|
| 848 | # Precomputed spatial interpolation if requested |
---|
| 849 | if interpolation_points is not None: |
---|
| 850 | #no longer true. sts files have spatial = True but |
---|
| 851 | #if self.spatial is False: |
---|
| 852 | # raise 'Triangles and vertex_coordinates must be specified' |
---|
[7276] | 853 | # |
---|
[6073] | 854 | try: |
---|
[7276] | 855 | self.interpolation_points = \ |
---|
[6073] | 856 | interpolation_points = ensure_numeric(interpolation_points) |
---|
| 857 | except: |
---|
[7276] | 858 | msg = 'Interpolation points must be an N x 2 numeric array ' \ |
---|
[6073] | 859 | 'or a list of points\n' |
---|
[7276] | 860 | msg += 'Got: %s.' %(str(self.interpolation_points)[:60] + '...') |
---|
[6073] | 861 | raise msg |
---|
| 862 | |
---|
[6621] | 863 | # Ensure 'mesh_boundary_polygon' is defined |
---|
| 864 | mesh_boundary_polygon = None |
---|
| 865 | |
---|
[6073] | 866 | if triangles is not None and vertex_coordinates is not None: |
---|
| 867 | # Check that all interpolation points fall within |
---|
| 868 | # mesh boundary as defined by triangles and vertex_coordinates. |
---|
| 869 | from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh |
---|
[7276] | 870 | from anuga.utilities.polygon import outside_polygon |
---|
[6073] | 871 | |
---|
| 872 | # Create temporary mesh object from mesh info passed |
---|
[7276] | 873 | # into this function. |
---|
[6073] | 874 | mesh = Mesh(vertex_coordinates, triangles) |
---|
| 875 | mesh_boundary_polygon = mesh.get_boundary_polygon() |
---|
| 876 | |
---|
| 877 | indices = outside_polygon(interpolation_points, |
---|
| 878 | mesh_boundary_polygon) |
---|
| 879 | |
---|
| 880 | # Record result |
---|
| 881 | #self.mesh_boundary_polygon = mesh_boundary_polygon |
---|
| 882 | self.indices_outside_mesh = indices |
---|
| 883 | |
---|
| 884 | # Report |
---|
| 885 | if len(indices) > 0: |
---|
| 886 | msg = 'Interpolation points in Interpolation function fall ' |
---|
| 887 | msg += 'outside specified mesh. Offending points:\n' |
---|
| 888 | out_interp_pts = [] |
---|
| 889 | for i in indices: |
---|
| 890 | msg += '%d: %s\n' % (i, interpolation_points[i]) |
---|
| 891 | out_interp_pts.append( |
---|
| 892 | ensure_numeric(interpolation_points[i])) |
---|
| 893 | |
---|
| 894 | if verbose is True: |
---|
| 895 | import sys |
---|
| 896 | if sys.platform == 'win32': |
---|
| 897 | # FIXME (Ole): Why only Windoze? |
---|
| 898 | from anuga.utilities.polygon import plot_polygons |
---|
| 899 | title = ('Interpolation points fall ' |
---|
| 900 | 'outside specified mesh') |
---|
| 901 | plot_polygons([mesh_boundary_polygon, |
---|
| 902 | interpolation_points, |
---|
| 903 | out_interp_pts], |
---|
| 904 | ['line', 'point', 'outside'], |
---|
| 905 | figname='points_boundary_out', |
---|
| 906 | label=title, |
---|
| 907 | verbose=verbose) |
---|
| 908 | |
---|
| 909 | # Joaquim Luis suggested this as an Exception, so |
---|
| 910 | # that the user can now what the problem is rather than |
---|
| 911 | # looking for NaN's. However, NANs are handy as they can |
---|
| 912 | # be ignored leaving good points for continued processing. |
---|
| 913 | if verbose: |
---|
[7317] | 914 | log.critical(msg) |
---|
[6073] | 915 | #raise Exception(msg) |
---|
[7276] | 916 | |
---|
[6073] | 917 | elif triangles is None and vertex_coordinates is not None: #jj |
---|
| 918 | #Dealing with sts file |
---|
| 919 | pass |
---|
| 920 | else: |
---|
| 921 | raise Exception('Sww file function requires both triangles and ' |
---|
| 922 | 'vertex_coordinates. sts file file function ' |
---|
| 923 | 'requires the latter.') |
---|
| 924 | |
---|
[6621] | 925 | # Plot boundary and interpolation points, |
---|
| 926 | # but only if if 'mesh_boundary_polygon' has data. |
---|
| 927 | if verbose is True and mesh_boundary_polygon is not None: |
---|
[6073] | 928 | import sys |
---|
| 929 | if sys.platform == 'win32': |
---|
| 930 | from anuga.utilities.polygon import plot_polygons |
---|
| 931 | title = ('Interpolation function: ' |
---|
| 932 | 'Polygon and interpolation points') |
---|
| 933 | plot_polygons([mesh_boundary_polygon, |
---|
| 934 | interpolation_points], |
---|
| 935 | ['line', 'point'], |
---|
| 936 | figname='points_boundary', |
---|
| 937 | label=title, |
---|
| 938 | verbose=verbose) |
---|
| 939 | |
---|
| 940 | m = len(self.interpolation_points) |
---|
| 941 | p = len(self.time) |
---|
| 942 | |
---|
[7276] | 943 | for name in quantity_names: |
---|
| 944 | self.precomputed_values[name] = num.zeros((p, m), num.float) |
---|
| 945 | |
---|
[6223] | 946 | if verbose is True: |
---|
[7317] | 947 | log.critical('Build interpolator') |
---|
[6223] | 948 | |
---|
[7276] | 949 | |
---|
[6073] | 950 | # Build interpolator |
---|
[6223] | 951 | if triangles is not None and vertex_coordinates is not None: |
---|
[7276] | 952 | if verbose: |
---|
[6073] | 953 | msg = 'Building interpolation matrix from source mesh ' |
---|
| 954 | msg += '(%d vertices, %d triangles)' \ |
---|
| 955 | % (vertex_coordinates.shape[0], |
---|
| 956 | triangles.shape[0]) |
---|
[7317] | 957 | log.critical(msg) |
---|
[7276] | 958 | |
---|
[6223] | 959 | # This one is no longer needed for STS files |
---|
| 960 | interpol = Interpolate(vertex_coordinates, |
---|
| 961 | triangles, |
---|
[7276] | 962 | verbose=verbose) |
---|
| 963 | |
---|
[6223] | 964 | elif triangles is None and vertex_coordinates is not None: |
---|
| 965 | if verbose: |
---|
[7317] | 966 | log.critical('Interpolation from STS file') |
---|
[6073] | 967 | |
---|
| 968 | |
---|
[6223] | 969 | |
---|
[6073] | 970 | if verbose: |
---|
[7317] | 971 | log.critical('Interpolating (%d interpolation points, %d timesteps).' |
---|
| 972 | % (self.interpolation_points.shape[0], self.time.shape[0])) |
---|
[7276] | 973 | |
---|
[6073] | 974 | if time_thinning > 1: |
---|
[7317] | 975 | log.critical('Timesteps were thinned by a factor of %d' |
---|
| 976 | % time_thinning) |
---|
[6073] | 977 | else: |
---|
[7317] | 978 | log.critical() |
---|
[6073] | 979 | |
---|
[7276] | 980 | for i, t in enumerate(self.time): |
---|
[6073] | 981 | # Interpolate quantities at this timestep |
---|
| 982 | if verbose and i%((p+10)/10) == 0: |
---|
[7317] | 983 | log.critical(' time step %d of %d' % (i, p)) |
---|
[7276] | 984 | |
---|
[6073] | 985 | for name in quantity_names: |
---|
| 986 | if len(quantities[name].shape) == 2: |
---|
| 987 | Q = quantities[name][i,:] # Quantities at timestep i |
---|
| 988 | else: |
---|
| 989 | Q = quantities[name][:] # No time dependency |
---|
| 990 | |
---|
| 991 | if verbose and i%((p+10)/10) == 0: |
---|
[7317] | 992 | log.critical(' quantity %s, size=%d' % (name, len(Q))) |
---|
[7276] | 993 | |
---|
| 994 | # Interpolate |
---|
[6073] | 995 | if triangles is not None and vertex_coordinates is not None: |
---|
| 996 | result = interpol.interpolate(Q, |
---|
| 997 | point_coordinates=\ |
---|
| 998 | self.interpolation_points, |
---|
[7673] | 999 | verbose=False, |
---|
[7675] | 1000 | output_centroids=output_centroids) |
---|
| 1001 | self.centroids = interpol.centroids |
---|
[6073] | 1002 | elif triangles is None and vertex_coordinates is not None: |
---|
[6185] | 1003 | result = interpolate_polyline(Q, |
---|
| 1004 | vertex_coordinates, |
---|
| 1005 | gauge_neighbour_id, |
---|
[6189] | 1006 | interpolation_points=\ |
---|
[6073] | 1007 | self.interpolation_points) |
---|
[7276] | 1008 | |
---|
[6073] | 1009 | #assert len(result), len(interpolation_points) |
---|
[7675] | 1010 | self.precomputed_values[name][i, :] = result |
---|
| 1011 | |
---|
[6073] | 1012 | # Report |
---|
| 1013 | if verbose: |
---|
[7675] | 1014 | log.critical(self.statistics()) |
---|
[6073] | 1015 | else: |
---|
| 1016 | # Store quantitites as is |
---|
| 1017 | for name in quantity_names: |
---|
| 1018 | self.precomputed_values[name] = quantities[name] |
---|
| 1019 | |
---|
| 1020 | ## |
---|
| 1021 | # @brief Override object representation method. |
---|
| 1022 | def __repr__(self): |
---|
| 1023 | # return 'Interpolation function (spatio-temporal)' |
---|
| 1024 | return self.statistics() |
---|
| 1025 | |
---|
| 1026 | ## |
---|
[7673] | 1027 | # @brief Evaluate interpolation function |
---|
[6073] | 1028 | # @param t Model time - must lie within existing timesteps. |
---|
| 1029 | # @param point_id Index of one of the preprocessed points. |
---|
| 1030 | # @param x ?? |
---|
| 1031 | # @param y ?? |
---|
| 1032 | # @return ?? |
---|
| 1033 | def __call__(self, t, point_id=None, x=None, y=None): |
---|
| 1034 | """Evaluate f(t) or f(t, point_id) |
---|
| 1035 | |
---|
[7276] | 1036 | Inputs: |
---|
| 1037 | t: time - Model time. Must lie within existing timesteps |
---|
| 1038 | point_id: index of one of the preprocessed points. |
---|
| 1039 | |
---|
| 1040 | If spatial info is present and all of point_id |
---|
| 1041 | are None an exception is raised |
---|
| 1042 | |
---|
[6073] | 1043 | If no spatial info is present, point_id arguments are ignored |
---|
| 1044 | making f a function of time only. |
---|
| 1045 | |
---|
[7276] | 1046 | FIXME: f(t, x, y) x, y could overrided location, point_id ignored |
---|
| 1047 | FIXME: point_id could also be a slice |
---|
| 1048 | FIXME: What if x and y are vectors? |
---|
[6073] | 1049 | FIXME: What about f(x,y) without t? |
---|
| 1050 | """ |
---|
| 1051 | |
---|
| 1052 | from math import pi, cos, sin, sqrt |
---|
[7276] | 1053 | from anuga.abstract_2d_finite_volumes.util import mean |
---|
[6073] | 1054 | |
---|
| 1055 | if self.spatial is True: |
---|
| 1056 | if point_id is None: |
---|
| 1057 | if x is None or y is None: |
---|
| 1058 | msg = 'Either point_id or x and y must be specified' |
---|
| 1059 | raise Exception(msg) |
---|
| 1060 | else: |
---|
| 1061 | if self.interpolation_points is None: |
---|
| 1062 | msg = 'Interpolation_function must be instantiated ' + \ |
---|
| 1063 | 'with a list of interpolation points before ' + \ |
---|
| 1064 | 'parameter point_id can be used' |
---|
| 1065 | raise Exception(msg) |
---|
| 1066 | |
---|
| 1067 | msg = 'Time interval [%.16f:%.16f]' % (self.time[0], self.time[-1]) |
---|
| 1068 | msg += ' does not match model time: %.16f\n' % t |
---|
| 1069 | if t < self.time[0]: raise Modeltime_too_early(msg) |
---|
| 1070 | if t > self.time[-1]: raise Modeltime_too_late(msg) |
---|
| 1071 | |
---|
| 1072 | oldindex = self.index #Time index |
---|
| 1073 | |
---|
| 1074 | # Find current time slot |
---|
| 1075 | while t > self.time[self.index]: self.index += 1 |
---|
| 1076 | while t < self.time[self.index]: self.index -= 1 |
---|
| 1077 | |
---|
| 1078 | if t == self.time[self.index]: |
---|
| 1079 | # Protect against case where t == T[-1] (last time) |
---|
| 1080 | # - also works in general when t == T[i] |
---|
| 1081 | ratio = 0 |
---|
| 1082 | else: |
---|
| 1083 | # t is now between index and index+1 |
---|
| 1084 | ratio = ((t - self.time[self.index]) / |
---|
| 1085 | (self.time[self.index+1] - self.time[self.index])) |
---|
| 1086 | |
---|
| 1087 | # Compute interpolated values |
---|
[7276] | 1088 | q = num.zeros(len(self.quantity_names), num.float) |
---|
| 1089 | for i, name in enumerate(self.quantity_names): |
---|
[6073] | 1090 | Q = self.precomputed_values[name] |
---|
| 1091 | |
---|
| 1092 | if self.spatial is False: |
---|
[7276] | 1093 | # If there is no spatial info |
---|
[6073] | 1094 | assert len(Q.shape) == 1 |
---|
| 1095 | |
---|
| 1096 | Q0 = Q[self.index] |
---|
| 1097 | if ratio > 0: Q1 = Q[self.index+1] |
---|
| 1098 | else: |
---|
| 1099 | if x is not None and y is not None: |
---|
| 1100 | # Interpolate to x, y |
---|
| 1101 | raise 'x,y interpolation not yet implemented' |
---|
| 1102 | else: |
---|
| 1103 | # Use precomputed point |
---|
| 1104 | Q0 = Q[self.index, point_id] |
---|
| 1105 | if ratio > 0: |
---|
| 1106 | Q1 = Q[self.index+1, point_id] |
---|
| 1107 | |
---|
[7276] | 1108 | # Linear temporal interpolation |
---|
[6073] | 1109 | if ratio > 0: |
---|
| 1110 | if Q0 == NAN and Q1 == NAN: |
---|
| 1111 | q[i] = Q0 |
---|
| 1112 | else: |
---|
| 1113 | q[i] = Q0 + ratio*(Q1 - Q0) |
---|
| 1114 | else: |
---|
| 1115 | q[i] = Q0 |
---|
| 1116 | |
---|
| 1117 | # Return vector of interpolated values |
---|
| 1118 | # FIXME: |
---|
| 1119 | if self.spatial is True: |
---|
| 1120 | return q |
---|
| 1121 | else: |
---|
| 1122 | # Replicate q according to x and y |
---|
| 1123 | # This is e.g used for Wind_stress |
---|
[7276] | 1124 | if x is None or y is None: |
---|
[6073] | 1125 | return q |
---|
| 1126 | else: |
---|
| 1127 | try: |
---|
| 1128 | N = len(x) |
---|
| 1129 | except: |
---|
| 1130 | return q |
---|
| 1131 | else: |
---|
| 1132 | # x is a vector - Create one constant column for each value |
---|
| 1133 | N = len(x) |
---|
| 1134 | assert len(y) == N, 'x and y must have same length' |
---|
| 1135 | res = [] |
---|
| 1136 | for col in q: |
---|
[7276] | 1137 | res.append(col*num.ones(N, num.float)) |
---|
| 1138 | |
---|
[6073] | 1139 | return res |
---|
| 1140 | |
---|
| 1141 | ## |
---|
| 1142 | # @brief Return model time as a vector of timesteps. |
---|
| 1143 | def get_time(self): |
---|
| 1144 | """Return model time as a vector of timesteps |
---|
| 1145 | """ |
---|
| 1146 | return self.time |
---|
| 1147 | |
---|
| 1148 | ## |
---|
| 1149 | # @brief Output statistics about interpolation_function. |
---|
| 1150 | # @return The statistics string. |
---|
| 1151 | def statistics(self): |
---|
| 1152 | """Output statistics about interpolation_function |
---|
| 1153 | """ |
---|
[7276] | 1154 | |
---|
[6073] | 1155 | vertex_coordinates = self.vertex_coordinates |
---|
[7276] | 1156 | interpolation_points = self.interpolation_points |
---|
[6073] | 1157 | quantity_names = self.quantity_names |
---|
| 1158 | #quantities = self.quantities |
---|
[7276] | 1159 | precomputed_values = self.precomputed_values |
---|
| 1160 | |
---|
[6073] | 1161 | x = vertex_coordinates[:,0] |
---|
[7276] | 1162 | y = vertex_coordinates[:,1] |
---|
[6073] | 1163 | |
---|
| 1164 | str = '------------------------------------------------\n' |
---|
| 1165 | str += 'Interpolation_function (spatio-temporal) statistics:\n' |
---|
| 1166 | str += ' Extent:\n' |
---|
| 1167 | str += ' x in [%f, %f], len(x) == %d\n'\ |
---|
| 1168 | %(min(x), max(x), len(x)) |
---|
| 1169 | str += ' y in [%f, %f], len(y) == %d\n'\ |
---|
| 1170 | %(min(y), max(y), len(y)) |
---|
| 1171 | str += ' t in [%f, %f], len(t) == %d\n'\ |
---|
| 1172 | %(min(self.time), max(self.time), len(self.time)) |
---|
| 1173 | str += ' Quantities:\n' |
---|
| 1174 | for name in quantity_names: |
---|
| 1175 | minq, maxq = self.quantities_range[name] |
---|
[7276] | 1176 | str += ' %s in [%f, %f]\n' %(name, minq, maxq) |
---|
| 1177 | #q = quantities[name][:].flatten() |
---|
[6073] | 1178 | #str += ' %s in [%f, %f]\n' %(name, min(q), max(q)) |
---|
| 1179 | |
---|
[7276] | 1180 | if interpolation_points is not None: |
---|
[6073] | 1181 | str += ' Interpolation points (xi, eta):'\ |
---|
| 1182 | ' number of points == %d\n' %interpolation_points.shape[0] |
---|
| 1183 | str += ' xi in [%f, %f]\n' %(min(interpolation_points[:,0]), |
---|
| 1184 | max(interpolation_points[:,0])) |
---|
| 1185 | str += ' eta in [%f, %f]\n' %(min(interpolation_points[:,1]), |
---|
| 1186 | max(interpolation_points[:,1])) |
---|
| 1187 | str += ' Interpolated quantities (over all timesteps):\n' |
---|
[7276] | 1188 | |
---|
[6073] | 1189 | for name in quantity_names: |
---|
[7276] | 1190 | q = precomputed_values[name][:].flatten() |
---|
[6073] | 1191 | str += ' %s at interpolation points in [%f, %f]\n'\ |
---|
| 1192 | %(name, min(q), max(q)) |
---|
| 1193 | str += '------------------------------------------------\n' |
---|
| 1194 | |
---|
| 1195 | return str |
---|
| 1196 | |
---|
| 1197 | |
---|
| 1198 | ## |
---|
| 1199 | # @brief ?? |
---|
| 1200 | # @param sww_file ?? |
---|
| 1201 | # @param time ?? |
---|
| 1202 | # @param interpolation_points ?? |
---|
| 1203 | # @param quantity_names ?? |
---|
| 1204 | # @param verbose ?? |
---|
| 1205 | # @note Obsolete. Use file_function() in utils. |
---|
| 1206 | def interpolate_sww(sww_file, time, interpolation_points, |
---|
| 1207 | quantity_names=None, verbose=False): |
---|
| 1208 | """ |
---|
| 1209 | obsolete. |
---|
| 1210 | use file_function in utils |
---|
| 1211 | """ |
---|
| 1212 | |
---|
| 1213 | #open sww file |
---|
| 1214 | x, y, volumes, time, quantities = read_sww(sww_file) |
---|
[7317] | 1215 | log.critical("x=%s" % str(x)) |
---|
| 1216 | log.critical("y=%s" % str(y)) |
---|
[7276] | 1217 | |
---|
[7317] | 1218 | log.critical("time=%s" % str(time)) |
---|
| 1219 | log.critical("quantities=%s" % str(quantities)) |
---|
[6073] | 1220 | |
---|
| 1221 | #Add the x and y together |
---|
[7276] | 1222 | vertex_coordinates = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]), |
---|
| 1223 | axis=1) |
---|
[6073] | 1224 | |
---|
| 1225 | #Will return the quantity values at the specified times and locations |
---|
| 1226 | interp = Interpolation_interface(time, |
---|
| 1227 | quantities, |
---|
[7276] | 1228 | quantity_names=quantity_names, |
---|
[6073] | 1229 | vertex_coordinates=vertex_coordinates, |
---|
| 1230 | triangles=volumes, |
---|
| 1231 | interpolation_points=interpolation_points, |
---|
| 1232 | verbose=verbose) |
---|
| 1233 | |
---|
| 1234 | |
---|
| 1235 | ## |
---|
| 1236 | # @brief ?? |
---|
| 1237 | # @param file_name Name of the .SWW file to read. |
---|
| 1238 | def read_sww(file_name): |
---|
| 1239 | """ |
---|
| 1240 | obsolete - Nothing should be calling this |
---|
[7276] | 1241 | |
---|
[6073] | 1242 | Read in an sww file. |
---|
[7276] | 1243 | |
---|
[6073] | 1244 | Input; |
---|
| 1245 | file_name - the sww file |
---|
[7276] | 1246 | |
---|
[6073] | 1247 | Output; |
---|
| 1248 | x - Vector of x values |
---|
| 1249 | y - Vector of y values |
---|
| 1250 | z - Vector of bed elevation |
---|
| 1251 | volumes - Array. Each row has 3 values, representing |
---|
| 1252 | the vertices that define the volume |
---|
| 1253 | time - Vector of the times where there is stage information |
---|
| 1254 | stage - array with respect to time and vertices (x,y) |
---|
| 1255 | """ |
---|
| 1256 | |
---|
| 1257 | msg = 'Function read_sww in interpolat.py is obsolete' |
---|
| 1258 | raise Exception, msg |
---|
[7276] | 1259 | |
---|
[6073] | 1260 | #FIXME Have this reader as part of data_manager? |
---|
[7276] | 1261 | |
---|
| 1262 | from Scientific.IO.NetCDF import NetCDFFile |
---|
[6073] | 1263 | import tempfile |
---|
| 1264 | import sys |
---|
| 1265 | import os |
---|
[7276] | 1266 | |
---|
[6073] | 1267 | #Check contents |
---|
| 1268 | #Get NetCDF |
---|
[7276] | 1269 | |
---|
[6073] | 1270 | # see if the file is there. Throw a QUIET IO error if it isn't |
---|
| 1271 | # I don't think I could implement the above |
---|
[7276] | 1272 | |
---|
[6073] | 1273 | #throws prints to screen if file not present |
---|
| 1274 | junk = tempfile.mktemp(".txt") |
---|
| 1275 | fd = open(junk,'w') |
---|
| 1276 | stdout = sys.stdout |
---|
| 1277 | sys.stdout = fd |
---|
[7276] | 1278 | fid = NetCDFFile(file_name, netcdf_mode_r) |
---|
[6073] | 1279 | sys.stdout = stdout |
---|
| 1280 | fd.close() |
---|
| 1281 | #clean up |
---|
[7276] | 1282 | os.remove(junk) |
---|
| 1283 | |
---|
[6073] | 1284 | # Get the variables |
---|
| 1285 | x = fid.variables['x'][:] |
---|
| 1286 | y = fid.variables['y'][:] |
---|
[7276] | 1287 | volumes = fid.variables['volumes'][:] |
---|
[6073] | 1288 | time = fid.variables['time'][:] |
---|
| 1289 | |
---|
| 1290 | keys = fid.variables.keys() |
---|
| 1291 | keys.remove("x") |
---|
| 1292 | keys.remove("y") |
---|
| 1293 | keys.remove("volumes") |
---|
| 1294 | keys.remove("time") |
---|
[7276] | 1295 | #Turn NetCDF objects into numeric arrays |
---|
[6073] | 1296 | quantities = {} |
---|
| 1297 | for name in keys: |
---|
| 1298 | quantities[name] = fid.variables[name][:] |
---|
[7276] | 1299 | |
---|
[6073] | 1300 | fid.close() |
---|
| 1301 | return x, y, volumes, time, quantities |
---|
| 1302 | |
---|
| 1303 | |
---|
| 1304 | #------------------------------------------------------------- |
---|
| 1305 | if __name__ == "__main__": |
---|
| 1306 | names = ["x","y"] |
---|
| 1307 | someiterable = [[1,2],[3,4]] |
---|
| 1308 | csvwriter = writer(file("some.csv", "wb")) |
---|
| 1309 | csvwriter.writerow(names) |
---|
| 1310 | for row in someiterable: |
---|
| 1311 | csvwriter.writerow(row) |
---|