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