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