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