Ignore:
Timestamp:
Jul 18, 2008, 4:59:46 PM (14 years ago)
Author:
ole
Message:

Work on permutations in urs_ext.c
It does not work yet, but the array is being passed in and the algorithm
has been described.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5531 r5537  
    61076107        return base_name, files
    61086108
    6109     def test_read_mux2_pyI(self):
    6110         """Constant stage,momentum at each gauge
     6109    def test_urs2sts_read_mux2_pyI(self):
     6110        """test_urs2sts_read_mux2_pyI(self):
     6111        Constant stage,momentum at each gauge
    61116112        """
    61126113        tide = 1
     
    61596160        assert allclose(yvelocity,va)
    61606161
    6161     def test_read_mux2_pyII(self):
     6162    def test_urs2sts_read_mux2_pyII(self):
    61626163        """Spatially varing stage
    61636164        """
     
    62176218        assert allclose(yvelocity,va)
    62186219
    6219     def test_read_mux2_pyIII(self):
     6220    def test_urs2sts_read_mux2_pyIII(self):
    62206221        """Varying start and finsh times
    62216222        """
     
    65876588        permutation = [3,0,2]
    65886589
    6589         # Do it wrongly at first and check that exception is being raised
     6590        _, ordering_filename = tempfile.mkstemp('')
     6591        order_fid = open(ordering_filename, 'w') 
     6592        order_fid.write('index, longitude, latitude\n')
     6593        for index in permutation:
     6594            order_fid.write('%d, %f, %f\n' %(index,
     6595                                             lat_long_points[index][1],
     6596                                             lat_long_points[index][0]))
     6597        order_fid.close()
     6598       
     6599           
     6600
     6601                                               
     6602        # Call urs2sts with multiple mux files
     6603        urs2sts([base_nameI, base_nameII],
     6604                basename_out=base_nameI,
     6605                ordering_filename=ordering_filename,
     6606                weights=[1.0,1.0],
     6607                mean_stage=tide,
     6608                verbose=False)
     6609
     6610        # now I want to check the sts file ...
     6611        sts_file = base_nameI + '.sts'
     6612
     6613        #Let's interrogate the sts file
     6614        # Note, the sts info is not gridded.  It is point data.
     6615        fid = NetCDFFile(sts_file)
     6616
     6617        # Make x and y absolute
     6618        x = fid.variables['x'][:]
     6619        y = fid.variables['y'][:]
     6620
     6621        geo_reference = Geo_reference(NetCDFObject=fid)
     6622        points = geo_reference.get_absolute(map(None, x, y))
     6623        points = ensure_numeric(points)
     6624
     6625        x = points[:,0]
     6626        y = points[:,1]
     6627
     6628        for i, index in enumerate(permutation):
     6629            # Check that STS points are stored in the correct order
     6630           
     6631            # Work out the UTM coordinates sts point i
     6632            zone, e, n = redfearn(lat_long_points[index][0],
     6633                                  lat_long_points[index][1])             
     6634           
     6635            assert allclose([x[i],y[i]], [e,n])
     6636           
     6637                       
     6638        # Check the time vector
     6639        times = fid.variables['time'][:]
     6640
     6641        times_actual = []
     6642        for i in range(time_step_count):
     6643            times_actual.append(time_step * i)
     6644
     6645        assert allclose(ensure_numeric(times),
     6646                        ensure_numeric(times_actual))
     6647                       
     6648
     6649        # Check sts values
     6650        stage = fid.variables['stage'][:]
     6651        xmomentum = fid.variables['xmomentum'][:]
     6652        ymomentum = fid.variables['ymomentum'][:]
     6653        elevation = fid.variables['elevation'][:]
     6654
     6655        # Set original data used to write mux file to be zero when gauges are
     6656        # not recdoring
     6657       
     6658        ha[0][0]=0.0
     6659        ha[0][time_step_count-1]=0.0
     6660        ha[2][0]=0.0
     6661        ua[0][0]=0.0
     6662        ua[0][time_step_count-1]=0.0
     6663        ua[2][0]=0.0
     6664        va[0][0]=0.0
     6665        va[0][time_step_count-1]=0.0
     6666        va[2][0]=0.0;
     6667
     6668       
     6669        # The stage stored in the .sts file should be the sum of the stage
     6670        # in the two mux2 files because both have weights = 1. In this case
     6671        # the mux2 files are the same so stage == 2.0 * ha
     6672        #print 2.0*transpose(ha) - stage
     6673       
     6674        ha_permutation = take(ha, permutation)
     6675        ua_permutation = take(ua, permutation)         
     6676        va_permutation = take(va, permutation)                 
     6677        gauge_depth_permutation = take(gauge_depth, permutation)                         
     6678
     6679       
     6680        assert allclose(2.0*transpose(ha_permutation), stage)  # Meters
     6681
     6682        #Check the momentums - ua
     6683        #momentum = velocity*(stage-elevation)
     6684        # elevation = - depth
     6685        #momentum = velocity_ua *(stage+depth)
     6686
     6687        depth=zeros((len(lat_long_points),time_step_count),Float)
     6688        for i in range(len(lat_long_points)):
     6689            depth[i]=gauge_depth[i]+tide+2.0*ha[i]
     6690            #2.0*ha necessary because using two files with weights=1 are used
     6691           
     6692        depth_permutation = take(depth, permutation)                     
     6693       
     6694
     6695        # The xmomentum stored in the .sts file should be the sum of the ua
     6696        # in the two mux2 files multiplied by the depth.
     6697        assert allclose(2.0*transpose(ua_permutation*depth_permutation), xmomentum)
     6698
     6699        #Check the momentums - va
     6700        #momentum = velocity*(stage-elevation)
     6701        # elevation = - depth
     6702        #momentum = velocity_va *(stage+depth)
     6703
     6704        # The ymomentum stored in the .sts file should be the sum of the va
     6705        # in the two mux2 files multiplied by the depth.
     6706        assert allclose(2.0*transpose(va_permutation*depth_permutation), ymomentum)
     6707
     6708        # check the elevation values.
     6709        # -ve since urs measures depth, sww meshers height,
     6710        assert allclose(-gauge_depth_permutation, elevation)  #Meters
     6711
     6712        fid.close()
     6713        self.delete_mux(filesI)
     6714        self.delete_mux(filesII)
     6715        os.remove(sts_file)
     6716       
     6717
     6718       
     6719       
     6720       
     6721    def test_urs2sts_ordering_exception(self):
     6722        """Test that inconsistent lats and lons in ordering file are caught.
     6723        """
     6724       
     6725        tide=0
     6726        time_step_count = 3
     6727        time_step = 2
     6728        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
     6729        n=len(lat_long_points)
     6730        first_tstep=ones(n,Int)
     6731        first_tstep[0]+=1
     6732        first_tstep[2]+=1
     6733        last_tstep=(time_step_count)*ones(n,Int)
     6734        last_tstep[0]-=1
     6735
     6736        gauge_depth=20*ones(n,Float)
     6737        ha=2*ones((n,time_step_count),Float)
     6738        ha[0]=arange(0,time_step_count)
     6739        ha[1]=arange(time_step_count,2*time_step_count)
     6740        ha[2]=arange(2*time_step_count,3*time_step_count)
     6741        ha[3]=arange(3*time_step_count,4*time_step_count)
     6742        ua=5*ones((n,time_step_count),Float)
     6743        va=-10*ones((n,time_step_count),Float)
     6744
     6745        # Create two identical mux files to be combined by urs2sts
     6746        base_nameI, filesI = self.write_mux2(lat_long_points,
     6747                                             time_step_count, time_step,
     6748                                             first_tstep, last_tstep,
     6749                                             depth=gauge_depth,
     6750                                             ha=ha,
     6751                                             ua=ua,
     6752                                             va=va)
     6753
     6754        base_nameII, filesII = self.write_mux2(lat_long_points,
     6755                                               time_step_count, time_step,
     6756                                               first_tstep, last_tstep,
     6757                                               depth=gauge_depth,
     6758                                               ha=ha,
     6759                                               ua=ua,
     6760                                               va=va)
     6761
     6762                                               
     6763        # Create ordering file
     6764        permutation = [3,0,2]
     6765
     6766        # Do it wrongly and check that exception is being raised
    65906767        _, ordering_filename = tempfile.mkstemp('')
    65916768        order_fid = open(ordering_filename, 'w') 
     
    66136790
    66146791       
    6615            
    6616         # Then do it right
    6617         _, ordering_filename = tempfile.mkstemp('')
    6618         order_fid = open(ordering_filename, 'w') 
    6619         order_fid.write('index, longitude, latitude\n')
    6620         for index in permutation:
    6621             order_fid.write('%d, %f, %f\n' %(index,
    6622                                              lat_long_points[index][1],
    6623                                              lat_long_points[index][0]))
    6624         order_fid.close()
    6625        
    6626            
    6627 
    6628                                                
    6629         # Call urs2sts with multiple mux files
    6630         urs2sts([base_nameI, base_nameII],
    6631                 basename_out=base_nameI,
    6632                 ordering_filename=ordering_filename,
    6633                 weights=[1.0,1.0],
    6634                 mean_stage=tide,
    6635                 verbose=False)
    6636 
    6637         # now I want to check the sts file ...
    6638         sts_file = base_nameI + '.sts'
    6639 
    6640         #Let's interrogate the sts file
    6641         # Note, the sts info is not gridded.  It is point data.
    6642         fid = NetCDFFile(sts_file)
    6643 
    6644         # Make x and y absolute
    6645         x = fid.variables['x'][:]
    6646         y = fid.variables['y'][:]
    6647 
    6648         geo_reference = Geo_reference(NetCDFObject=fid)
    6649         points = geo_reference.get_absolute(map(None, x, y))
    6650         points = ensure_numeric(points)
    6651 
    6652         x = points[:,0]
    6653         y = points[:,1]
    6654 
    6655         for i, index in enumerate(permutation):
    6656             # Check that STS points are stored in the correct order
    6657            
    6658             # Work out the UTM coordinates sts point i
    6659             zone, e, n = redfearn(lat_long_points[index][0],
    6660                                   lat_long_points[index][1])             
    6661            
    6662             assert allclose([x[i],y[i]], [e,n])
    6663            
    6664                        
    6665         # Check the time vector
    6666         times = fid.variables['time'][:]
    6667 
    6668         times_actual = []
    6669         for i in range(time_step_count):
    6670             times_actual.append(time_step * i)
    6671 
    6672         assert allclose(ensure_numeric(times),
    6673                         ensure_numeric(times_actual))
    6674                        
    6675 
    6676         # Check sts values
    6677         stage = fid.variables['stage'][:]
    6678         xmomentum = fid.variables['xmomentum'][:]
    6679         ymomentum = fid.variables['ymomentum'][:]
    6680         elevation = fid.variables['elevation'][:]
    6681 
    6682         # Set original data used to write mux file to be zero when gauges are
    6683         # not recdoring
    6684        
    6685         ha[0][0]=0.0
    6686         ha[0][time_step_count-1]=0.0
    6687         ha[2][0]=0.0
    6688         ua[0][0]=0.0
    6689         ua[0][time_step_count-1]=0.0
    6690         ua[2][0]=0.0
    6691         va[0][0]=0.0
    6692         va[0][time_step_count-1]=0.0
    6693         va[2][0]=0.0;
    6694 
    6695        
    6696         # The stage stored in the .sts file should be the sum of the stage
    6697         # in the two mux2 files because both have weights = 1. In this case
    6698         # the mux2 files are the same so stage == 2.0 * ha
    6699         #print 2.0*transpose(ha) - stage
    6700        
    6701         ha_permutation = take(ha, permutation)
    6702         ua_permutation = take(ua, permutation)         
    6703         va_permutation = take(va, permutation)                 
    6704         gauge_depth_permutation = take(gauge_depth, permutation)                         
    6705 
    6706        
    6707         assert allclose(2.0*transpose(ha_permutation), stage)  # Meters
    6708 
    6709         #Check the momentums - ua
    6710         #momentum = velocity*(stage-elevation)
    6711         # elevation = - depth
    6712         #momentum = velocity_ua *(stage+depth)
    6713 
    6714         depth=zeros((len(lat_long_points),time_step_count),Float)
    6715         for i in range(len(lat_long_points)):
    6716             depth[i]=gauge_depth[i]+tide+2.0*ha[i]
    6717             #2.0*ha necessary because using two files with weights=1 are used
    6718            
    6719         depth_permutation = take(depth, permutation)                     
    6720        
    6721 
    6722         # The xmomentum stored in the .sts file should be the sum of the ua
    6723         # in the two mux2 files multiplied by the depth.
    6724         assert allclose(2.0*transpose(ua_permutation*depth_permutation), xmomentum)
    6725 
    6726         #Check the momentums - va
    6727         #momentum = velocity*(stage-elevation)
    6728         # elevation = - depth
    6729         #momentum = velocity_va *(stage+depth)
    6730 
    6731         # The ymomentum stored in the .sts file should be the sum of the va
    6732         # in the two mux2 files multiplied by the depth.
    6733         assert allclose(2.0*transpose(va_permutation*depth_permutation), ymomentum)
    6734 
    6735         # check the elevation values.
    6736         # -ve since urs measures depth, sww meshers height,
    6737         assert allclose(-gauge_depth_permutation, elevation)  #Meters
    6738 
    6739         fid.close()
    6740         self.delete_mux(filesI)
    6741         self.delete_mux(filesII)
    6742         os.remove(sts_file)
    6743        
    6744 
    6745     def test_file_boundary_stsI(self):
     6792       
     6793       
     6794    def test_urs2sts_file_boundary_stsI(self):
    67466795        from anuga.shallow_water import Domain
    67476796        from anuga.shallow_water import Reflective_boundary
     
    68186867        os.remove(meshname)
    68196868
    6820     def test_file_boundary_stsII(self):
     6869    def test_urs2sts_file_boundary_stsII(self):
    68216870        """ mux2 file has points not included in boundary
    68226871         mux2 gauges are not stored with the same order as they are
     
    69016950        os.remove(meshname)
    69026951
    6903     def test_file_boundary_stsIII(self):
     6952    def test_urs2sts_file_boundary_stsIII_ordering(self):
    69046953        """Read correct points from ordering file
    69056954        """
     
    88918940if __name__ == "__main__":
    88928941
    8893     suite = unittest.makeSuite(Test_Data_Manager,'test')
    8894     #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts0')
    8895     #suite = unittest.makeSuite(Test_Data_Manager,'test_export_gridII')
     8942    #suite = unittest.makeSuite(Test_Data_Manager,'test')
     8943    suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_read_mux2_pyI')
     8944    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts')
    88968945    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
    88978946    #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
     
    89068955    else:
    89078956        pass
    8908     runner = unittest.TextTestRunner() # verbosity=2)
     8957    runner = unittest.TextTestRunner(verbosity=2)
    89098958    runner.run(suite)
    89108959
Note: See TracChangeset for help on using the changeset viewer.