Changeset 3720

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:

• data_manager.py (modified) (3 diffs)
• test_data_manager.py (modified) (4 diffs)

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -60 +60 @@
 import csv
 import os
+from struct import unpack
+import array as p_array
 
 from Numeric import concatenate, array, Float, Int, Int32, resize, sometrue, \
-     searchsorted, zeros, allclose, around
+     searchsorted, zeros, allclose, around, reshape
+from Scientific.IO.NetCDF import NetCDFFile
 
 
@@ -68 +71 @@
 from anuga.geospatial_data.geospatial_data import Geospatial_data
 from anuga.config import minimum_storable_height as default_minimum_storable_height
+from anuga.utilities.numerical_tools import ensure_numeric
+from anuga.caching.caching import myhash
+from anuga.utilities.anuga_exceptions import ANUGAError
 
 def make_filename(s):
@@ -4043 +4049 @@
     return
 
+    ####  URS 2 SWW  ###
+
+lon_name = 'LON'
+lat_name = 'LAT'
+time_name = 'TIME'
+precision = Float # So if we want to change the precision its done here        
+class Write_nc:
+    """
+    Write an nc file.
+    
+    Note, this should be checked to meet cdc netcdf conventions for gridded
+    data. http://www.cdc.noaa.gov/cdc/conventions/cdc_netcdf_standard.shtml
+    
+    """
+    def __init__(self,
+                 quantity_name,
+                 file_name,
+                 time_step_count,
+                 time_step,
+                 lon,
+                 lat):
+        """
+        time_step_count is the number of time steps.
+        time_step is the time step size
+        
+        pre-condition: quantity_name must be 'HA' 'UA'or 'VA'.
+        """
+        self.quantity_name = quantity_name
+        quantity_units= {'HA':'METERS',
+                              'UA':'METERS/SECOND',
+                              'VA':'METERS/SECOND'}
+        
+        #self.file_name = file_name
+        self.time_step_count = time_step_count
+        self.time_step = time_step
+
+        # NetCDF file definition
+        self.outfile = NetCDFFile(file_name, 'w')
+        outfile = self.outfile       
+
+        #Create new file
+        nc_lon_lat_header(outfile, lon, lat)
+    
+        # TIME
+        outfile.createDimension(time_name, None)
+        outfile.createVariable(time_name, precision, (time_name,))
+
+        #QUANTITY
+        outfile.createVariable(self.quantity_name, precision,
+                               (time_name, lat_name, lon_name))
+        outfile.variables[self.quantity_name].missing_value=-1.e+034
+        outfile.variables[self.quantity_name].units= \
+                                 quantity_units[self.quantity_name]
+        outfile.variables[lon_name][:]= ensure_numeric(lon)
+        outfile.variables[lat_name][:]= ensure_numeric(lat)
+
+        #Assume no one will be wanting to read this, while we are writing
+        #outfile.close()
+        
+    def store_timestep(self, quantity_slice):
+        """
+        quantity_slice is the data to be stored at this time step
+        """
+        
+        outfile = self.outfile
+        
+        # Get the variables
+        time = outfile.variables[time_name]
+        quantity = outfile.variables[self.quantity_name]
+            
+        i = len(time)
+
+        #Store time
+        time[i] = i*self.time_step #self.domain.time
+        quantity[i,:] = quantity_slice*100 # To convert from m to cm
+                                           # And m/s to cm/sec
+        
+    def close(self):
+        self.outfile.close()
+
+def urs2sww(basename_in='o', basename_out=None, verbose=False,
+            remove_nc_files=True,
+            minlat=None, maxlat=None,
+            minlon= None, maxlon=None,
+            mint=None, maxt=None,
+            mean_stage=0,
+            zscale=1,
+            fail_on_NaN=True,
+            NaN_filler=0,
+            elevation=None):
+    """
+    Convert URS C binary format for wave propagation to
+    sww format native to abstract_2d_finite_volumes.
+
+    Specify only basename_in and read files of the form
+    basefilename-z-mux, basefilename-e-mux and basefilename-n-mux containing
+    relative height, x-velocity and y-velocity, respectively.
+
+    Also convert latitude and longitude to UTM. All coordinates are
+    assumed to be given in the GDA94 datum. The latitude and longitude
+    information is for  a grid.
+
+    min's and max's: If omitted - full extend is used.
+    To include a value min may equal it, while max must exceed it.
+    Lat and lon are assuemd to be in decimal degrees
+
+    URS C binary format has data orgainised as TIME, LONGITUDE, LATITUDE
+    which means that latitude is the fastest
+    varying dimension (row major order, so to speak)
+
+    In URS C binary the latitudes and longitudes are in assending order.
+    """
+    
+    if basename_out == None:
+        basename_out = basename_in
+    files_out = urs2nc(basename_in, basename_out)
+    ferret2sww(basename_out,
+               minlat=minlat,
+               maxlat=maxlat,
+               minlon=minlon,
+               mint=mint,
+               maxt=maxt,
+               mean_stage=mean_stage,
+               zscale=zscale,
+               fail_on_NaN=fail_on_NaN,
+               NaN_filler=NaN_filler,
+               inverted_bathymetry=True,
+               verbose=verbose)
+    #print "files_out",files_out
+    if remove_nc_files:
+        for file_out in files_out:
+            os.remove(file_out)
+    
+def urs2nc(basename_in = 'o', basename_out = 'urs'):
+    files_in = [basename_in+'-z-mux',
+                basename_in+'-e-mux',
+                basename_in+'-n-mux']
+    files_out = [basename_out+'_ha.nc',
+                 basename_out+'_ua.nc',
+                 basename_out+'_va.nc']
+    quantities = ['HA','UA','VA']
+
+    hashed_elevation = None
+    for file_in, file_out, quantity in map(None, files_in,
+                                           files_out,
+                                           quantities):
+        lonlatdep, lon, lat = _binary_c2nc(file_in,
+                                         file_out,
+                                         quantity)
+        #print "lon",lon
+        #print "lat",lat 
+        if hashed_elevation == None:
+            elevation_file = basename_out+'_e.nc'
+            write_elevation_sww(elevation_file,
+                                lon,
+                                lat,
+                                lonlatdep[:,2])
+            hashed_elevation = myhash(lonlatdep)
+        else:
+            msg = "The elevation information in the mux files is inconsistent"
+            assert hashed_elevation == myhash(lonlatdep), msg
+    files_out.append(elevation_file)
+    return files_out
+    
+def _binary_c2nc(file_in, file_out, quantity):
+    """
+
+    return the depth info, so it can be written to a file
+    """
+    columns = 3 # long, lat , depth
+    f = open(file_in, 'rb')
+    
+    # Number of points/stations
+    (points_num,)= unpack('i',f.read(4))
+
+    # nt, int - Number of time steps
+    (time_step_count,)= unpack('i',f.read(4))
+
+    #dt, float - time step, seconds
+    (time_step,) = unpack('f', f.read(4))
+    
+    msg = "Bad data in the mux file."
+    if points_num < 0:
+        raise ANUGAError, msg
+    if time_step_count < 0:
+        raise ANUGAError, msg
+    if time_step < 0:
+        raise ANUGAError, msg
+    
+    lonlatdep = p_array.array('f')
+    lonlatdep.read(f, columns * points_num)
+    lonlatdep = array(lonlatdep, typecode=Float)    
+    lonlatdep = reshape(lonlatdep, ( points_num, columns))
+    
+    lon, lat = lon_lat2grid(lonlatdep)
+    nc_file = Write_nc(quantity,
+                       file_out,
+                        time_step_count,
+                        time_step,
+                        lon,
+                        lat)
+
+    for i in range(time_step_count):
+        #Read in a time slice    
+        hz_p_array = p_array.array('f')
+        hz_p_array.read(f, points_num)
+        hz_p = array(hz_p_array, typecode=Float)
+        hz_p = reshape(hz_p, ( len(lat), len(lon)))
+        
+        nc_file.store_timestep(hz_p)
+        
+    nc_file.close()
+
+    return lonlatdep, lon, lat
+    
+
+def write_elevation_sww(file_out, lon, lat, depth_vector):
+      
+    # NetCDF file definition
+    outfile = NetCDFFile(file_out, 'w')
+
+    #Create new file
+    nc_lon_lat_header(outfile, lon, lat)
+    
+    # ELEVATION
+    zname = 'ELEVATION'
+    outfile.createVariable(zname, precision, (lat_name, lon_name))
+    outfile.variables[zname].units='CENTIMETERS'
+    outfile.variables[zname].missing_value=-1.e+034
+
+    outfile.variables[lon_name][:]= ensure_numeric(lon)
+    outfile.variables[lat_name][:]= ensure_numeric(lat)
+
+    depth = reshape(depth_vector, ( len(lat), len(lon)))
+    #print "depth",depth 
+    outfile.variables[zname][:]= depth
+    
+    outfile.close()
+    
+def nc_lon_lat_header(outfile, lon, lat):
+    
+    outfile.institution = 'Geoscience Australia'
+    outfile.description = 'Converted from URS binary C'
+    
+    # Longitude
+    outfile.createDimension(lon_name, len(lon))
+    outfile.createVariable(lon_name, precision, (lon_name,))
+    outfile.variables[lon_name].point_spacing='uneven'
+    outfile.variables[lon_name].units='degrees_east'
+    outfile.variables[lon_name].assignValue(lon)
+
+
+    # Latitude
+    outfile.createDimension(lat_name, len(lat))
+    outfile.createVariable(lat_name, precision, (lat_name,))
+    outfile.variables[lat_name].point_spacing='uneven'
+    outfile.variables[lat_name].units='degrees_north'
+    outfile.variables[lat_name].assignValue(lat)
+
+
+    
+def lon_lat2grid(long_lat_dep):
+    """
+    given a list of points that are assumed to be an a grid,
+    return the long's and lat's of the grid.
+    long_lat_dep is an array where each row is a position.
+    The first column is longitudes.
+    The second column is latitudes.
+    """
+    LONG = 0
+    LAT = 1
+    points_num = len(long_lat_dep)
+    lat = [long_lat_dep[0][LAT]]
+    this_rows_long = long_lat_dep[0][LONG]
+    i = 1 # Index of long_lat_dep 
+
+    #Working out the lat's
+    while long_lat_dep[i][LONG] == this_rows_long and i < points_num:
+            lat.append(long_lat_dep[i][LAT])
+            i += 1
+            
+    lats_per_long = i        
+    long = [long_lat_dep[0][LONG]]
+    #Working out the longs
+    while i < points_num:
+        msg = 'Data is not gridded.  It must be for this operation'
+        assert long_lat_dep[i][LAT] == lat[i%lats_per_long], msg
+        if i%lats_per_long == 0:
+            long.append(long_lat_dep[i][LONG])
+        i += 1
+    
+    msg = 'Our of range latitudes/longitudes'
+    for l in lat:assert -90 < l < 90 , msg
+    for l in long:assert -180 < l < 180 , msg
+        
+    return long, lat
+
+    ####  END URS 2 SWW  ###     
 #-------------------------------------------------------------
 if __name__ == "__main__":

