Changeset 2201 for inundation/fit_interpolate
- Timestamp:
- Jan 11, 2006, 5:32:23 PM (19 years ago)
- Location:
- inundation/fit_interpolate
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/fit_interpolate/interpolate.py
r2192 r2201 15 15 * remove points outside the mesh ?(in interpolate_block)? 16 16 * geo-ref (in interpolate_block) 17 * add functional interpolate interface - in mesh and points, out interp data 17 18 """ 18 19 … … 37 38 vertex_coordinates, 38 39 triangles, 39 #point_coordinates = None,40 40 mesh_origin = None, 41 41 verbose=False, 42 max_ points_per_cell=30):42 max_vertices_per_cell=30): 43 43 44 44 … … 49 49 50 50 vertex_coordinates: List of coordinate pairs [xi, eta] of 51 points constituting mesh (or aan m x 2 Numeric array)52 Points may appear multiple times53 (e.g. if vertices have discontinuities)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 54 55 55 triangles: List of 3-tuples (or a Numeric array) of 56 integers representing indices of all vertices in the mesh. 57 58 point_coordinates: List of coordinate pairs [x, y] of 59 data points (or an nx2 Numeric array) 60 If point_coordinates is absent, only smoothing matrix will 61 be built 62 63 alpha: Smoothing parameter 64 65 data_origin and mesh_origin are 3-tuples consisting of 66 UTM zone, easting and northing. If specified 67 point coordinates and vertex coordinates are assumed to be 68 relative to their respective origins. 69 70 """ 71 72 from pyvolution.util import ensure_numeric 73 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 74 68 # Initialise variabels 75 self. A_can_be_reused = False76 self. point_coordinates = None69 self._A_can_be_reused = False 70 self._point_coordinates = None 77 71 78 72 #Convert input to Numeric arrays … … 102 96 103 97 self.root = build_quadtree(self.mesh, 104 max_points_per_cell = max_ points_per_cell)98 max_points_per_cell = max_vertices_per_cell) 105 99 106 100 … … 148 142 #For each data_coordinate point 149 143 if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n) 150 151 144 x = point_coordinates[i] 152 153 145 element_found, sigma0, sigma1, sigma2, k = \ 154 146 search_tree_of_vertices(self.root, self.mesh, x) … … 168 160 else: 169 161 print 'Could not find triangle for point', x 170 171 162 return A 172 163 … … 281 272 per underlying matrix-matrix multiplication 282 273 point_coordinates: Interpolate mesh data to these positions. 274 List of coordinate pairs [x, y] of 275 data points (or an nx2 Numeric array) 276 If point_coordinates is absent, the points inputted last time 277 this method was called are used, if possible. 283 278 start_blocking_len: If the # of points is more or greater than this, 284 279 start blocking … … 290 285 # Can I interpolate, based on previous point_coordinates? 291 286 if point_coordinates is None: 292 if self. A_can_be_reused is True and \293 len(self. point_coordinates) < start_blocking_len:287 if self._A_can_be_reused is True and \ 288 len(self._point_coordinates) < start_blocking_len: 294 289 z = self._get_point_data_z(f, 295 290 verbose=verbose) 296 elif self. point_coordinates is not None:291 elif self._point_coordinates is not None: 297 292 # if verbose, give warning 298 293 if verbose: 299 294 print 'WARNING: Recalculating A matrix, due to blocking.' 300 point_coordinates = self. point_coordinates295 point_coordinates = self._point_coordinates 301 296 else: 302 297 #There are no good point_coordinates. import sys; sys.exit() … … 305 300 306 301 if point_coordinates is not None: 307 self. point_coordinates = point_coordinates302 self._point_coordinates = point_coordinates 308 303 if len(point_coordinates) < start_blocking_len or \ 309 304 start_blocking_len == 0: 310 self. A_can_be_reused = True305 self._A_can_be_reused = True 311 306 z = self.interpolate_block(f, point_coordinates, 312 307 verbose=verbose) 313 308 else: 314 309 #Handle blocking 315 self. A_can_be_reused = False310 self._A_can_be_reused = False 316 311 start=0 317 312 z = self.interpolate_block(f, point_coordinates[0:0]) … … 337 332 """ 338 333 if point_coordinates is not None: 339 self. A =self._build_interpolation_matrix_A(point_coordinates,334 self._A =self._build_interpolation_matrix_A(point_coordinates, 340 335 verbose=verbose) 341 336 return self._get_point_data_z(f) 342 337 343 338 def _get_point_data_z(self, f, verbose=False): 344 return self. A * f339 return self._A * f 345 340 346 341 #------------------------------------------------------------- -
inundation/fit_interpolate/search_functions.py
r2190 r2201 7 7 """ 8 8 from Numeric import dot 9 10 11 12 import time13 14 from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \15 ArrayType, allclose, take16 17 from pyvolution.mesh import Mesh18 from pyvolution.sparse import Sparse, Sparse_CSR19 from pyvolution.cg_solve import conjugate_gradient, VectorShapeError20 from coordinate_transforms.geo_reference import Geo_reference21 from pyvolution.quad import build_quadtree22 from utilities.numerical_tools import ensure_numeric23 from utilities.polygon import inside_polygon24 9 25 10 -
inundation/fit_interpolate/test_interpolate.py
r2190 r2201 67 67 data = [ [4,4] ] 68 68 interp = Interpolate(points, triangles, 69 max_ points_per_cell = 4)69 max_vertices_per_cell = 4) 70 70 #print "PDSG - interp.get_A()", interp.get_A() 71 71 answer = [ [ 0.06666667, 0.46666667, 0.46666667, 0., … … 507 507 #print "answer",answer 508 508 assert allclose(z, answer) 509 assert allclose(interp. A_can_be_reused, True)509 assert allclose(interp._A_can_be_reused, True) 510 510 511 511 z = interp.interpolate(f) … … 515 515 z = interp.interpolate(f, start_blocking_len=10) 516 516 assert allclose(z, answer) 517 assert allclose(interp. A_can_be_reused, False)517 assert allclose(interp._A_can_be_reused, False) 518 518 519 519 #A is recalculated 520 520 z = interp.interpolate(f) 521 521 assert allclose(z, answer) 522 assert allclose(interp. A_can_be_reused, True)523 522 assert allclose(interp._A_can_be_reused, True) 523 524 524 interp = Interpolate(vertices, triangles) 525 525 #Must raise an exception, no points specified
Note: See TracChangeset
for help on using the changeset viewer.