Ignore:
Timestamp:
Oct 10, 2006, 11:02:18 AM (17 years ago)
Author:
duncan
Message:

added urs2sww to convert urs tsunami info to an sww boundary

File:
1 edited

Legend:

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

    r3694 r3720  
    99import os
    1010from Scientific.IO.NetCDF import NetCDFFile
     11from struct import pack
    1112
    1213from anuga.shallow_water.data_manager import *
     
    1617# This is needed to run the tests of local functions
    1718import data_manager
    18 
     19from anuga.coordinate_transforms.redfearn import redfearn
    1920from anuga.coordinate_transforms.geo_reference import Geo_reference
    2021
     
    7172        self.F = bed
    7273
    73 
    74 
    75 
    7674        #Write A testfile (not realistic. Values aren't realistic)
    77 
    7875        self.test_MOST_file = 'most_small'
    7976
     
    49844981
    49854982        os.remove(csv_file)
     4983
     4984    #### TESTS FOR URS 2 SWW  ###     
     4985   
     4986    def create_mux(self, points_num=None):
     4987        # write all the mux stuff.
     4988        time_step_count = 3
     4989        time_step = 0.5
     4990       
     4991        longitudes = [150.66667, 150.83334, 151., 151.16667]
     4992        latitudes = [-34.5, -34.33333, -34.16667, -34]
     4993
     4994        if points_num == None:
     4995            points_num = len(longitudes) * len(latitudes)
     4996
     4997        lonlatdeps = []
     4998        quantities = ['HA','UA','VA']
     4999        mux_names = ['-z-mux','-e-mux','-n-mux']
     5000        quantities_init = [[],[],[]]
     5001        for i,lon in enumerate(longitudes):
     5002            for j,lat in enumerate(latitudes):
     5003                _ , e, n = redfearn(lat, lon)
     5004                lonlatdeps.append([lon, lat, n])
     5005                quantities_init[0].append(e) # HA
     5006                quantities_init[1].append(n ) # UA
     5007                quantities_init[2].append(e) # VA
     5008        #print "lonlatdeps",lonlatdeps
     5009        base_name = str(id(self))
     5010        files = []       
     5011        for i,q in enumerate(quantities):
     5012            quantities_init[i] = ensure_numeric(quantities_init[i])
     5013            #print "HA_init", HA_init
     5014            q_time = zeros((time_step_count, points_num), Float)
     5015            for time in range(time_step_count):
     5016                q_time[time,:] = quantities_init[i] #* time * 4
     5017           
     5018            #Write C files
     5019            columns = 3 # long, lat , depth
     5020            file = base_name + mux_names[i]
     5021            f = open(file, 'wb')
     5022            files.append(file)
     5023            f.write(pack('i',points_num))
     5024            f.write(pack('i',time_step_count))
     5025            f.write(pack('f',time_step))
     5026
     5027            #write lat/long info
     5028            for lonlatdep in lonlatdeps:
     5029                for float in lonlatdep:
     5030                    f.write(pack('f',float))
     5031                   
     5032            # Write quantity info
     5033            for time in  range(time_step_count):
     5034                for i in range(points_num):
     5035                    f.write(pack('f',q_time[time,i]))
     5036            f.close()
     5037        return base_name, files
     5038       
     5039   
     5040    def delete_mux(self, files):
     5041        for file in files:
     5042            os.remove(file)
     5043           
     5044    def test_urs2sww_test_fail(self):
     5045        points_num = -100
     5046        time_step_count = 45
     5047        time_step = -7
     5048        base_name = str(id(self))
     5049        files = []
     5050        quantities = ['HA','UA','VA']
     5051        mux_names = ['-z-mux','-e-mux','-n-mux']
     5052        for i,q in enumerate(quantities):
     5053            #Write C files
     5054            columns = 3 # long, lat , depth
     5055            file = base_name + mux_names[i]
     5056            f = open(file, 'wb')
     5057            files.append(file)
     5058            f.write(pack('i',points_num))
     5059            f.write(pack('i',time_step_count))
     5060            f.write(pack('f',time_step))
     5061
     5062            f.close()   
     5063        tide = 1
     5064        try:
     5065            urs2sww(base_name, remove_nc_files=True, mean_stage=tide)       
     5066        except ANUGAError:
     5067            pass
     5068        else:
     5069            self.delete_mux(files)
     5070            msg = 'Should have raised exception'
     5071            raise msg
     5072        sww_file = base_name + '.sww'
     5073       
     5074        self.delete_mux(files)
     5075       
     5076    def test_urs2sww(self):
     5077        tide = 1
     5078        base_name, files = self.create_mux()
     5079        urs2sww(base_name, remove_nc_files=True, mean_stage=tide)
     5080        sww_file = base_name + '.sww'
     5081       
     5082        #Let's interigate the sww file
     5083        fid = NetCDFFile(sww_file)
     5084
     5085        x = fid.variables['x'][:]
     5086        y = fid.variables['y'][:]
     5087        geo_reference = Geo_reference(NetCDFObject=fid)
     5088       
     5089        #Check that first coordinate is correctly represented       
     5090        #Work out the UTM coordinates for first point
     5091        zone, e, n = redfearn(-34.5, 150.66667)       
     5092       
     5093        assert allclose(geo_reference.get_absolute([[x[0],y[0]]]), [e,n])
     5094
     5095        #Check first value
     5096        stage = fid.variables['stage'][:]
     5097        xmomentum = fid.variables['xmomentum'][:]
     5098        ymomentum = fid.variables['ymomentum'][:]
     5099
     5100        #print ymomentum
     5101        assert allclose(stage[0,0], e +tide)  #Meters
     5102
     5103        #Check the momentums - ua
     5104        #momentum = velocity*(stage-elevation)
     5105        #momentum = velocity*(stage+elevation)
     5106        # -(-elevation) since elevation is inverted in mux files
     5107        # = n*(e+tide+n) based on how I'm writing these files
     5108        answer = n*(e+tide+n)
     5109        actual = xmomentum[0,0]
     5110        #print "answer",answer
     5111        #print "actual",actual
     5112        assert allclose(answer, actual)  #Meters
     5113       
     5114        fid.close()
     5115
     5116
     5117        self.delete_mux(files)
     5118        os.remove(sww_file)
     5119       
     5120    def bad_test_assuming_theres_a_mux_file(self):
     5121        # these mux files aren't in the repository, plus they are bad!
     5122        base_name = 'o-z-mux'
     5123        base_name = 'o-e-mux'
     5124        file_name = base_name + '.nc'
     5125        lonlatdep_numeric, lon, lat = _binary_c2nc(base_name, file_name, 'HA')
     5126       
     5127        #os.remove(file_name)
     5128       
     5129    def test_lon_lat2grid(self):
     5130        lonlatdep = [
     5131            [ 113.06700134  ,  -26.06669998 ,   0.        ] ,
     5132            [ 113.06700134  ,  -26.33329964 ,   0.        ] ,
     5133            [ 113.19999695  ,  -26.06669998 ,   0.        ] ,
     5134            [ 113.19999695  ,  -26.33329964 ,   0.        ] ]
     5135           
     5136        long, lat = lon_lat2grid(lonlatdep)
     5137
     5138        for i, result in enumerate(lat):
     5139            assert lonlatdep [i][1] == result
     5140        assert len(lat) == 2
     5141
     5142        for i, result in enumerate(long):
     5143            assert lonlatdep [i*2][0] == result
     5144        assert len(long) == 2
     5145
     5146    def test_lon_lat2grid_bad(self):
     5147        lonlatdep  = [
     5148            [ -26.06669998,  113.06700134,    1.        ],
     5149            [ -26.06669998 , 113.19999695 ,   2.        ],
     5150            [ -26.06669998 , 113.33300018,    3.        ],
     5151            [ -26.06669998 , 113.43299866   , 4.        ],
     5152            [ -26.20000076 , 113.06700134,    5.        ],
     5153            [ -26.20000076 , 113.19999695 ,   6.        ],
     5154            [ -26.20000076 , 113.33300018  ,  7.        ],
     5155            [ -26.20000076 , 113.43299866   , 8.        ],
     5156            [ -26.33329964 , 113.06700134,    9.        ],
     5157            [ -26.33329964 , 113.19999695 ,   10.        ],
     5158            [ -26.33329964 , 113.33300018  ,  11.        ],
     5159            [ -26.33329964 , 113.43299866 ,   12.        ],
     5160            [ -26.43330002 , 113.06700134 ,   13        ],
     5161            [ -26.43330002 , 113.19999695 ,   14.        ],
     5162            [ -26.43330002 , 113.33300018,    15.        ],
     5163            [ -26.43330002 , 113.43299866,    16.        ]]
     5164        try:
     5165            long, lat = lon_lat2grid(lonlatdep)
     5166        except AssertionError:
     5167            pass
     5168        else:
     5169            msg = 'Should have raised exception'
     5170            raise msg
     5171       
     5172    def test_lon_lat2gridII(self):
     5173        lonlatdep = [
     5174            [ 113.06700134  ,  -26.06669998 ,   0.        ] ,
     5175            [ 113.06700134  ,  -26.33329964 ,   0.        ] ,
     5176            [ 113.19999695  ,  -26.06669998 ,   0.        ] ,
     5177            [ 113.19999695  ,  -26.344329964 ,   0.        ] ]
     5178        try:
     5179            long, lat = lon_lat2grid(lonlatdep)
     5180        except AssertionError:
     5181            pass
     5182        else:
     5183            msg = 'Should have raised exception'
     5184            raise msg
     5185       
     5186    def trial_loading(self):
     5187        basename_in = 'karratha'
     5188        basename_out = basename_in
     5189        urs2sww(basename_in, basename_out)
     5190    #### END TESTS FOR URS 2 SWW  ###
     5191
    49865192       
    49875193#-------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.