Changeset 5361
 Timestamp:
 May 24, 2008, 11:55:52 AM (16 years ago)
 Location:
 anuga_core/source/anuga
 Files:

 4 edited
Legend:
 Unmodified
 Added
 Removed

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py
r5349 r5361 19 19 from anuga.geospatial_data.geospatial_data import ensure_absolute 20 20 from math import sqrt, atan, degrees 21 from exceptions import IOError 21 22 22 23 23 … … 110 110 # Instead we pass in those attributes that are needed (and return them 111 111 # if modified) 112 113 112 kwargs = {'quantities': quantities, 114 113 'interpolation_points': interpolation_points, … … 179 178 try: 180 179 fid = open(filename) 181 except IOError, e:180 except Exception, e: 182 181 msg = 'File "%s" could not be opened: Error="%s"'\ 183 182 %(filename, e) 184 raise IOError, msg # So IOErrors can be caught183 raise msg 185 184 186 185 line = fid.readline() … … 244 243 raise Exception, msg 245 244 246 245 247 246 if interpolation_points is not None: 248 247 interpolation_points = ensure_absolute(interpolation_points) … … 250 249 assert interpolation_points.shape[1] == 2, msg 251 250 252 if verbose:253 print 'File_function using quantities %s from file %s' %(str(quantity_names), filename)254 251 255 252 # Now assert that requested quantitites (and the independent ones) … … 273 270 274 271 if filename[3:] == 'tms' and spatial is True: 275 msg = 'Files of type tms must not contain spatial information'272 msg = 'Files of type tms must not contain spatial information' 276 273 raise msg 277 274 278 275 if filename[3:] == 'sww' and spatial is False: 279 276 msg = 'Files of type sww must contain spatial information' 280 raise msg 277 raise msg 278 279 if filename[3:] == 'sts' and spatial is False: 280 #What if mux file only contains one point 281 msg = 'Files of type sts must contain spatial information' 282 raise msg 281 283 282 284 # Get first timestep … … 300 302 x = fid.variables['x'][:] 301 303 y = fid.variables['y'][:] 302 triangles = fid.variables['volumes'][:] 304 if filename[3:] == 'sww': 305 triangles = fid.variables['volumes'][:] 303 306 304 307 x = reshape(x, (len(x),1)) … … 308 311 if interpolation_points is not None: 309 312 # Adjust for georef 310 311 # FIXME (Ole): Use geo_reference.get_relative312 313 interpolation_points[:,0] = xllcorner 313 314 interpolation_points[:,1] = yllcorner 314 315 316 317 315 318 316 if domain_starttime is not None: … … 350 348 351 349 if not spatial: 352 vertex_coordinates = triangles = interpolation_points = None 353 350 vertex_coordinates = triangles = interpolation_points = None 351 if filename[3:] == 'sts':#added 352 triangles = None 353 #vertex coordinates is position of urs gauges 354 354 355 355 # Return Interpolation_function instance as well as … … 782 782 #print 'label', label 783 783 leg_label.append(label) 784 785 784 786 785 … … 2145 2144 2146 2145 callable_sww = file_function(sww_file, 2147 2148 2149 2150 2146 quantities=core_quantities, 2147 interpolation_points=points_array, 2148 verbose=verbose, 2149 use_cache=use_cache) 2151 2150 2152 2151 gauge_file='gauge_' 
anuga_core/source/anuga/fit_interpolate/general_fit_interpolate.py
r5321 r5361 92 92 max_vertices_per_cell = MAX_VERTICES_PER_CELL 93 93 if mesh is None: 94 # Fixme (DSG) Throw errors if triangles or vertex_coordinates 95 # are None 94 if vertex_coordinates is not None and triangles is not None: 95 # Fixme (DSG) Throw errors if triangles or vertex_coordinates 96 # are None 96 97 97 #Convert input to Numeric arrays98 triangles = ensure_numeric(triangles, Int)99 vertex_coordinates = ensure_absolute(vertex_coordinates,98 #Convert input to Numeric arrays 99 triangles = ensure_numeric(triangles, Int) 100 vertex_coordinates = ensure_absolute(vertex_coordinates, 100 101 geo_reference = mesh_origin) 101 102 102 if verbose: print 'FitInterpolate: Building mesh' 103 self.mesh = Mesh(vertex_coordinates, triangles) 104 #self.mesh.check_integrity() # Time consuming 103 if verbose: print 'FitInterpolate: Building mesh' 104 self.mesh = Mesh(vertex_coordinates, triangles) 105 #self.mesh.check_integrity() # Time consuming 106 else: 107 self.mesh = None 105 108 else: 106 109 self.mesh = mesh 110 111 if self.mesh is not None: 112 if verbose: print 'FitInterpolate: Building quad tree' 113 #This stores indices of vertices 114 t0 = time.time() 115 #print "self.mesh.get_extent(absolute=True)", \ 116 #self.mesh.get_extent(absolute=True) 117 self.root = build_quadtree(self.mesh, 118 max_points_per_cell = max_vertices_per_cell) 119 #print "self.root",self.root.show() 107 120 108 if verbose: print 'FitInterpolate: Building quad tree' 109 # This stores indices of vertices 110 t0 = time.time() 111 #print "self.mesh.get_extent(absolute=True)", \ 112 #self.mesh.get_extent(absolute=True) 113 self.root = build_quadtree(self.mesh, 114 max_points_per_cell = max_vertices_per_cell) 115 #print "self.root",self.root.show() 116 117 build_quadtree_time = time.time()t0 118 set_last_triangle() 121 build_quadtree_time = time.time()t0 122 set_last_triangle() 119 123 120 124 def __repr__(self): 
anuga_core/source/anuga/fit_interpolate/interpolate.py
r5321 r5361 89 89 verbose=verbose, 90 90 max_vertices_per_cell=max_vertices_per_cell) 91 92 def interpolate_polyline(self, 93 f, 94 vertex_coordinates, 95 point_coordinates=None, 96 start_blocking_len=500000, 97 verbose=False): 98 99 if isinstance(point_coordinates, Geospatial_data): 100 point_coordinates = point_coordinates.get_data_points( \ 101 absolute = True) 102 103 from utilities.polygon import point_on_line,point_on_line_py 104 from Numeric import ones 105 z=ones(len(point_coordinates),Float) 106 107 msg='point coordinates are not given (interpolate.py)' 108 assert point_coordinates is not None, msg 109 msg='function value must be specified at every interpolation node' 110 assert f.shape[0]==vertex_coordinates.shape[0],msg 111 msg='Must define function value at one or more nodes' 112 assert f.shape[0]>0,msg 113 114 #print point_on_line_py(point_coordinates[3],[vertex_coordinates[0],vertex_coordinates[1]]) 115 #print vertex_coordinates[0],vertex_coordinates[1] 116 117 #print point_coordinates 118 #print vertex_coordinates 119 n=f.shape[0] 120 if n==1: 121 z=f*z 122 elif n>1: 123 index=0#index of segment on which last point was found 124 k=0#number of points on line 125 for i in range(len(point_coordinates)): 126 found = False 127 #find the segment the cetnroid lies on 128 #Start with the segment the last point was found on 129 #For n points there will be n1 segments 130 for j in range(n1): 131 #print 'searcing segment', index+j 132 #reached last segment look at first segment 133 if (index+j)==n2: 134 index=0 135 if point_on_line_py(point_coordinates[i],[vertex_coordinates[(j)+index],vertex_coordinates[(j+1)+index]]): 136 found=True 137 x0=vertex_coordinates[j][0];y0=vertex_coordinates[j][1] 138 x1=vertex_coordinates[j+1][0];y1=vertex_coordinates[j+1][1] 139 x2=point_coordinates[i][0];y2=point_coordinates[i][1] 140 141 segment_len=sqrt((x1x0)**2+(y1y0)**2) 142 dist=sqrt((x2x0)**2+(y2y0)**2) 143 z[i]=(f[j+1]f[j])/segment_len*dist+f[j] 144 #print 'element found on segment',j+index 145 index=j 146 break 147 148 149 if not found: 150 z[i]=0.0 151 #print 'point not on urs boundary' 152 else: 153 k+=1 154 #print 'number of midpoints on urs boundary',k 155 #print z 156 return z 91 157 92 158 # FIXME: What is a good start_blocking_len value? … … 529 595 vertex_coordinates = ensure_absolute(vertex_coordinates) 530 596 531 assert triangles is not None, 'Triangles array must be specified'532 triangles = ensure_numeric(triangles)533 self.spatial = True 597 if triangles is not None: 598 triangles = ensure_numeric(triangles) 599 self.spatial = True 534 600 535 601 # Thin timesteps if needed … … 557 623 # Precomputed spatial interpolation if requested 558 624 if interpolation_points is not None: 559 if self.spatial is False: 560 raise 'Triangles and vertex_coordinates must be specified' 561 625 #no longer true. sts files have spatial = True but 626 #if self.spatial is False: 627 # raise 'Triangles and vertex_coordinates must be specified' 628 # 562 629 try: 563 630 self.interpolation_points = interpolation_points = ensure_numeric(interpolation_points) … … 569 636 raise msg 570 637 571 572 # Check that all interpolation points fall within 573 # mesh boundary as defined by triangles and vertex_coordinates. 574 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh 575 from anuga.utilities.polygon import outside_polygon 576 577 # Create temporary mesh object from mesh info passed 578 # into this function. 579 mesh = Mesh(vertex_coordinates, triangles) 580 mesh_boundary_polygon = mesh.get_boundary_polygon() 581 582 583 indices = outside_polygon(interpolation_points, 584 mesh_boundary_polygon) 585 586 # Record result 587 #self.mesh_boundary_polygon = mesh_boundary_polygon 588 self.indices_outside_mesh = indices 589 590 # Report 591 if len(indices) > 0: 592 msg = 'Interpolation points in Interpolation function fall ' 593 msg += 'outside specified mesh. ' 594 msg += 'Offending points:\n' 595 out_interp_pts = [] 596 for i in indices: 597 msg += '%d: %s\n' %(i, interpolation_points[i]) 598 out_interp_pts.append(ensure_numeric(interpolation_points[i])) 599 600 if verbose is True: 601 import sys 602 if sys.platform == 'win32': 603 from anuga.utilities.polygon import plot_polygons 604 #out_interp_pts = take(interpolation_points,[indices]) 605 title = 'Interpolation points fall outside specified mesh' 606 plot_polygons([mesh_boundary_polygon,interpolation_points,out_interp_pts], 607 ['line','point','outside'],figname='points_boundary_out',label=title,verbose=verbose) 608 609 # Joaquim Luis suggested this as an Exception, so 610 # that the user can now what the problem is rather than 611 # looking for NaN's. However, NANs are handy as they can 612 # be ignored leaving good points for continued processing. 613 if verbose: 614 print msg 615 #raise Exception(msg) 638 if triangles is not None and vertex_coordinates is not None: 639 # Check that all interpolation points fall within 640 # mesh boundary as defined by triangles and vertex_coordinates. 641 from anuga.abstract_2d_finite_volumes.neighbour_mesh import Mesh 642 from anuga.utilities.polygon import outside_polygon 643 644 # Create temporary mesh object from mesh info passed 645 # into this function. 646 mesh = Mesh(vertex_coordinates, triangles) 647 mesh_boundary_polygon = mesh.get_boundary_polygon() 648 649 650 indices = outside_polygon(interpolation_points, 651 mesh_boundary_polygon) 652 653 # Record result 654 #self.mesh_boundary_polygon = mesh_boundary_polygon 655 self.indices_outside_mesh = indices 656 657 # Report 658 if len(indices) > 0: 659 msg = 'Interpolation points in Interpolation function fall ' 660 msg += 'outside specified mesh. ' 661 msg += 'Offending points:\n' 662 out_interp_pts = [] 663 for i in indices: 664 msg += '%d: %s\n' %(i, interpolation_points[i]) 665 out_interp_pts.append(ensure_numeric(interpolation_points[i])) 666 667 if verbose is True: 668 import sys 669 if sys.platform == 'win32': 670 from anuga.utilities.polygon import plot_polygons 671 #out_interp_pts = take(interpolation_points,[indices]) 672 title = 'Interpolation points fall outside specified mesh' 673 plot_polygons([mesh_boundary_polygon,interpolation_points,out_interp_pts], 674 ['line','point','outside'],figname='points_boundary_out',label=title,verbose=verbose) 675 676 # Joaquim Luis suggested this as an Exception, so 677 # that the user can now what the problem is rather than 678 # looking for NaN's. However, NANs are handy as they can 679 # be ignored leaving good points for continued processing. 680 if verbose: 681 print msg 682 #raise Exception(msg) 683 elif triangles is None and vertex_coordinates is not None:#jj 684 #Dealing with sts file 685 pass 686 else: 687 msg = 'Sww file function requires both triangles and vertex_coordinates. sts file file function requires the later.' 688 raise Exception(msg) 616 689 617 690 # Plot boundary and interpolation points … … 632 705 # Build interpolator 633 706 if verbose: 634 msg = 'Building interpolation matrix from source mesh ' 635 msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0], 636 triangles.shape[0]) 707 if triangles is not None and vertex_coordinates is not None: 708 msg = 'Building interpolation matrix from source mesh ' 709 msg += '(%d vertices, %d triangles)' %(vertex_coordinates.shape[0], 710 triangles.shape[0]) 711 elif triangles is None and vertex_coordinates is not None: 712 msg = 'Building interpolation matrix from source points' 713 637 714 print msg 638 715 … … 665 742 print ' quantity %s, size=%d' %(name, len(Q)) 666 743 667 # Interpolate 668 result = interpol.interpolate(Q, 669 point_coordinates=\ 670 self.interpolation_points, 671 verbose=False) # Don't clutter 744 # Interpolate 745 if triangles is not None and vertex_coordinates is not None: 746 result = interpol.interpolate(Q, 747 point_coordinates=\ 748 self.interpolation_points, 749 verbose=False) # Don't clutter 750 elif triangles is None and vertex_coordinates is not None: 751 result=interpol.interpolate_polyline(Q,vertex_coordinates,point_coordinates=self.interpolation_points) 672 752 673 753 #assert len(result), len(interpolation_points) 
anuga_core/source/anuga/shallow_water/test_data_manager.py
r5358 r5361 6289 6289 6290 6290 def test_urs2sts(self): 6291 """tide = 16292 time_step_count = 36293 time_step = 26294 lat_long = [[21.5,114.5],[21,114.5],[21,115]]6295 depth=206296 ha=26297 ua=56298 va=10 #ve added to take into account mux file format where south6299 # is positive.6300 """6301 #tide = 16302 6291 tide=0 6303 6292 time_step_count = 3 … … 6407 6396 self.delete_mux(files) 6408 6397 os.remove(sts_file) 6398 6399 def test_file_boundary_sts(self): 6400 from anuga.shallow_water import Domain 6401 from anuga.shallow_water import Reflective_boundary 6402 from anuga.shallow_water import Dirichlet_boundary 6403 from anuga.shallow_water import File_boundary 6404 from anuga.pmesh.mesh_interface import create_mesh_from_regions 6405 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance 6406 bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]] 6407 tide=0. 6408 time_step_count = 5 6409 time_step = 2 6410 lat_long_points =bounding_polygon[0:3] 6411 n=len(lat_long_points) 6412 first_tstep=ones(n,Int) 6413 last_tstep=(time_step_count)*ones(n,Int) 6414 gauge_depth=20*ones(n,Float) 6415 ha=2*ones((n,time_step_count),Float) 6416 ua=10*ones((n,time_step_count),Float) 6417 va=10*ones((n,time_step_count),Float) 6418 base_name, files = self.write_mux2(lat_long_points, 6419 time_step_count, time_step, 6420 first_tstep, last_tstep, 6421 depth=gauge_depth, 6422 ha=ha, 6423 ua=ua, 6424 va=va) 6425 6426 sts_file=base_name 6427 urs2sts(base_name,sts_file,mean_stage=tide,verbose=False) 6428 self.delete_mux(files) 6429 6430 #print 'start create mesh from regions' 6431 for i in range(len(bounding_polygon)): 6432 zone,bounding_polygon[i][0],bounding_polygon[i][1]=redfearn(bounding_polygon[i][0],bounding_polygon[i][1]) 6433 extent_res=1000000 6434 meshname = 'urs_test_mesh' + '.tsh' 6435 interior_regions=None 6436 boundary_tags={'ocean': [0,1], 'otherocean': [2,3,4]} 6437 create_mesh_from_regions(bounding_polygon,boundary_tags=boundary_tags, 6438 maximum_triangle_area=extent_res,filename=meshname, 6439 interior_regions=interior_regions,verbose=False) 6440 6441 domain_fbound = pmesh_to_domain_instance(meshname, Domain) 6442 domain_fbound.set_quantity('stage', tide) 6443 Bf = File_boundary(sts_file+'.sts', domain_fbound) 6444 Br = Reflective_boundary(domain_fbound) 6445 Bd=Dirichlet_boundary([2.0,220,220]) 6446 domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br}) 6447 finaltime=time_step*(time_step_count1) 6448 yieldstep=time_step 6449 temp_fbound=zeros(int(finaltime/yieldstep)+1,Float) 6450 i=0 6451 for t in domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime, 6452 skip_initial_step = False): 6453 temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2] 6454 i+=1 6455 6456 domain_drchlt = pmesh_to_domain_instance(meshname, Domain) 6457 domain_drchlt.set_quantity('stage', tide) 6458 Br = Reflective_boundary(domain_drchlt) 6459 Bd=Dirichlet_boundary([2.0,220,220]) 6460 domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br}) 6461 temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float) 6462 i=0 6463 for t in domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime, 6464 skip_initial_step = False): 6465 temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2] 6466 i+=1 6467 6468 assert temp_fboundtemp_drchlt<epsilon 6469 os.remove(sts_file+'.sts') 6470 os.remove(meshname) 6409 6471 6410 6472 def test_lon_lat2grid(self):
Note: See TracChangeset
for help on using the changeset viewer.