Changeset 5418
- Timestamp:
- Jun 24, 2008, 11:42:50 AM (15 years ago)
- Location:
- anuga_core/source/anuga/shallow_water
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
anuga_core/source/anuga/shallow_water/data_manager.py
r5412 r5418 5081 5081 outfile.close() 5082 5082 5083 def create_sts_boundary(order_file,stsname): 5084 """ 5085 Create boundary segments from .sts file. Points can be stored in 5086 arbitrary order within the .sts file. The order in which the .sts points 5087 make up the boundary are given in order.txt file 5088 """ 5089 5090 if not order_file[-3:] == 'txt': 5091 msg = 'Order file must be a txt file' 5092 raise Exception, msg 5093 5094 try: 5095 fid=open(order_file,'r') 5096 file_header=fid.readline() 5097 lines=fid.readlines() 5098 fid.close() 5099 except: 5100 msg = 'Cannot open %s'%filename 5101 raise msg 5102 5103 header="index,longitude,latitude\n" 5104 5105 if not file_header==header: 5106 msg = 'File must contain header\n'+header+"\n" 5107 raise Exception, msg 5108 5109 try: 5110 fid = NetCDFFile(stsname+'.sts', 'r') 5111 except: 5112 msg = 'Cannot open %s'%filename+'.sts' 5113 5114 5115 xllcorner = fid.xllcorner[0] 5116 yllcorner = fid.yllcorner[0] 5117 #Points stored in sts file are normailsed to [xllcorner,yllcorner] but 5118 #we cannot assume that boundary polygon will be. At least the 5119 #additional points specified by the user after this function is called 5120 x = fid.variables['x'][:]+xllcorner 5121 y = fid.variables['y'][:]+yllcorner 5122 5123 x = reshape(x, (len(x),1)) 5124 y = reshape(y, (len(y),1)) 5125 sts_points = concatenate((x,y), axis=1) 5126 5127 xllcorner = fid.xllcorner[0] 5128 yllcorner = fid.yllcorner[0] 5129 5130 boundary_polygon=[] 5131 for line in lines: 5132 fields = line.split(',') 5133 index=int(fields[0]) 5134 zone,easting,northing=redfearn(float(fields[1]),float(fields[2])) 5135 if not zone == fid.zone[0]: 5136 msg = 'Inconsitent zones' 5137 raise Exception, msg 5138 if not allclose(array([easting,northing]),sts_points[index]): 5139 msg = "Points in sts file do not match those in the"+\ 5140 ".txt file spcifying the order of the boundary points" 5141 raise Exception, msg 5142 boundary_polygon.append(sts_points[index].tolist()) 5143 5144 return boundary_polygon 5145 5083 5146 class Write_sww: 5084 5147 from anuga.shallow_water.shallow_water_domain import Domain -
anuga_core/source/anuga/shallow_water/test_data_manager.py
r5416 r5418 6536 6536 os.remove(sts_file) 6537 6537 6538 def test_file_boundary_sts (self):6538 def test_file_boundary_stsI(self): 6539 6539 from anuga.shallow_water import Domain 6540 6540 from anuga.shallow_water import Reflective_boundary … … 6609 6609 os.remove(meshname) 6610 6610 6611 def test_file_boundary_sts2(self): 6612 # mux2 file has points not included in boundary 6613 # mux2 gauges are not stored with the same order as they are 6614 # found in bounding_polygon 6611 def test_file_boundary_stsII(self): 6612 """ mux2 file has points not included in boundary 6613 mux2 gauges are not stored with the same order as they are 6614 found in bounding_polygon. This does not matter as long as bounding 6615 polygon passed to file_function contains the mux2 points needed (in 6616 the correct order). 6617 """ 6615 6618 from anuga.shallow_water import Domain 6616 6619 from anuga.shallow_water import Reflective_boundary … … 6661 6664 domain_fbound.set_quantity('stage', tide) 6662 6665 Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=bounding_polygon) 6666 Br = Reflective_boundary(domain_fbound) 6667 Bd=Dirichlet_boundary([2.0,220,-220]) 6668 domain_fbound.set_boundary({'ocean': Bf,'otherocean': Br}) 6669 finaltime=time_step*(time_step_count-1) 6670 yieldstep=time_step 6671 temp_fbound=zeros(int(finaltime/yieldstep)+1,Float) 6672 i=0 6673 for t in domain_fbound.evolve(yieldstep=yieldstep,finaltime=finaltime, 6674 skip_initial_step = False): 6675 temp_fbound[i]=domain_fbound.quantities['stage'].centroid_values[2] 6676 i+=1 6677 6678 domain_drchlt = pmesh_to_domain_instance(meshname, Domain) 6679 domain_drchlt.set_quantity('stage', tide) 6680 Br = Reflective_boundary(domain_drchlt) 6681 Bd=Dirichlet_boundary([2.0,220,-220]) 6682 domain_drchlt.set_boundary({'ocean': Bd,'otherocean': Br}) 6683 temp_drchlt=zeros(int(finaltime/yieldstep)+1,Float) 6684 i=0 6685 for t in domain_drchlt.evolve(yieldstep=yieldstep,finaltime=finaltime, 6686 skip_initial_step = False): 6687 temp_drchlt[i]=domain_drchlt.quantities['stage'].centroid_values[2] 6688 i+=1 6689 6690 assert allclose(temp_fbound,temp_drchlt) 6691 os.remove(sts_file+'.sts') 6692 os.remove(meshname) 6693 6694 def test_file_boundary_stsIII(self): 6695 """Read correct points from order file 6696 """ 6697 from anuga.shallow_water import Domain 6698 from anuga.shallow_water import Reflective_boundary 6699 from anuga.shallow_water import Dirichlet_boundary 6700 from anuga.shallow_water import File_boundary 6701 from anuga.pmesh.mesh_interface import create_mesh_from_regions 6702 from anuga.abstract_2d_finite_volumes.pmesh2domain import pmesh_to_domain_instance 6703 6704 lat_long_points=[[6.01,97.0],[6.02,97.0],[6.05,96.9],[6.0,97.0]] 6705 bounding_polygon=[[6.0,97.0],[6.01,97.0],[6.02,97.0],[6.02,97.02],[6.00,97.02]] 6706 tide=0. 6707 time_step_count = 5 6708 time_step = 2 6709 n=len(lat_long_points) 6710 first_tstep=ones(n,Int) 6711 last_tstep=(time_step_count)*ones(n,Int) 6712 gauge_depth=20*ones(n,Float) 6713 ha=2*ones((n,time_step_count),Float) 6714 ua=10*ones((n,time_step_count),Float) 6715 va=-10*ones((n,time_step_count),Float) 6716 base_name, files = self.write_mux2(lat_long_points, 6717 time_step_count, time_step, 6718 first_tstep, last_tstep, 6719 depth=gauge_depth, 6720 ha=ha, 6721 ua=ua, 6722 va=va) 6723 6724 #Write order file 6725 file_handle, order_base_name = tempfile.mkstemp("") 6726 os.close(file_handle) 6727 os.remove(order_base_name) 6728 d="," 6729 order_file=order_base_name+'order.txt' 6730 fid=open(order_file,'w') 6731 #Write Header 6732 header="index,longitude,latitude\n" 6733 fid.write(header) 6734 indices=[3,0,1] 6735 for i in indices: 6736 line=str(i)+d+str(lat_long_points[i][0])+d+\ 6737 str(lat_long_points[i][1])+"\n" 6738 fid.write(line) 6739 fid.close() 6740 6741 sts_file=base_name 6742 urs2sts(base_name,sts_file,mean_stage=tide,verbose=False) 6743 self.delete_mux(files) 6744 6745 boundary_polygon = create_sts_boundary(order_file,base_name) 6746 6747 os.remove(order_file) 6748 6749 #Append the remaining part of the boundary polygon to be defined by 6750 #the user 6751 bounding_polygon_utm=[] 6752 for point in bounding_polygon: 6753 zone,easting,northing=redfearn(point[0],point[1]) 6754 bounding_polygon_utm.append([easting,northing]) 6755 6756 boundary_polygon.append(bounding_polygon_utm[3]) 6757 boundary_polygon.append(bounding_polygon_utm[4]) 6758 6759 plot=False 6760 if plot: 6761 from pylab import plot,show,axis 6762 boundary_polygon=ensure_numeric(boundary_polygon) 6763 bounding_polygon_utm=ensure_numeric(bounding_polygon_utm) 6764 #plot(lat_long_points[:,0],lat_long_points[:,1],'o') 6765 plot(boundary_polygon[:,0], boundary_polygon[:,1]) 6766 plot(bounding_polygon_utm[:,0],bounding_polygon_utm[:,1]) 6767 show() 6768 6769 assert allclose(bounding_polygon_utm,boundary_polygon) 6770 6771 6772 extent_res=1000000 6773 meshname = 'urs_test_mesh' + '.tsh' 6774 interior_regions=None 6775 boundary_tags={'ocean': [0,4], 'otherocean': [1,2,3]} 6776 #have to change boundary tags from last example because now bounding 6777 #polygon starts in different place. 6778 create_mesh_from_regions(boundary_polygon,boundary_tags=boundary_tags, 6779 maximum_triangle_area=extent_res,filename=meshname, 6780 interior_regions=interior_regions,verbose=False) 6781 6782 domain_fbound = pmesh_to_domain_instance(meshname, Domain) 6783 domain_fbound.set_quantity('stage', tide) 6784 Bf = File_boundary(sts_file+'.sts', domain_fbound, boundary_polygon=boundary_polygon) 6663 6785 Br = Reflective_boundary(domain_fbound) 6664 6786 Bd=Dirichlet_boundary([2.0,220,-220]) … … 8554 8676 #suite = unittest.makeSuite(Test_Data_Manager,'test_export_gridII') 8555 8677 # suite = unittest.makeSuite(Test_Data_Manager,'test_screen_catcher') 8678 <<<<<<< .mine 8679 suite = unittest.makeSuite(Test_Data_Manager,'test') 8680 ======= 8556 8681 #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2stsII') 8682 >>>>>>> .r5417 8557 8683 #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo') 8558 8684 #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
Note: See TracChangeset
for help on using the changeset viewer.