Changeset 5418


Ignore:
Timestamp:
Jun 24, 2008, 11:42:50 AM (16 years ago)
Author:
jakeman
Message:

added ability to read in a file with sts indices to enable user to create boundary

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  
    50815081    outfile.close()
    50825082
     5083def 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
    50835146class Write_sww:
    50845147    from anuga.shallow_water.shallow_water_domain import Domain
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5416 r5418  
    65366536        os.remove(sts_file)
    65376537
    6538     def test_file_boundary_sts(self):
     6538    def test_file_boundary_stsI(self):
    65396539        from anuga.shallow_water import Domain
    65406540        from anuga.shallow_water import Reflective_boundary
     
    66096609        os.remove(meshname)
    66106610
    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         """
    66156618        from anuga.shallow_water import Domain
    66166619        from anuga.shallow_water import Reflective_boundary
     
    66616664        domain_fbound.set_quantity('stage', tide)
    66626665        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)
    66636785        Br = Reflective_boundary(domain_fbound)
    66646786        Bd=Dirichlet_boundary([2.0,220,-220])
     
    85548676    #suite = unittest.makeSuite(Test_Data_Manager,'test_export_gridII')
    85558677#    suite = unittest.makeSuite(Test_Data_Manager,'test_screen_catcher')
     8678<<<<<<< .mine
     8679    suite = unittest.makeSuite(Test_Data_Manager,'test')
     8680=======
    85568681    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2stsII')
     8682>>>>>>> .r5417
    85578683    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
    85588684    #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
Note: See TracChangeset for help on using the changeset viewer.