Changeset 5589


Ignore:
Timestamp:
Jul 31, 2008, 7:13:45 PM (11 years ago)
Author:
sexton
Message:

unit tests for comparing urs gauges to sts outputs

Location:
anuga_core/source/anuga/shallow_water
Files:
30 added
2 edited

Legend:

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

    r5586 r5589  
    51675167    # Create new file
    51685168    sts = Write_sts()
     5169
    51695170    sts.store_header(outfile,
    51705171                     times+starttime,
     
    56685669        outfile.institution = 'Geoscience Australia'
    56695670        outfile.description = description
    5670        
     5671
    56715672        try:
    56725673            revision_number = get_revision_number()
  • anuga_core/source/anuga/shallow_water/test_data_manager.py

    r5583 r5589  
    65416541        self.delete_mux(filesII)
    65426542        os.remove(sts_file)
    6543        
     6543
     6544    def test_urs2sts_individual_sources(self):   
     6545        """Test that individual sources compare to actual urs output
     6546           Test that the first recording time is the smallest
     6547           over waveheight, easting and northing velocity
     6548        """
     6549        from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
     6550            compress,zeros,fabs,take,size
     6551       
     6552        ordering_filename='thinned_bound_order_test.txt'
     6553        dir = 'urs_test_data'
     6554        sources = ['1-z.grd','2-z.grd','3-z.grd']
     6555        source_number = 0
     6556       
     6557        # make sts file for each source
     6558        for source in sources:
     6559            source_number += 1
     6560            urs_filenames = os.path.join(dir,source)
     6561            weights = [1.]
     6562            sts_name_out = 'test'
     6563           
     6564            urs2sts(urs_filenames,
     6565                    basename_out=sts_name_out,
     6566                    ordering_filename=ordering_filename,
     6567                    weights=weights,
     6568                    mean_stage=0.0,
     6569                    verbose=False)
     6570
     6571            # read in sts file for first source
     6572            fid = NetCDFFile(sts_name_out+'.sts', 'r')    #Open existing file for read
     6573            x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
     6574            y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
     6575            points=transpose(asarray([x.tolist(),y.tolist()]))
     6576            time=fid.variables['time'][:]+fid.starttime
     6577
     6578            # make sure start time is 9.5 which is the minimum time across all quantities, i.e. starting time
     6579
     6580            # get quantity data from sts file
     6581            quantity_names=['stage','xmomentum','ymomentum']
     6582            quantities = {}
     6583            for i, name in enumerate(quantity_names):
     6584                quantities[name] = fid.variables[name][:]
     6585
     6586            # for each station, compare urs2sts output to known urs output
     6587           
     6588            # taken manually from header files from urs output
     6589            time_start_z = [[10.0,11.5,13,14.5,17.7],
     6590                            [9.8,11.2,12.7,14.2,17.4],
     6591                            [9.5,10.9,12.4,13.9,17.1]]
     6592            delta_t = 0.1
     6593            #time_start_e = time_start_z
     6594            #time_start_n = time_start_z
     6595           
     6596            for j in range(len(x)):
     6597                index_start_urs = 0
     6598                index_end_urs = 0
     6599                index_start = 0
     6600                index_end = 0
     6601                count = 0
     6602                urs_file_name = 'z_'+str(source_number)+'_'+str(j)+'.csv'
     6603                dict = csv2array(os.path.join(dir,urs_file_name))
     6604                urs_stage = dict['urs_stage']
     6605                for i in range(len(urs_stage)):
     6606                    if urs_stage[i] == 0.0:
     6607                        index_start_urs = i+1
     6608                    if int(urs_stage[i]) == 99 and count <> 1:
     6609                        count +=1
     6610                        index_end_urs = i
     6611
     6612                if count == 0: index_end_urs = len(urs_stage)
     6613
     6614                start_times = time_start_z[source_number-1]
     6615                # check that actual start time matches header information
     6616                msg = 'start time from urs file is not the same as the header file for source %i and station %i' %(source_number,j)
     6617                assert allclose(index_start_urs,start_times[j]/delta_t), msg
     6618
     6619                index_start = 0
     6620                index_end = 0
     6621                count = 0
     6622                sts_stage = quantities['stage'][:,j]
     6623                for i in range(len(sts_stage)):
     6624                    if sts_stage[i] <> 0.0 and count <> 1:
     6625                        count += 1
     6626                        index_start = i
     6627                    if int(sts_stage[i]) == 99 and count <> 1:
     6628                        count += 1
     6629                        index_end = i
     6630
     6631                index_end = index_start + len(urs_stage[index_start_urs:index_end_urs])
     6632
     6633                # check that the lengths of urs stage and sts stage are the same
     6634                msg = 'Length of urs stage is not equal to the length of sts stage for station %i' %j
     6635                assert allclose(len(urs_stage[index_start_urs:index_end_urs]), len(sts_stage[index_start:index_end])), msg
     6636
     6637                #print urs_stage[index_start_urs:index_end_urs]
     6638                #print sts_stage[index_start:index_end]
     6639                # check that urs stage and sts stage are the same
     6640                msg = 'urs stage is not equal to sts stage for station %i' %j
     6641                max_error = max(urs_stage[index_start_urs:index_end_urs] - sts_stage[index_start:index_end])
     6642                min_error = min(urs_stage[index_start_urs:index_end_urs] - sts_stage[index_start:index_end])
     6643                assert max_error < 1.e-5, msg
     6644                assert abs(min_error) < 1.e5, msg
     6645               
     6646            fid.close()
     6647
     6648        #from os import sys
     6649        #sys.exit()
     6650        os.remove(sts_name_out+'.sts')
     6651
     6652    def test_urs2sts_combined_sources(self):   
     6653        """Test that combined sources compare to actual urs output
     6654           Test that the first recording time is the smallest
     6655           over waveheight, easting and northing velocity
     6656        """
     6657        from Numeric import asarray,transpose,sqrt,argmax,argmin,arange,Float,\
     6658            compress,zeros,fabs,take,size
     6659       
     6660        # make sts file for combined sources
     6661        weights = [1., 2., 3.]
     6662        ordering_filename='thinned_bound_order_test.txt'
     6663        dir = 'urs_test_data'
     6664        urs_filenames = [os.path.join(dir,'1-z.grd'),os.path.join(dir,'2-z.grd'),os.path.join(dir,'3-z.grd')]
     6665        sts_name_out = 'test'
     6666       
     6667        urs2sts(urs_filenames,
     6668                basename_out=sts_name_out,
     6669                ordering_filename=ordering_filename,
     6670                weights=weights,
     6671                mean_stage=0.0,
     6672                verbose=False)
     6673       
     6674        # read in sts file for third source
     6675        fid = NetCDFFile(sts_name_out+'.sts', 'r')    #Open existing file for read
     6676        x = fid.variables['x'][:]+fid.xllcorner   #x-coordinates of vertices
     6677        y = fid.variables['y'][:]+fid.yllcorner   #y-coordinates of vertices
     6678        points=transpose(asarray([x.tolist(),y.tolist()]))
     6679        time=fid.variables['time'][:]+fid.starttime
     6680
     6681        # get quantity data from sts file
     6682        quantity_names=['stage','xmomentum','ymomentum']
     6683        quantities = {}
     6684        for i, name in enumerate(quantity_names):
     6685            quantities[name] = fid.variables[name][:]
     6686
     6687        fid.close()
     6688       
     6689        os.remove(sts_name_out+'.sts')
    65446690       
    65456691    def test_urs2sts_ordering(self):
     
    97669912    #suite = unittest.makeSuite(Test_Data_Manager,'test_get_flow_through_cross_section_with_geo')
    97679913    #suite = unittest.makeSuite(Test_Data_Manager,'covered_')
     9914    #suite = unittest.makeSuite(Test_Data_Manager,'test_urs2sts_individual_sources')
    97689915
    97699916   
Note: See TracChangeset for help on using the changeset viewer.