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