Changeset 3720


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

added urs2sww to convert urs tsunami info to an sww boundary

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

Legend:

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

    r3694 r3720  
    6060import csv
    6161import os
     62from struct import unpack
     63import array as p_array
    6264
    6365from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \
    64      searchsorted, zeros, allclose, around
     66     searchsorted, zeros, allclose, around, reshape
     67from Scientific.IO.NetCDF import NetCDFFile
    6568
    6669
     
    6871from anuga.geospatial_data.geospatial_data import Geospatial_data
    6972from anuga.config import minimum_storable_height as default_minimum_storable_height
     73from anuga.utilities.numerical_tools import ensure_numeric
     74from anuga.caching.caching import myhash
     75from anuga.utilities.anuga_exceptions import ANUGAError
    7076
    7177def make_filename(s):
     
    40434049    return
    40444050
     4051    ####  URS 2 SWW  ###
     4052
     4053lon_name = 'LON'
     4054lat_name = 'LAT'
     4055time_name = 'TIME'
     4056precision = Float # So if we want to change the precision its done here       
     4057class Write_nc:
     4058    """
     4059    Write an nc file.
     4060   
     4061    Note, this should be checked to meet cdc netcdf conventions for gridded
     4062    data. http://www.cdc.noaa.gov/cdc/conventions/cdc_netcdf_standard.shtml
     4063   
     4064    """
     4065    def __init__(self,
     4066                 quantity_name,
     4067                 file_name,
     4068                 time_step_count,
     4069                 time_step,
     4070                 lon,
     4071                 lat):
     4072        """
     4073        time_step_count is the number of time steps.
     4074        time_step is the time step size
     4075       
     4076        pre-condition: quantity_name must be 'HA' 'UA'or 'VA'.
     4077        """
     4078        self.quantity_name = quantity_name
     4079        quantity_units= {'HA':'METERS',
     4080                              'UA':'METERS/SECOND',
     4081                              'VA':'METERS/SECOND'}
     4082       
     4083        #self.file_name = file_name
     4084        self.time_step_count = time_step_count
     4085        self.time_step = time_step
     4086
     4087        # NetCDF file definition
     4088        self.outfile = NetCDFFile(file_name, 'w')
     4089        outfile = self.outfile       
     4090
     4091        #Create new file
     4092        nc_lon_lat_header(outfile, lon, lat)
     4093   
     4094        # TIME
     4095        outfile.createDimension(time_name, None)
     4096        outfile.createVariable(time_name, precision, (time_name,))
     4097
     4098        #QUANTITY
     4099        outfile.createVariable(self.quantity_name, precision,
     4100                               (time_name, lat_name, lon_name))
     4101        outfile.variables[self.quantity_name].missing_value=-1.e+034
     4102        outfile.variables[self.quantity_name].units= \
     4103                                 quantity_units[self.quantity_name]
     4104        outfile.variables[lon_name][:]= ensure_numeric(lon)
     4105        outfile.variables[lat_name][:]= ensure_numeric(lat)
     4106
     4107        #Assume no one will be wanting to read this, while we are writing
     4108        #outfile.close()
     4109       
     4110    def store_timestep(self, quantity_slice):
     4111        """
     4112        quantity_slice is the data to be stored at this time step
     4113        """
     4114       
     4115        outfile = self.outfile
     4116       
     4117        # Get the variables
     4118        time = outfile.variables[time_name]
     4119        quantity = outfile.variables[self.quantity_name]
     4120           
     4121        i = len(time)
     4122
     4123        #Store time
     4124        time[i] = i*self.time_step #self.domain.time
     4125        quantity[i,:] = quantity_slice*100 # To convert from m to cm
     4126                                           # And m/s to cm/sec
     4127       
     4128    def close(self):
     4129        self.outfile.close()
     4130
     4131def urs2sww(basename_in='o', basename_out=None, verbose=False,
     4132            remove_nc_files=True,
     4133            minlat=None, maxlat=None,
     4134            minlon= None, maxlon=None,
     4135            mint=None, maxt=None,
     4136            mean_stage=0,
     4137            zscale=1,
     4138            fail_on_NaN=True,
     4139            NaN_filler=0,
     4140            elevation=None):
     4141    """
     4142    Convert URS C binary format for wave propagation to
     4143    sww format native to abstract_2d_finite_volumes.
     4144
     4145    Specify only basename_in and read files of the form
     4146    basefilename-z-mux, basefilename-e-mux and basefilename-n-mux containing
     4147    relative height, x-velocity and y-velocity, respectively.
     4148
     4149    Also convert latitude and longitude to UTM. All coordinates are
     4150    assumed to be given in the GDA94 datum. The latitude and longitude
     4151    information is for  a grid.
     4152
     4153    min's and max's: If omitted - full extend is used.
     4154    To include a value min may equal it, while max must exceed it.
     4155    Lat and lon are assuemd to be in decimal degrees
     4156
     4157    URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE
     4158    which means that latitude is the fastest
     4159    varying dimension (row major order, so to speak)
     4160
     4161    In URS C binary the latitudes and longitudes are in assending order.
     4162    """
     4163   
     4164    if basename_out == None:
     4165        basename_out = basename_in
     4166    files_out = urs2nc(basename_in, basename_out)
     4167    ferret2sww(basename_out,
     4168               minlat=minlat,
     4169               maxlat=maxlat,
     4170               minlon=minlon,
     4171               mint=mint,
     4172               maxt=maxt,
     4173               mean_stage=mean_stage,
     4174               zscale=zscale,
     4175               fail_on_NaN=fail_on_NaN,
     4176               NaN_filler=NaN_filler,
     4177               inverted_bathymetry=True,
     4178               verbose=verbose)
     4179    #print "files_out",files_out
     4180    if remove_nc_files:
     4181        for file_out in files_out:
     4182            os.remove(file_out)
     4183   
     4184def urs2nc(basename_in = 'o', basename_out = 'urs'):
     4185    files_in = [basename_in+'-z-mux',
     4186                basename_in+'-e-mux',
     4187                basename_in+'-n-mux']
     4188    files_out = [basename_out+'_ha.nc',
     4189                 basename_out+'_ua.nc',
     4190                 basename_out+'_va.nc']
     4191    quantities = ['HA','UA','VA']
     4192
     4193    hashed_elevation = None
     4194    for file_in, file_out, quantity in map(None, files_in,
     4195                                           files_out,
     4196                                           quantities):
     4197        lonlatdep, lon, lat = _binary_c2nc(file_in,
     4198                                         file_out,
     4199                                         quantity)
     4200        #print "lon",lon
     4201        #print "lat",lat
     4202        if hashed_elevation == None:
     4203            elevation_file = basename_out+'_e.nc'
     4204            write_elevation_sww(elevation_file,
     4205                                lon,
     4206                                lat,
     4207                                lonlatdep[:,2])
     4208            hashed_elevation = myhash(lonlatdep)
     4209        else:
     4210            msg = "The elevation information in the mux files is inconsistent"
     4211            assert hashed_elevation == myhash(lonlatdep), msg
     4212    files_out.append(elevation_file)
     4213    return files_out
     4214   
     4215def _binary_c2nc(file_in, file_out, quantity):
     4216    """
     4217
     4218    return the depth info, so it can be written to a file
     4219    """
     4220    columns = 3 # long, lat , depth
     4221    f = open(file_in, 'rb')
     4222   
     4223    # Number of points/stations
     4224    (points_num,)= unpack('i',f.read(4))
     4225
     4226    # nt, int - Number of time steps
     4227    (time_step_count,)= unpack('i',f.read(4))
     4228
     4229    #dt, float - time step, seconds
     4230    (time_step,) = unpack('f', f.read(4))
     4231   
     4232    msg = "Bad data in the mux file."
     4233    if points_num < 0:
     4234        raise ANUGAError, msg
     4235    if time_step_count < 0:
     4236        raise ANUGAError, msg
     4237    if time_step < 0:
     4238        raise ANUGAError, msg
     4239   
     4240    lonlatdep = p_array.array('f')
     4241    lonlatdep.read(f, columns * points_num)
     4242    lonlatdep = array(lonlatdep, typecode=Float)   
     4243    lonlatdep = reshape(lonlatdep, ( points_num, columns))
     4244   
     4245    lon, lat = lon_lat2grid(lonlatdep)
     4246    nc_file = Write_nc(quantity,
     4247                       file_out,
     4248                        time_step_count,
     4249                        time_step,
     4250                        lon,
     4251                        lat)
     4252
     4253    for i in range(time_step_count):
     4254        #Read in a time slice   
     4255        hz_p_array = p_array.array('f')
     4256        hz_p_array.read(f, points_num)
     4257        hz_p = array(hz_p_array, typecode=Float)
     4258        hz_p = reshape(hz_p, ( len(lat), len(lon)))
     4259       
     4260        nc_file.store_timestep(hz_p)
     4261       
     4262    nc_file.close()
     4263
     4264    return lonlatdep, lon, lat
     4265   
     4266
     4267def write_elevation_sww(file_out, lon, lat, depth_vector):
     4268     
     4269    # NetCDF file definition
     4270    outfile = NetCDFFile(file_out, 'w')
     4271
     4272    #Create new file
     4273    nc_lon_lat_header(outfile, lon, lat)
     4274   
     4275    # ELEVATION
     4276    zname = 'ELEVATION'
     4277    outfile.createVariable(zname, precision, (lat_name, lon_name))
     4278    outfile.variables[zname].units='CENTIMETERS'
     4279    outfile.variables[zname].missing_value=-1.e+034
     4280
     4281    outfile.variables[lon_name][:]= ensure_numeric(lon)
     4282    outfile.variables[lat_name][:]= ensure_numeric(lat)
     4283
     4284    depth = reshape(depth_vector, ( len(lat), len(lon)))
     4285    #print "depth",depth
     4286    outfile.variables[zname][:]= depth
     4287   
     4288    outfile.close()
     4289   
     4290def nc_lon_lat_header(outfile, lon, lat):
     4291   
     4292    outfile.institution = 'Geoscience Australia'
     4293    outfile.description = 'Converted from URS binary C'
     4294   
     4295    # Longitude
     4296    outfile.createDimension(lon_name, len(lon))
     4297    outfile.createVariable(lon_name, precision, (lon_name,))
     4298    outfile.variables[lon_name].point_spacing='uneven'
     4299    outfile.variables[lon_name].units='degrees_east'
     4300    outfile.variables[lon_name].assignValue(lon)
     4301
     4302
     4303    # Latitude
     4304    outfile.createDimension(lat_name, len(lat))
     4305    outfile.createVariable(lat_name, precision, (lat_name,))
     4306    outfile.variables[lat_name].point_spacing='uneven'
     4307    outfile.variables[lat_name].units='degrees_north'
     4308    outfile.variables[lat_name].assignValue(lat)
     4309
     4310
     4311   
     4312def lon_lat2grid(long_lat_dep):
     4313    """
     4314    given a list of points that are assumed to be an a grid,
     4315    return the long's and lat's of the grid.
     4316    long_lat_dep is an array where each row is a position.
     4317    The first column is longitudes.
     4318    The second column is latitudes.
     4319    """
     4320    LONG = 0
     4321    LAT = 1
     4322    points_num = len(long_lat_dep)
     4323    lat = [long_lat_dep[0][LAT]]
     4324    this_rows_long = long_lat_dep[0][LONG]
     4325    i = 1 # Index of long_lat_dep
     4326
     4327    #Working out the lat's
     4328    while long_lat_dep[i][LONG] == this_rows_long and i < points_num:
     4329            lat.append(long_lat_dep[i][LAT])
     4330            i += 1
     4331           
     4332    lats_per_long = i       
     4333    long = [long_lat_dep[0][LONG]]
     4334    #Working out the longs
     4335    while i < points_num:
     4336        msg = 'Data is not gridded.  It must be for this operation'
     4337        assert long_lat_dep[i][LAT] == lat[i%lats_per_long], msg
     4338        if i%lats_per_long == 0:
     4339            long.append(long_lat_dep[i][LONG])
     4340        i += 1
     4341   
     4342    msg = 'Our of range latitudes/longitudes'
     4343    for l in lat:assert -90 < l < 90 , msg
     4344    for l in long:assert -180 < l < 180 , msg
     4345       
     4346    return long, lat
     4347
     4348    ####  END URS 2 SWW  ###     
    40454349#-------------------------------------------------------------
    40464350if __name__ == "__main__":
  • 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.