anuga_core/source/anuga/shallow_water/test_data_manager.py

@@ -9 +9 @@
 import os
 from Scientific.IO.NetCDF import NetCDFFile
+from struct import pack
 
 from anuga.shallow_water.data_manager import *
@@ -16 +17 @@
 # This is needed to run the tests of local functions
 import data_manager 
-
+from anuga.coordinate_transforms.redfearn import redfearn
 from anuga.coordinate_transforms.geo_reference import Geo_reference
 
@@ -71 +72 @@
         self.F = bed
 
-
-
-
         #Write A testfile (not realistic. Values aren't realistic)
-
         self.test_MOST_file = 'most_small'
 
@@ -4984 +4981 @@
 
         os.remove(csv_file)
+
+    #### TESTS FOR URS 2 SWW  ###     
+    
+    def create_mux(self, points_num=None):
+        # write all the mux stuff.
+        time_step_count = 3
+        time_step = 0.5
+        
+        longitudes = [150.66667, 150.83334, 151., 151.16667]
+        latitudes = [-34.5, -34.33333, -34.16667, -34]
+
+        if points_num == None:
+            points_num = len(longitudes) * len(latitudes)
+
+        lonlatdeps = []
+        quantities = ['HA','UA','VA']
+        mux_names = ['-z-mux','-e-mux','-n-mux']
+        quantities_init = [[],[],[]]
+        for i,lon in enumerate(longitudes):
+            for j,lat in enumerate(latitudes):
+                _ , e, n = redfearn(lat, lon)
+                lonlatdeps.append([lon, lat, n])
+                quantities_init[0].append(e) # HA
+                quantities_init[1].append(n ) # UA
+                quantities_init[2].append(e) # VA
+        #print "lonlatdeps",lonlatdeps 
+        base_name = str(id(self))
+        files = []        
+        for i,q in enumerate(quantities): 
+            quantities_init[i] = ensure_numeric(quantities_init[i])
+            #print "HA_init", HA_init
+            q_time = zeros((time_step_count, points_num), Float)
+            for time in range(time_step_count):
+                q_time[time,:] = quantities_init[i] #* time * 4
+            
+            #Write C files
+            columns = 3 # long, lat , depth
+            file = base_name + mux_names[i]
+            f = open(file, 'wb')
+            files.append(file)
+            f.write(pack('i',points_num))
+            f.write(pack('i',time_step_count))
+            f.write(pack('f',time_step))
+
+            #write lat/long info
+            for lonlatdep in lonlatdeps:
+                for float in lonlatdep:
+                    f.write(pack('f',float))
+                    
+            # Write quantity info
+            for time in  range(time_step_count):
+                for i in range(points_num):
+                    f.write(pack('f',q_time[time,i]))
+            f.close()
+        return base_name, files
+        
+    
+    def delete_mux(self, files):
+        for file in files:
+            os.remove(file)
+            
+    def test_urs2sww_test_fail(self):
+        points_num = -100
+        time_step_count = 45
+        time_step = -7
+        base_name = str(id(self))
+        files = []
+        quantities = ['HA','UA','VA']
+        mux_names = ['-z-mux','-e-mux','-n-mux']
+        for i,q in enumerate(quantities): 
+            #Write C files
+            columns = 3 # long, lat , depth
+            file = base_name + mux_names[i]
+            f = open(file, 'wb')
+            files.append(file)
+            f.write(pack('i',points_num))
+            f.write(pack('i',time_step_count))
+            f.write(pack('f',time_step))
+
+            f.close()   
+        tide = 1
+        try:
+            urs2sww(base_name, remove_nc_files=True, mean_stage=tide)        
+        except ANUGAError:
+            pass
+        else:
+            self.delete_mux(files)
+            msg = 'Should have raised exception'
+            raise msg
+        sww_file = base_name + '.sww'
+        
+        self.delete_mux(files)
+        
+    def test_urs2sww(self):
+        tide = 1
+        base_name, files = self.create_mux()
+        urs2sww(base_name, remove_nc_files=True, mean_stage=tide)
+        sww_file = base_name + '.sww'
+        
+        #Let's interigate the sww file 
+        fid = NetCDFFile(sww_file)
+
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_reference = Geo_reference(NetCDFObject=fid)
+        
+        #Check that first coordinate is correctly represented       
+        #Work out the UTM coordinates for first point
+        zone, e, n = redfearn(-34.5, 150.66667)       
+       
+        assert allclose(geo_reference.get_absolute([[x[0],y[0]]]), [e,n])
+
+        #Check first value
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        ymomentum = fid.variables['ymomentum'][:]
+
+        #print ymomentum
+        assert allclose(stage[0,0], e +tide)  #Meters
+
+        #Check the momentums - ua
+        #momentum = velocity*(stage-elevation)
+        #momentum = velocity*(stage+elevation)
+        # -(-elevation) since elevation is inverted in mux files
+        # = n*(e+tide+n) based on how I'm writing these files
+        answer = n*(e+tide+n)
+        actual = xmomentum[0,0]
+        #print "answer",answer
+        #print "actual",actual 
+        assert allclose(answer, actual)  #Meters
+        
+        fid.close()
+
+
+        self.delete_mux(files)
+        os.remove(sww_file)
+        
+    def bad_test_assuming_theres_a_mux_file(self):
+        # these mux files aren't in the repository, plus they are bad!
+        base_name = 'o-z-mux' 
+        base_name = 'o-e-mux'
+        file_name = base_name + '.nc'
+        lonlatdep_numeric, lon, lat = _binary_c2nc(base_name, file_name, 'HA')
+        
+        #os.remove(file_name)
+        
+    def test_lon_lat2grid(self):
+        lonlatdep = [
+            [ 113.06700134  ,  -26.06669998 ,   0.        ] ,
+            [ 113.06700134  ,  -26.33329964 ,   0.        ] ,
+            [ 113.19999695  ,  -26.06669998 ,   0.        ] ,
+            [ 113.19999695  ,  -26.33329964 ,   0.        ] ]
+            
+        long, lat = lon_lat2grid(lonlatdep)
+
+        for i, result in enumerate(lat):
+            assert lonlatdep [i][1] == result
+        assert len(lat) == 2 
+
+        for i, result in enumerate(long):
+            assert lonlatdep [i*2][0] == result
+        assert len(long) == 2
+
+    def test_lon_lat2grid_bad(self):
+        lonlatdep  = [
+            [ -26.06669998,  113.06700134,    1.        ],
+            [ -26.06669998 , 113.19999695 ,   2.        ],
+            [ -26.06669998 , 113.33300018,    3.        ],
+            [ -26.06669998 , 113.43299866   , 4.        ],
+            [ -26.20000076 , 113.06700134,    5.        ],
+            [ -26.20000076 , 113.19999695 ,   6.        ],
+            [ -26.20000076 , 113.33300018  ,  7.        ],
+            [ -26.20000076 , 113.43299866   , 8.        ],
+            [ -26.33329964 , 113.06700134,    9.        ],
+            [ -26.33329964 , 113.19999695 ,   10.        ],
+            [ -26.33329964 , 113.33300018  ,  11.        ],
+            [ -26.33329964 , 113.43299866 ,   12.        ],
+            [ -26.43330002 , 113.06700134 ,   13        ],
+            [ -26.43330002 , 113.19999695 ,   14.        ],
+            [ -26.43330002 , 113.33300018,    15.        ],
+            [ -26.43330002 , 113.43299866,    16.        ]]
+        try:
+            long, lat = lon_lat2grid(lonlatdep)
+        except AssertionError:
+            pass
+        else:
+            msg = 'Should have raised exception'
+            raise msg
+       
+    def test_lon_lat2gridII(self):
+        lonlatdep = [
+            [ 113.06700134  ,  -26.06669998 ,   0.        ] ,
+            [ 113.06700134  ,  -26.33329964 ,   0.        ] ,
+            [ 113.19999695  ,  -26.06669998 ,   0.        ] ,
+            [ 113.19999695  ,  -26.344329964 ,   0.        ] ]
+        try:
+            long, lat = lon_lat2grid(lonlatdep)
+        except AssertionError:
+            pass
+        else:
+            msg = 'Should have raised exception'
+            raise msg
+        
+    def trial_loading(self):
+        basename_in = 'karratha'
+        basename_out = basename_in
+        urs2sww(basename_in, basename_out)
+    #### END TESTS FOR URS 2 SWW  ###
+
         
 #-------------------------------------------------------------