Changeset 6410 for branches/numpy/anuga/fit_interpolate/interpolate.py
- Timestamp:
- Feb 25, 2009, 9:37:22 AM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/fit_interpolate/interpolate.py
r6304 r6410 68 68 max_vertices_per_cell=None, 69 69 start_blocking_len=500000, 70 use_cache=False, 71 verbose=False): 70 use_cache=False, 71 verbose=False): 72 72 """Interpolate vertex_values to interpolation points. 73 73 74 74 Inputs (mandatory): 75 75 76 76 77 77 vertex_coordinates: List of coordinate pairs [xi, eta] of 78 points constituting a mesh 78 points constituting a mesh 79 79 (or an m x 2 numeric array or 80 80 a geospatial object) … … 83 83 84 84 triangles: List of 3-tuples (or a numeric array) of 85 integers representing indices of all vertices 85 integers representing indices of all vertices 86 86 in the mesh. 87 87 88 88 vertex_values: Vector or array of data at the mesh vertices. 89 89 If array, interpolation will be done for each column as 90 90 per underlying matrix-matrix multiplication 91 91 92 92 interpolation_points: Interpolate mesh data to these positions. 93 93 List of coordinate pairs [x, y] of 94 94 data points or an nx2 numeric array or a 95 95 Geospatial_data object 96 96 97 97 Inputs (optional) 98 98 99 99 mesh_origin: A geo_reference object or 3-tuples consisting of 100 100 UTM zone, easting and northing. … … 108 108 object and a mesh origin, since geospatial has its 109 109 own mesh origin. 110 110 111 111 start_blocking_len: If the # of points is more or greater than this, 112 start blocking 113 112 start blocking 113 114 114 use_cache: True or False 115 115 116 116 117 117 Output: 118 118 119 119 Interpolated values at specified point_coordinates 120 120 121 Note: This function is a simple shortcut for case where 121 Note: This function is a simple shortcut for case where 122 122 interpolation matrix is unnecessary 123 Note: This function does not take blocking into account, 123 Note: This function does not take blocking into account, 124 124 but allows caching. 125 125 126 126 """ 127 127 128 128 # FIXME(Ole): Probably obsolete since I is precomputed and 129 129 # interpolate_block caches 130 130 131 131 from anuga.caching import cache 132 132 133 133 # Create interpolation object with matrix 134 args = (ensure_numeric(vertex_coordinates, num.float), 134 args = (ensure_numeric(vertex_coordinates, num.float), 135 135 ensure_numeric(triangles)) 136 136 kwargs = {'mesh_origin': mesh_origin, 137 137 'max_vertices_per_cell': max_vertices_per_cell, 138 138 'verbose': verbose} 139 139 140 140 if use_cache is True: 141 141 if sys.platform != 'win32': … … 149 149 else: 150 150 I = apply(Interpolate, args, kwargs) 151 151 152 152 # Call interpolate method with interpolation points 153 153 result = I.interpolate_block(vertex_values, interpolation_points, 154 154 use_cache=use_cache, 155 155 verbose=verbose) 156 156 157 157 return result 158 158 159 159 160 160 ## 161 # @brief 161 # @brief 162 162 class Interpolate (FitInterpolate): 163 163 … … 181 181 Inputs: 182 182 vertex_coordinates: List of coordinate pairs [xi, eta] of 183 183 points constituting a mesh (or an m x 2 numeric array or 184 184 a geospatial object) 185 185 Points may appear multiple times … … 202 202 203 203 # FIXME (Ole): Need an input check 204 204 205 205 FitInterpolate.__init__(self, 206 206 vertex_coordinates=vertex_coordinates, … … 209 209 verbose=verbose, 210 210 max_vertices_per_cell=max_vertices_per_cell) 211 211 212 212 # Initialise variables 213 213 self._A_can_be_reused = False # FIXME (Ole): Probably obsolete … … 215 215 self.interpolation_matrices = {} # Store precomputed matrices 216 216 217 217 218 218 219 219 ## … … 243 243 point_coordinates: Interpolate mesh data to these positions. 244 244 List of coordinate pairs [x, y] of 245 246 247 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 248 this method was called are used, if possible. 249 249 250 250 start_blocking_len: If the # of points is more or greater than this, 251 start blocking 252 253 254 251 start blocking 252 253 Output: 254 Interpolated values at inputted points (z). 255 255 """ 256 256 … … 262 262 # This has now been addressed through an attempt in interpolate_block 263 263 264 if verbose: print 'Build intepolation object' 264 if verbose: print 'Build intepolation object' 265 265 if isinstance(point_coordinates, Geospatial_data): 266 266 point_coordinates = point_coordinates.get_data_points(absolute=True) … … 280 280 msg = 'ERROR (interpolate.py): No point_coordinates inputted' 281 281 raise Exception(msg) 282 282 283 283 if point_coordinates is not None: 284 284 self._point_coordinates = point_coordinates … … 293 293 start = 0 294 294 # creating a dummy array to concatenate to. 295 295 296 296 f = ensure_numeric(f, num.float) 297 297 if len(f.shape) > 1: … … 299 299 else: 300 300 z = num.zeros((0,), num.int) #array default# 301 301 302 302 for end in range(start_blocking_len, 303 303 len(point_coordinates), … … 313 313 z = num.concatenate((z, t)) 314 314 return z 315 315 316 316 317 317 ## … … 322 322 # @param verbose True if this function is verbose. 323 323 # @return ?? 324 def interpolate_block(self, f, point_coordinates, 324 def interpolate_block(self, f, point_coordinates, 325 325 use_cache=False, verbose=False): 326 326 """ … … 329 329 330 330 Return the point data, z. 331 331 332 332 See interpolate for doc info. 333 333 """ 334 334 335 335 # FIXME (Ole): I reckon we should change the interface so that 336 # the user can specify the interpolation matrix instead of the 336 # the user can specify the interpolation matrix instead of the 337 337 # interpolation points to save time. 338 338 … … 342 342 # Convert lists to numeric arrays if necessary 343 343 point_coordinates = ensure_numeric(point_coordinates, num.float) 344 f = ensure_numeric(f, num.float) 344 f = ensure_numeric(f, num.float) 345 345 346 346 from anuga.caching import myhash 347 import sys 348 347 import sys 348 349 349 if use_cache is True: 350 350 if sys.platform != 'win32': 351 351 # FIXME (Ole): (Why doesn't this work on windoze?) 352 352 # Still absolutely fails on Win 24 Oct 2008 353 353 354 354 X = cache(self._build_interpolation_matrix_A, 355 355 args=(point_coordinates,), 356 kwargs={'verbose': verbose}, 357 verbose=verbose) 356 kwargs={'verbose': verbose}, 357 verbose=verbose) 358 358 else: 359 359 # FIXME … … 361 361 # This will work on Linux as well if we want to use it there. 362 362 key = myhash(point_coordinates) 363 364 reuse_A = False 365 363 364 reuse_A = False 365 366 366 if self.interpolation_matrices.has_key(key): 367 X, stored_points = self.interpolation_matrices[key] 367 X, stored_points = self.interpolation_matrices[key] 368 368 if num.alltrue(stored_points == point_coordinates): 369 reuse_A = True 370 369 reuse_A = True # Reuse interpolation matrix 370 371 371 if reuse_A is False: 372 372 X = self._build_interpolation_matrix_A(point_coordinates, … … 375 375 else: 376 376 X = self._build_interpolation_matrix_A(point_coordinates, 377 verbose=verbose) 378 379 # Unpack result 377 verbose=verbose) 378 379 # Unpack result 380 380 self._A, self.inside_poly_indices, self.outside_poly_indices = X 381 381 … … 389 389 msg += ' I got %d points and %d matrix rows.' \ 390 390 % (point_coordinates.shape[0], self._A.shape[0]) 391 assert point_coordinates.shape[0] == self._A.shape[0], msg 391 assert point_coordinates.shape[0] == self._A.shape[0], msg 392 392 393 393 msg = 'The number of columns in matrix A must be the same as the ' 394 394 msg += 'number of mesh vertices.' 395 395 msg += ' I got %d vertices and %d matrix columns.' \ 396 % (f.shape[0], self._A.shape[1]) 396 % (f.shape[0], self._A.shape[1]) 397 397 assert self._A.shape[1] == f.shape[0], msg 398 398 … … 409 409 """ 410 410 Return the point data, z. 411 411 412 412 Precondition: The _A matrix has been created 413 413 """ … … 416 416 417 417 # Taking into account points outside the mesh. 418 for i in self.outside_poly_indices: 418 for i in self.outside_poly_indices: 419 419 z[i] = NAN 420 420 return z … … 456 456 self.mesh.get_boundary_polygon(), 457 457 closed=True, verbose=verbose) 458 458 459 459 # Build n x m interpolation matrix 460 460 if verbose and len(outside_poly_indices) > 0: 461 461 print '\n WARNING: Points outside mesh boundary. \n' 462 462 463 463 # Since you can block, throw a warning, not an error. 464 464 if verbose and 0 == len(inside_poly_indices): 465 465 print '\n WARNING: No points within the mesh! \n' 466 466 467 467 m = self.mesh.number_of_nodes # Nbr of basis functions (1/vertex) 468 468 n = point_coordinates.shape[0] # Nbr of data points … … 474 474 475 475 n = len(inside_poly_indices) 476 476 477 477 # Compute matrix elements for points inside the mesh 478 478 if verbose: print 'Building interpolation matrix from %d points' %n … … 485 485 element_found, sigma0, sigma1, sigma2, k = \ 486 486 search_tree_of_vertices(self.root, self.mesh, x) 487 488 487 488 # Update interpolation matrix A if necessary 489 489 if element_found is True: 490 490 # Assign values to matrix A … … 499 499 A[i, j] = sigmas[j] 500 500 else: 501 msg = 'Could not find triangle for point', x 501 msg = 'Could not find triangle for point', x 502 502 raise Exception(msg) 503 503 return A, inside_poly_indices, outside_poly_indices 504 504 505 505 506 507 508 509 506 507 508 509 510 510 ## 511 511 # @brief ?? … … 527 527 List of coordinate pairs [x, y] of 528 528 data points or an nx2 numeric array or a Geospatial_data object 529 529 530 530 No test for this yet. 531 531 Note, this has no time the input data has no time dimension. Which is 532 532 different from most of the data we interpolate, eg sww info. 533 533 534 534 Output: 535 535 Interpolated values at inputted points. … … 537 537 538 538 interp = Interpolate(vertices, 539 triangles, 539 triangles, 540 540 max_vertices_per_cell=max_points_per_cell, 541 541 mesh_origin=mesh_origin) 542 542 543 543 calc = interp.interpolate(vertex_attributes, 544 544 points, … … 554 554 # @param velocity_y_file Name of the output y velocity file. 555 555 # @param stage_file Name of the output stage file. 556 # @param froude_file 556 # @param froude_file 557 557 # @param time_thinning Time thinning step to use. 558 558 # @param verbose True if this function is to be verbose. … … 604 604 time_thinning=time_thinning, 605 605 use_cache=use_cache) 606 606 607 607 depth_writer = writer(file(depth_file, "wb")) 608 608 velocity_x_writer = writer(file(velocity_x_file, "wb")) … … 620 620 velocity_y_writer.writerow(heading) 621 621 if stage_file is not None: 622 stage_writer.writerow(heading) 622 stage_writer.writerow(heading) 623 623 if froude_file is not None: 624 froude_writer.writerow(heading) 625 624 froude_writer.writerow(heading) 625 626 626 for time in callable_sww.get_time(): 627 627 depths = [time] 628 628 velocity_xs = [time] 629 629 velocity_ys = [time] 630 if stage_file is not None: 631 stages = [time] 632 if froude_file is not None: 633 froudes = [time] 630 if stage_file is not None: 631 stages = [time] 632 if froude_file is not None: 633 froudes = [time] 634 634 for point_i, point in enumerate(points): 635 635 quantities = callable_sww(time,point_i) 636 636 637 637 w = quantities[0] 638 638 z = quantities[1] 639 639 momentum_x = quantities[2] 640 640 momentum_y = quantities[3] 641 depth = w - z 642 641 depth = w - z 642 643 643 if w == NAN or z == NAN or momentum_x == NAN: 644 644 velocity_x = NAN … … 677 677 678 678 if stage_file is not None: 679 stage_writer.writerow(stages) 679 stage_writer.writerow(stages) 680 680 if froude_file is not None: 681 froude_writer.writerow(froudes) 681 froude_writer.writerow(froudes) 682 682 683 683 684 684 ## 685 # @brief 685 # @brief 686 686 class Interpolation_function: 687 687 """Interpolation_interface - creates callable object f(t, id) or f(t,x,y) … … 695 695 Mandatory input 696 696 time: px1 array of monotonously increasing times (float) 697 quantities: Dictionary of arrays or 1 array (float) 697 quantities: Dictionary of arrays or 1 array (float) 698 698 The arrays must either have dimensions pxm or mx1. 699 699 The resulting function will be time dependent in 700 700 the former case while it will be constant with 701 701 respect to time in the latter case. 702 702 703 703 Optional input: 704 704 quantity_names: List of keys into the quantities dictionary for … … 706 706 vertex_coordinates: mx2 array of coordinates (float) 707 707 triangles: nx3 array of indices into vertex_coordinates (Int) 708 interpolation_points: Nx2 array of coordinates to be interpolated to 708 interpolation_points: Nx2 array of coordinates to be interpolated to 709 709 verbose: Level of reporting 710 710 … … 742 742 time, 743 743 quantities, 744 quantity_names=None, 744 quantity_names=None, 745 745 vertex_coordinates=None, 746 746 triangles=None, … … 762 762 print 'Interpolation_function: input checks' 763 763 764 764 # Check temporal info 765 765 time = ensure_numeric(time) 766 766 if not num.alltrue(time[1:] - time[:-1] >= 0): … … 791 791 if triangles is not None: 792 792 triangles = ensure_numeric(triangles) 793 self.spatial = True 793 self.spatial = True 794 794 795 795 if verbose is True: 796 796 print 'Interpolation_function: thinning by %d' % time_thinning 797 797 798 798 799 799 # Thin timesteps if needed 800 800 # Note array() is used to make the thinned arrays contiguous in memory 801 self.time = num.array(time[::time_thinning]) 801 self.time = num.array(time[::time_thinning]) 802 802 for name in quantity_names: 803 803 if len(quantities[name].shape) == 2: … … 806 806 if verbose is True: 807 807 print 'Interpolation_function: precomputing' 808 808 809 809 # Save for use with statistics 810 810 self.quantities_range = {} … … 812 812 q = quantities[name][:].flatten() 813 813 self.quantities_range[name] = [min(q), max(q)] 814 815 self.quantity_names = quantity_names 816 self.vertex_coordinates = vertex_coordinates 814 815 self.quantity_names = quantity_names 816 self.vertex_coordinates = vertex_coordinates 817 817 self.interpolation_points = interpolation_points 818 818 … … 825 825 #if self.spatial is False: 826 826 # raise 'Triangles and vertex_coordinates must be specified' 827 # 827 # 828 828 try: 829 829 self.interpolation_points = \ 830 830 interpolation_points = ensure_numeric(interpolation_points) 831 831 except: 832 832 msg = 'Interpolation points must be an N x 2 numeric array ' \ 833 833 'or a list of points\n' 834 834 msg += 'Got: %s.' %(str(self.interpolation_points)[:60] + '...') 835 835 raise msg 836 836 … … 839 839 # mesh boundary as defined by triangles and vertex_coordinates. 840 840 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh 841 from anuga.utilities.polygon import outside_polygon 841 from anuga.utilities.polygon import outside_polygon 842 842 843 843 # Create temporary mesh object from mesh info passed 844 # into this function. 844 # into this function. 845 845 mesh = Mesh(vertex_coordinates, triangles) 846 846 mesh_boundary_polygon = mesh.get_boundary_polygon() … … 868 868 # FIXME (Ole): Why only Windoze? 869 869 from anuga.utilities.polygon import plot_polygons 870 #out_interp_pts = num.take(interpolation_points,[indices])871 870 title = ('Interpolation points fall ' 872 871 'outside specified mesh') … … 886 885 print msg 887 886 #raise Exception(msg) 888 887 889 888 elif triangles is None and vertex_coordinates is not None: #jj 890 889 #Dealing with sts file … … 911 910 m = len(self.interpolation_points) 912 911 p = len(self.time) 913 914 912 913 for name in quantity_names: 915 914 self.precomputed_values[name] = num.zeros((p, m), num.float) 916 915 … … 918 917 print 'Build interpolator' 919 918 920 919 921 920 # Build interpolator 922 921 if triangles is not None and vertex_coordinates is not None: 923 if verbose: 922 if verbose: 924 923 msg = 'Building interpolation matrix from source mesh ' 925 924 msg += '(%d vertices, %d triangles)' \ … … 927 926 triangles.shape[0]) 928 927 print msg 929 928 930 929 # This one is no longer needed for STS files 931 930 interpol = Interpolate(vertex_coordinates, 932 931 triangles, 933 verbose=verbose) 934 932 verbose=verbose) 933 935 934 elif triangles is None and vertex_coordinates is not None: 936 935 if verbose: … … 943 942 print 'Interpolating (%d interpolation points, %d timesteps).' \ 944 943 %(self.interpolation_points.shape[0], self.time.shape[0]), 945 944 946 945 if time_thinning > 1: 947 946 print 'Timesteps were thinned by a factor of %d' \ … … 950 949 print 951 950 952 951 for i, t in enumerate(self.time): 953 952 # Interpolate quantities at this timestep 954 953 if verbose and i%((p+10)/10) == 0: 955 954 print ' time step %d of %d' %(i, p) 956 955 957 956 for name in quantity_names: 958 957 if len(quantities[name].shape) == 2: … … 963 962 if verbose and i%((p+10)/10) == 0: 964 963 print ' quantity %s, size=%d' %(name, len(Q)) 965 966 # Interpolate 964 965 # Interpolate 967 966 if triangles is not None and vertex_coordinates is not None: 968 967 result = interpol.interpolate(Q, … … 976 975 interpolation_points=\ 977 976 self.interpolation_points) 978 977 979 978 #assert len(result), len(interpolation_points) 980 979 self.precomputed_values[name][i, :] = result … … 1003 1002 def __call__(self, t, point_id=None, x=None, y=None): 1004 1003 """Evaluate f(t) or f(t, point_id) 1005 1006 1007 1008 1009 1010 1011 are None an exception is raised 1012 1004 1005 Inputs: 1006 t: time - Model time. Must lie within existing timesteps 1007 point_id: index of one of the preprocessed points. 1008 1009 If spatial info is present and all of point_id 1010 are None an exception is raised 1011 1013 1012 If no spatial info is present, point_id arguments are ignored 1014 1013 making f a function of time only. 1015 1014 1016 FIXME: f(t, x, y) x, y could overrided location, point_id ignored 1017 1018 1015 FIXME: f(t, x, y) x, y could overrided location, point_id ignored 1016 FIXME: point_id could also be a slice 1017 FIXME: What if x and y are vectors? 1019 1018 FIXME: What about f(x,y) without t? 1020 1019 """ 1021 1020 1022 1021 from math import pi, cos, sin, sqrt 1023 from anuga.abstract_2d_finite_volumes.util import mean 1022 from anuga.abstract_2d_finite_volumes.util import mean 1024 1023 1025 1024 if self.spatial is True: … … 1057 1056 # Compute interpolated values 1058 1057 q = num.zeros(len(self.quantity_names), num.float) 1059 1058 for i, name in enumerate(self.quantity_names): 1060 1059 Q = self.precomputed_values[name] 1061 1060 1062 1061 if self.spatial is False: 1063 # If there is no spatial info 1062 # If there is no spatial info 1064 1063 assert len(Q.shape) == 1 1065 1064 … … 1076 1075 Q1 = Q[self.index+1, point_id] 1077 1076 1078 # Linear temporal interpolation 1077 # Linear temporal interpolation 1079 1078 if ratio > 0: 1080 1079 if Q0 == NAN and Q1 == NAN: … … 1092 1091 # Replicate q according to x and y 1093 1092 # This is e.g used for Wind_stress 1094 if x is None or y is None: 1093 if x is None or y is None: 1095 1094 return q 1096 1095 else: … … 1106 1105 for col in q: 1107 1106 res.append(col*num.ones(N, num.float)) 1108 1107 1109 1108 return res 1110 1109 … … 1122 1121 """Output statistics about interpolation_function 1123 1122 """ 1124 1123 1125 1124 vertex_coordinates = self.vertex_coordinates 1126 interpolation_points = self.interpolation_points 1125 interpolation_points = self.interpolation_points 1127 1126 quantity_names = self.quantity_names 1128 1127 #quantities = self.quantities 1129 precomputed_values = self.precomputed_values 1130 1128 precomputed_values = self.precomputed_values 1129 1131 1130 x = vertex_coordinates[:,0] 1132 y = vertex_coordinates[:,1] 1131 y = vertex_coordinates[:,1] 1133 1132 1134 1133 str = '------------------------------------------------\n' … … 1144 1143 for name in quantity_names: 1145 1144 minq, maxq = self.quantities_range[name] 1146 str += ' %s in [%f, %f]\n' %(name, minq, maxq) 1145 str += ' %s in [%f, %f]\n' %(name, minq, maxq) 1147 1146 #q = quantities[name][:].flatten() 1148 1147 #str += ' %s in [%f, %f]\n' %(name, min(q), max(q)) 1149 1148 1150 if interpolation_points is not None: 1149 if interpolation_points is not None: 1151 1150 str += ' Interpolation points (xi, eta):'\ 1152 1151 ' number of points == %d\n' %interpolation_points.shape[0] … … 1156 1155 max(interpolation_points[:,1])) 1157 1156 str += ' Interpolated quantities (over all timesteps):\n' 1158 1157 1159 1158 for name in quantity_names: 1160 1159 q = precomputed_values[name][:].flatten() … … 1185 1184 print "x",x 1186 1185 print "y",y 1187 1186 1188 1187 print "time", time 1189 1188 print "quantities", quantities … … 1195 1194 interp = Interpolation_interface(time, 1196 1195 quantities, 1197 quantity_names=quantity_names, 1196 quantity_names=quantity_names, 1198 1197 vertex_coordinates=vertex_coordinates, 1199 1198 triangles=volumes, … … 1208 1207 """ 1209 1208 obsolete - Nothing should be calling this 1210 1209 1211 1210 Read in an sww file. 1212 1211 1213 1212 Input; 1214 1213 file_name - the sww file 1215 1214 1216 1215 Output; 1217 1216 x - Vector of x values … … 1226 1225 msg = 'Function read_sww in interpolat.py is obsolete' 1227 1226 raise Exception, msg 1228 1227 1229 1228 #FIXME Have this reader as part of data_manager? 1230 1231 from Scientific.IO.NetCDF import NetCDFFile 1229 1230 from Scientific.IO.NetCDF import NetCDFFile 1232 1231 import tempfile 1233 1232 import sys 1234 1233 import os 1235 1234 1236 1235 #Check contents 1237 1236 #Get NetCDF 1238 1237 1239 1238 # see if the file is there. Throw a QUIET IO error if it isn't 1240 1239 # I don't think I could implement the above 1241 1240 1242 1241 #throws prints to screen if file not present 1243 1242 junk = tempfile.mktemp(".txt") … … 1245 1244 stdout = sys.stdout 1246 1245 sys.stdout = fd 1247 fid = NetCDFFile(file_name, netcdf_mode_r) 1246 fid = NetCDFFile(file_name, netcdf_mode_r) 1248 1247 sys.stdout = stdout 1249 1248 fd.close() 1250 1249 #clean up 1251 os.remove(junk) 1252 1250 os.remove(junk) 1251 1253 1252 # Get the variables 1254 1253 x = fid.variables['x'][:] 1255 1254 y = fid.variables['y'][:] 1256 volumes = fid.variables['volumes'][:] 1255 volumes = fid.variables['volumes'][:] 1257 1256 time = fid.variables['time'][:] 1258 1257 … … 1266 1265 for name in keys: 1267 1266 quantities[name] = fid.variables[name][:] 1268 1267 1269 1268 fid.close() 1270 1269 return x, y, volumes, time, quantities
Note: See TracChangeset
for help on using the changeset viewer.
