Ignore:
Timestamp:
Jul 4, 2008, 11:48:37 AM (16 years ago)
Author:
ole
Message:

John and Ole implemented urs2sts using an ordering file creating smaller and pre-ordered sts files for use with boundary conditions.

File:
1 edited

Legend:

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

    r5442 r5462  
    63196319                                      va=va)
    63206320
    6321         urs2sts(base_name,mean_stage=tide,verbose=False)
     6321        urs2sts(base_name,
     6322                basename_out=base_name,
     6323                mean_stage=tide,verbose=False)
    63226324
    63236325        # now I want to check the sts file ...
     
    64416443
    64426444        # Call urs2sts with multiple mux files
    6443         urs2sts([base_nameI, base_nameII], weights=[1.0,1.0],
     6445        urs2sts([base_nameI, base_nameII],
     6446                basename_out=base_nameI,
     6447                weights=[1.0,1.0],
    64446448                mean_stage=tide,
    64456449                verbose=False)
     
    65346538        self.delete_mux(filesII)
    65356539        os.remove(sts_file)
     6540       
     6541       
     6542    def test_urs2sts_ordering(self):
     6543        """Test multiple sources with ordering file
     6544        """
     6545       
     6546        tide=0
     6547        time_step_count = 3
     6548        time_step = 2
     6549        lat_long_points =[(-21.5,114.5),(-21,114.5),(-21.5,115), (-21.,115.)]
     6550        n=len(lat_long_points)
     6551        first_tstep=ones(n,Int)
     6552        first_tstep[0]+=1
     6553        first_tstep[2]+=1
     6554        last_tstep=(time_step_count)*ones(n,Int)
     6555        last_tstep[0]-=1
     6556
     6557        gauge_depth=20*ones(n,Float)
     6558        ha=2*ones((n,time_step_count),Float)
     6559        ha[0]=arange(0,time_step_count)
     6560        ha[1]=arange(time_step_count,2*time_step_count)
     6561        ha[2]=arange(2*time_step_count,3*time_step_count)
     6562        ha[3]=arange(3*time_step_count,4*time_step_count)
     6563        ua=5*ones((n,time_step_count),Float)
     6564        va=-10*ones((n,time_step_count),Float)
     6565
     6566        # Create two identical mux files to be combined by urs2sts
     6567        base_nameI, filesI = self.write_mux2(lat_long_points,
     6568                                             time_step_count, time_step,
     6569                                             first_tstep, last_tstep,
     6570                                             depth=gauge_depth,
     6571                                             ha=ha,
     6572                                             ua=ua,
     6573                                             va=va)
     6574
     6575        base_nameII, filesII = self.write_mux2(lat_long_points,
     6576                                               time_step_count, time_step,
     6577                                               first_tstep, last_tstep,
     6578                                               depth=gauge_depth,
     6579                                               ha=ha,
     6580                                               ua=ua,
     6581                                               va=va)
     6582
     6583                                               
     6584        # Create ordering file
     6585        permutation = [3,0,2]
     6586
     6587        # Do it wrongly at first and check that exception is being raised
     6588        _, ordering_filename = tempfile.mkstemp('')
     6589        order_fid = open(ordering_filename, 'w') 
     6590        order_fid.write('index, longitude, latitude\n')
     6591        for index in permutation:
     6592            order_fid.write('%d, %f, %f\n' %(index,
     6593                                             lat_long_points[index][0],
     6594                                             lat_long_points[index][1]))
     6595        order_fid.close()
     6596       
     6597        try:
     6598            urs2sts([base_nameI, base_nameII],
     6599                    basename_out=base_nameI,
     6600                    ordering_filename=ordering_filename,
     6601                    weights=[1.0,1.0],
     6602                    mean_stage=tide,
     6603                    verbose=False) 
     6604            os.remove(ordering_filename)           
     6605        except:
     6606            pass
     6607        else:
     6608            msg = 'Should have caught wrong lat longs'
     6609            raise Exception, msg
     6610
     6611
     6612       
     6613           
     6614        # Then do it right
     6615        _, ordering_filename = tempfile.mkstemp('')
     6616        order_fid = open(ordering_filename, 'w') 
     6617        order_fid.write('index, longitude, latitude\n')
     6618        for index in permutation:
     6619            order_fid.write('%d, %f, %f\n' %(index,
     6620                                             lat_long_points[index][1],
     6621                                             lat_long_points[index][0]))
     6622        order_fid.close()
     6623       
     6624           
     6625
     6626                                               
     6627        # Call urs2sts with multiple mux files
     6628        urs2sts([base_nameI, base_nameII],
     6629                basename_out=base_nameI,
     6630                ordering_filename=ordering_filename,
     6631                weights=[1.0,1.0],
     6632                mean_stage=tide,
     6633                verbose=False)
     6634
     6635        # now I want to check the sts file ...
     6636        sts_file = base_nameI + '.sts'
     6637
     6638        #Let's interrogate the sts file
     6639        # Note, the sts info is not gridded.  It is point data.
     6640        fid = NetCDFFile(sts_file)
     6641
     6642        # Make x and y absolute
     6643        x = fid.variables['x'][:]
     6644        y = fid.variables['y'][:]
     6645
     6646        geo_reference = Geo_reference(NetCDFObject=fid)
     6647        points = geo_reference.get_absolute(map(None, x, y))
     6648        points = ensure_numeric(points)
     6649
     6650        x = points[:,0]
     6651        y = points[:,1]
     6652
     6653        for i, index in enumerate(permutation):
     6654            # Check that STS points are stored in the correct order
     6655           
     6656            # Work out the UTM coordinates sts point i
     6657            zone, e, n = redfearn(lat_long_points[index][0],
     6658                                  lat_long_points[index][1])             
     6659           
     6660            assert allclose([x[i],y[i]], [e,n])
     6661           
     6662                       
     6663        # Check the time vector
     6664        times = fid.variables['time'][:]
     6665
     6666        times_actual = []
     6667        for i in range(time_step_count):
     6668            times_actual.append(time_step * i)
     6669
     6670        assert allclose(ensure_numeric(times),
     6671                        ensure_numeric(times_actual))
     6672                       
     6673
     6674        # Check sts values
     6675        stage = fid.variables['stage'][:]
     6676        xmomentum = fid.variables['xmomentum'][:]
     6677        ymomentum = fid.variables['ymomentum'][:]
     6678        elevation = fid.variables['elevation'][:]
     6679
     6680        # Set original data used to write mux file to be zero when gauges are
     6681        # not recdoring
     6682       
     6683        ha[0][0]=0.0
     6684        ha[0][time_step_count-1]=0.0
     6685        ha[2][0]=0.0
     6686        ua[0][0]=0.0
     6687        ua[0][time_step_count-1]=0.0
     6688        ua[2][0]=0.0
     6689        va[0][0]=0.0
     6690        va[0][time_step_count-1]=0.0
     6691        va[2][0]=0.0;
     6692
     6693       
     6694        # The stage stored in the .sts file should be the sum of the stage
     6695        # in the two mux2 files because both have weights = 1. In this case
     6696        # the mux2 files are the same so stage == 2.0 * ha
     6697        #print 2.0*transpose(ha) - stage
     6698       
     6699        ha_permutation = take(ha, permutation)
     6700        ua_permutation = take(ua, permutation)         
     6701        va_permutation = take(va, permutation)                 
     6702        gauge_depth_permutation = take(gauge_depth, permutation)                         
     6703
     6704       
     6705        assert allclose(2.0*transpose(ha_permutation), stage)  # Meters
     6706
     6707        #Check the momentums - ua
     6708        #momentum = velocity*(stage-elevation)
     6709        # elevation = - depth
     6710        #momentum = velocity_ua *(stage+depth)
     6711
     6712        depth=zeros((len(lat_long_points),time_step_count),Float)
     6713        for i in range(len(lat_long_points)):
     6714            depth[i]=gauge_depth[i]+tide+2.0*ha[i]
     6715            #2.0*ha necessary because using two files with weights=1 are used
     6716           
     6717        depth_permutation = take(depth, permutation)                     
     6718       
     6719
     6720        # The xmomentum stored in the .sts file should be the sum of the ua
     6721        # in the two mux2 files multiplied by the depth.
     6722        assert allclose(2.0*transpose(ua_permutation*depth_permutation), xmomentum)
     6723
     6724        #Check the momentums - va
     6725        #momentum = velocity*(stage-elevation)
     6726        # elevation = - depth
     6727        #momentum = velocity_va *(stage+depth)
     6728
     6729        # The ymomentum stored in the .sts file should be the sum of the va
     6730        # in the two mux2 files multiplied by the depth.
     6731        assert allclose(2.0*transpose(va_permutation*depth_permutation), ymomentum)
     6732
     6733        # check the elevation values.
     6734        # -ve since urs measures depth, sww meshers height,
     6735        assert allclose(-gauge_depth_permutation, elevation)  #Meters
     6736
     6737        fid.close()
     6738        self.delete_mux(filesI)
     6739        self.delete_mux(filesII)
     6740        os.remove(sts_file)
     6741       
    65366742
    65376743    def test_file_boundary_stsI(self):
     
    65636769
    65646770        sts_file=base_name
    6565         urs2sts(base_name,sts_file,mean_stage=tide,verbose=False)
     6771        urs2sts(base_name,
     6772                sts_file,
     6773                mean_stage=tide,verbose=False)
    65666774        self.delete_mux(files)
    65676775
Note: See TracChangeset for help on using the changeset viewer.