Changeset 7276 for anuga_core/source/anuga/fit_interpolate/interpolate.py
- Timestamp:
- Jun 30, 2009, 2:07:41 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/fit_interpolate/interpolate.py
r6621 r7276 41 41 42 42 43 import Numericas num43 import numpy as num 44 44 45 45 … … 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 79 (or an m x 2 Numeric array or78 points constituting a mesh 79 (or an m x 2 numeric array or 80 80 a geospatial object) 81 81 Points may appear multiple times 82 82 (e.g. if vertices have discontinuities) 83 83 84 triangles: List of 3-tuples (or a Numeric array) of85 integers representing indices of all vertices 84 triangles: List of 3-tuples (or a numeric array) of 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 data points or an nx2 Numeric array or a94 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 points constituting a mesh (or an m x 2 Numeric array or183 points constituting a mesh (or an m x 2 numeric array or 184 184 a geospatial object) 185 185 Points may appear multiple times 186 186 (e.g. if vertices have discontinuities) 187 187 188 triangles: List of 3-tuples (or a Numeric array) of188 triangles: List of 3-tuples (or a numeric array) of 189 189 integers representing indices of all vertices in the mesh. 190 190 … … 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 data points or an nx2 Numeric array or a Geospatial_data object246 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 296 f = ensure_numeric(f, num. Float)295 296 f = ensure_numeric(f, num.float) 297 297 if len(f.shape) > 1: 298 z = num.zeros((0, f.shape[1]) )298 z = num.zeros((0, f.shape[1]), num.int) #array default# 299 299 else: 300 z = num.zeros((0,) )301 300 z = num.zeros((0,), num.int) #array default# 301 302 302 for end in range(start_blocking_len, 303 303 len(point_coordinates), … … 305 305 t = self.interpolate_block(f, point_coordinates[start:end], 306 306 verbose=verbose) 307 z = num.concatenate((z, t) )307 z = num.concatenate((z, t), axis=0) #??default# 308 308 start = end 309 309 … … 311 311 t = self.interpolate_block(f, point_coordinates[start:end], 312 312 verbose=verbose) 313 z = num.concatenate((z, t) )313 z = num.concatenate((z, t), axis=0) #??default# 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 … … 340 340 point_coordinates = point_coordinates.get_data_points(absolute=True) 341 341 342 # Convert lists to Numeric arrays if necessary343 point_coordinates = ensure_numeric(point_coordinates, num. Float)344 f = ensure_numeric(f, num. Float)342 # Convert lists to numeric arrays if necessary 343 point_coordinates = ensure_numeric(point_coordinates, 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 … … 447 447 if verbose: print 'Building interpolation matrix' 448 448 449 # Convert point_coordinates to Numeric arrays, in case it was a list.450 point_coordinates = ensure_numeric(point_coordinates, num. Float)449 # Convert point_coordinates to numeric arrays, in case it was a list. 450 point_coordinates = ensure_numeric(point_coordinates, num.float) 451 451 452 452 if verbose: print 'Getting indices inside mesh boundary' … … 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 … … 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 ?? … … 526 526 points: Interpolate mesh data to these positions. 527 527 List of coordinate pairs [x, y] of 528 data points or an nx2 Numeric array or a Geospatial_data object529 528 data points or an nx2 numeric array or a Geospatial_data object 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) … … 694 694 695 695 Mandatory input 696 time: px1 array of monotonously increasing times ( Float)697 quantities: Dictionary of arrays or 1 array ( Float)696 time: px1 array of monotonously increasing times (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 705 705 imposing a particular order on the output vector. 706 vertex_coordinates: mx2 array of coordinates ( Float)707 triangles: nx3 array of indices into vertex_coordinates ( Int)708 interpolation_points: Nx2 array of coordinates to be interpolated to 706 vertex_coordinates: mx2 array of coordinates (float) 707 triangles: nx3 array of indices into vertex_coordinates (int) 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 = {} 811 811 for name in quantity_names: 812 q = quantities[name][:].flat 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 msg = 'Interpolation points must be an N x 2 Numeric array ' \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 … … 842 842 # mesh boundary as defined by triangles and vertex_coordinates. 843 843 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh 844 from anuga.utilities.polygon import outside_polygon 844 from anuga.utilities.polygon import outside_polygon 845 845 846 846 # Create temporary mesh object from mesh info passed 847 # into this function. 847 # into this function. 848 848 mesh = Mesh(vertex_coordinates, triangles) 849 849 mesh_boundary_polygon = mesh.get_boundary_polygon() … … 871 871 # FIXME (Ole): Why only Windoze? 872 872 from anuga.utilities.polygon import plot_polygons 873 #out_interp_pts = num.take(interpolation_points,[indices])874 873 title = ('Interpolation points fall ' 875 874 'outside specified mesh') … … 889 888 print msg 890 889 #raise Exception(msg) 891 890 892 891 elif triangles is None and vertex_coordinates is not None: #jj 893 892 #Dealing with sts file … … 915 914 m = len(self.interpolation_points) 916 915 p = len(self.time) 917 918 919 self.precomputed_values[name] = num.zeros((p, m), num. Float)916 917 for name in quantity_names: 918 self.precomputed_values[name] = num.zeros((p, m), num.float) 920 919 921 920 if verbose is True: 922 921 print 'Build interpolator' 923 922 924 923 925 924 # Build interpolator 926 925 if triangles is not None and vertex_coordinates is not None: 927 if verbose: 926 if verbose: 928 927 msg = 'Building interpolation matrix from source mesh ' 929 928 msg += '(%d vertices, %d triangles)' \ … … 931 930 triangles.shape[0]) 932 931 print msg 933 932 934 933 # This one is no longer needed for STS files 935 934 interpol = Interpolate(vertex_coordinates, 936 935 triangles, 937 verbose=verbose) 938 936 verbose=verbose) 937 939 938 elif triangles is None and vertex_coordinates is not None: 940 939 if verbose: … … 947 946 print 'Interpolating (%d interpolation points, %d timesteps).' \ 948 947 %(self.interpolation_points.shape[0], self.time.shape[0]), 949 948 950 949 if time_thinning > 1: 951 950 print 'Timesteps were thinned by a factor of %d' \ … … 954 953 print 955 954 956 955 for i, t in enumerate(self.time): 957 956 # Interpolate quantities at this timestep 958 957 if verbose and i%((p+10)/10) == 0: 959 958 print ' time step %d of %d' %(i, p) 960 959 961 960 for name in quantity_names: 962 961 if len(quantities[name].shape) == 2: … … 967 966 if verbose and i%((p+10)/10) == 0: 968 967 print ' quantity %s, size=%d' %(name, len(Q)) 969 970 # Interpolate 968 969 # Interpolate 971 970 if triangles is not None and vertex_coordinates is not None: 972 971 result = interpol.interpolate(Q, … … 980 979 interpolation_points=\ 981 980 self.interpolation_points) 982 981 983 982 #assert len(result), len(interpolation_points) 984 983 self.precomputed_values[name][i, :] = result … … 1007 1006 def __call__(self, t, point_id=None, x=None, y=None): 1008 1007 """Evaluate f(t) or f(t, point_id) 1009 1010 1011 1012 1013 1014 1015 are None an exception is raised 1016 1008 1009 Inputs: 1010 t: time - Model time. Must lie within existing timesteps 1011 point_id: index of one of the preprocessed points. 1012 1013 If spatial info is present and all of point_id 1014 are None an exception is raised 1015 1017 1016 If no spatial info is present, point_id arguments are ignored 1018 1017 making f a function of time only. 1019 1018 1020 FIXME: f(t, x, y) x, y could overrided location, point_id ignored 1021 1022 1019 FIXME: f(t, x, y) x, y could overrided location, point_id ignored 1020 FIXME: point_id could also be a slice 1021 FIXME: What if x and y are vectors? 1023 1022 FIXME: What about f(x,y) without t? 1024 1023 """ 1025 1024 1026 1025 from math import pi, cos, sin, sqrt 1027 from anuga.abstract_2d_finite_volumes.util import mean 1026 from anuga.abstract_2d_finite_volumes.util import mean 1028 1027 1029 1028 if self.spatial is True: … … 1060 1059 1061 1060 # Compute interpolated values 1062 q = num.zeros(len(self.quantity_names), num. Float)1063 1061 q = num.zeros(len(self.quantity_names), num.float) 1062 for i, name in enumerate(self.quantity_names): 1064 1063 Q = self.precomputed_values[name] 1065 1064 1066 1065 if self.spatial is False: 1067 # If there is no spatial info 1066 # If there is no spatial info 1068 1067 assert len(Q.shape) == 1 1069 1068 … … 1080 1079 Q1 = Q[self.index+1, point_id] 1081 1080 1082 # Linear temporal interpolation 1081 # Linear temporal interpolation 1083 1082 if ratio > 0: 1084 1083 if Q0 == NAN and Q1 == NAN: … … 1096 1095 # Replicate q according to x and y 1097 1096 # This is e.g used for Wind_stress 1098 if x is None or y is None: 1097 if x is None or y is None: 1099 1098 return q 1100 1099 else: … … 1109 1108 res = [] 1110 1109 for col in q: 1111 res.append(col*num.ones(N, num. Float))1112 1110 res.append(col*num.ones(N, num.float)) 1111 1113 1112 return res 1114 1113 … … 1126 1125 """Output statistics about interpolation_function 1127 1126 """ 1128 1127 1129 1128 vertex_coordinates = self.vertex_coordinates 1130 interpolation_points = self.interpolation_points 1129 interpolation_points = self.interpolation_points 1131 1130 quantity_names = self.quantity_names 1132 1131 #quantities = self.quantities 1133 precomputed_values = self.precomputed_values 1134 1132 precomputed_values = self.precomputed_values 1133 1135 1134 x = vertex_coordinates[:,0] 1136 y = vertex_coordinates[:,1] 1135 y = vertex_coordinates[:,1] 1137 1136 1138 1137 str = '------------------------------------------------\n' … … 1148 1147 for name in quantity_names: 1149 1148 minq, maxq = self.quantities_range[name] 1150 str += ' %s in [%f, %f]\n' %(name, minq, maxq) 1151 #q = quantities[name][:].flat 1149 str += ' %s in [%f, %f]\n' %(name, minq, maxq) 1150 #q = quantities[name][:].flatten() 1152 1151 #str += ' %s in [%f, %f]\n' %(name, min(q), max(q)) 1153 1152 1154 if interpolation_points is not None: 1153 if interpolation_points is not None: 1155 1154 str += ' Interpolation points (xi, eta):'\ 1156 1155 ' number of points == %d\n' %interpolation_points.shape[0] … … 1160 1159 max(interpolation_points[:,1])) 1161 1160 str += ' Interpolated quantities (over all timesteps):\n' 1162 1161 1163 1162 for name in quantity_names: 1164 q = precomputed_values[name][:].flat 1163 q = precomputed_values[name][:].flatten() 1165 1164 str += ' %s at interpolation points in [%f, %f]\n'\ 1166 1165 %(name, min(q), max(q)) … … 1189 1188 print "x",x 1190 1189 print "y",y 1191 1190 1192 1191 print "time", time 1193 1192 print "quantities", quantities 1194 1193 1195 1194 #Add the x and y together 1196 vertex_coordinates = num.concatenate((x[:,num.NewAxis], y[:,num.NewAxis]),axis=1) 1195 vertex_coordinates = num.concatenate((x[:,num.newaxis], y[:,num.newaxis]), 1196 axis=1) 1197 1197 1198 1198 #Will return the quantity values at the specified times and locations 1199 1199 interp = Interpolation_interface(time, 1200 1200 quantities, 1201 quantity_names=quantity_names, 1201 quantity_names=quantity_names, 1202 1202 vertex_coordinates=vertex_coordinates, 1203 1203 triangles=volumes, … … 1212 1212 """ 1213 1213 obsolete - Nothing should be calling this 1214 1214 1215 1215 Read in an sww file. 1216 1216 1217 1217 Input; 1218 1218 file_name - the sww file 1219 1219 1220 1220 Output; 1221 1221 x - Vector of x values … … 1230 1230 msg = 'Function read_sww in interpolat.py is obsolete' 1231 1231 raise Exception, msg 1232 1232 1233 1233 #FIXME Have this reader as part of data_manager? 1234 1235 from Scientific.IO.NetCDF import NetCDFFile 1234 1235 from Scientific.IO.NetCDF import NetCDFFile 1236 1236 import tempfile 1237 1237 import sys 1238 1238 import os 1239 1239 1240 1240 #Check contents 1241 1241 #Get NetCDF 1242 1242 1243 1243 # see if the file is there. Throw a QUIET IO error if it isn't 1244 1244 # I don't think I could implement the above 1245 1245 1246 1246 #throws prints to screen if file not present 1247 1247 junk = tempfile.mktemp(".txt") … … 1249 1249 stdout = sys.stdout 1250 1250 sys.stdout = fd 1251 fid = NetCDFFile(file_name, netcdf_mode_r) 1251 fid = NetCDFFile(file_name, netcdf_mode_r) 1252 1252 sys.stdout = stdout 1253 1253 fd.close() 1254 1254 #clean up 1255 os.remove(junk) 1256 1255 os.remove(junk) 1256 1257 1257 # Get the variables 1258 1258 x = fid.variables['x'][:] 1259 1259 y = fid.variables['y'][:] 1260 volumes = fid.variables['volumes'][:] 1260 volumes = fid.variables['volumes'][:] 1261 1261 time = fid.variables['time'][:] 1262 1262 … … 1266 1266 keys.remove("volumes") 1267 1267 keys.remove("time") 1268 #Turn NetCDF objects into Numeric arrays1268 #Turn NetCDF objects into numeric arrays 1269 1269 quantities = {} 1270 1270 for name in keys: 1271 1271 quantities[name] = fid.variables[name][:] 1272 1272 1273 1273 fid.close() 1274 1274 return x, y, volumes, time, quantities
Note: See TracChangeset
for help on using the changeset viewer.
