Changeset 2003

Timestamp:

Nov 4, 2005, 5:17:50 PM (18 years ago)

Author:

duncan

Message:

working on gippsland to sww

Location:

inundation/pyvolution

Files:

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

inundation/pyvolution/data_manager.py

@@ -48 +48 @@
 
 """
-
-from Numeric import concatenate, array, Float
+import os
+
+from Numeric import concatenate, array, Float, resize
 
 from coordinate_transforms.geo_reference import Geo_reference
@@ -3056 +3057 @@
     fid.close()
 
-def asc_csiro2sww(bath_file_in, verbose=False):
-    pass
-
+def asc_csiro2sww(bath_dir, elevation_dir, sww_file, verbose=False):
+    """
+    
+    assume:
+    All files are in ersi asc format
+
+    4 types of information
+    bathymetry
+    elevation
+    u velocity
+    v velocity
+
+    Assume each type is in a seperate directory
+    assume one bath file with extention .000
+    """
+    from Scientific.IO.NetCDF import NetCDFFile
+    from Numeric import Float, Int, Int32, searchsorted, zeros, array
+    from Numeric import allclose, around
+        
+    from coordinate_transforms.redfearn import redfearn
+    
+    precision = Float # So if we want to change the precision its done here
+    
+    # go in to the bath dir and load the 000 file,
+    bath_files = os.listdir(bath_dir)
+    #print "bath_files",bath_files 
+
+    #fixme: make more general?
+    bath_file = bath_files[0]
+    bath_dir_file =  bath_dir + os.sep + bath_file
+    bath_metadata,bath_grid =  _read_asc(bath_dir_file)
+    #print "bath_metadata",bath_metadata
+    #print "bath_grid",bath_grid 
+
+    #Use the date.time of the bath file as a basis for
+    #the start time for other files
+    base_start = bath_file[-12:]
+    #print "base_start",base_start
+    
+    #go into the elevation dir and load the 000 file
+    elevation_dir_file = elevation_dir  + os.sep + 'el' + base_start
+    elevation_metadata,elevation_grid =  _read_asc(elevation_dir_file)
+    #print "elevation_dir_file",elevation_dir_file 
+    #print "elevation_grid", elevation_grid 
+    
+    assert bath_metadata == elevation_metadata
+
+
+    
+    
+    number_of_latitudes = bath_grid.shape[0] 
+    number_of_longitudes = bath_grid.shape[1]
+    number_of_times = len(os.listdir(elevation_dir))
+    number_of_points = number_of_latitudes*number_of_longitudes
+    number_of_volumes = (number_of_latitudes-1)*(number_of_longitudes-1)*2
+
+    longitudes = [bath_metadata['xllcorner']+x*bath_metadata['cellsize'] for x in range(number_of_longitudes)]
+
+    
+    latitudes = [bath_metadata['yllcorner']+y*bath_metadata['cellsize'] for y in range(number_of_latitudes)]
+
+    # reverse order of lat, so the fist lat represents the first grid row
+    latitudes.reverse()
+    
+    #print "number_of_latitudes", number_of_latitudes
+    #print "number_of_longitudes", number_of_longitudes 
+    #print "number_of_times", number_of_times 
+    #print "latitudes", latitudes
+    #print "longitudes", longitudes
+    
+    ######### WRITE THE SWW FILE #############
+    # NetCDF file definition
+    outfile = NetCDFFile(sww_file, 'w')
+
+    #Create new file
+    outfile.institution = 'Geoscience Australia'
+    outfile.description = 'Converted from XXX'
+
+
+    #For sww compatibility
+    outfile.smoothing = 'Yes'
+    outfile.order = 1
+    
+
+    # dimension definitions
+    outfile.createDimension('number_of_volumes', number_of_volumes)
+
+    outfile.createDimension('number_of_vertices', 3)
+    outfile.createDimension('number_of_points', number_of_points)
+    outfile.createDimension('number_of_timesteps', number_of_times)
+
+    # variable definitions
+    outfile.createVariable('x', precision, ('number_of_points',))
+    outfile.createVariable('y', precision, ('number_of_points',))
+    outfile.createVariable('elevation', precision, ('number_of_points',))
+
+    #FIXME: Backwards compatibility
+    outfile.createVariable('z', precision, ('number_of_points',))
+    #################################
+
+    outfile.createVariable('volumes', Int, ('number_of_volumes',
+                                            'number_of_vertices'))
+
+    outfile.createVariable('time', precision,
+                           ('number_of_timesteps',))
+
+    outfile.createVariable('stage', precision,
+                           ('number_of_timesteps',
+                            'number_of_points'))
+
+    outfile.createVariable('xmomentum', precision,
+                           ('number_of_timesteps',
+                            'number_of_points'))
+
+    outfile.createVariable('ymomentum', precision,
+                           ('number_of_timesteps',
+                            'number_of_points'))
+
+    #Store
+    from coordinate_transforms.redfearn import redfearn
+    x = zeros(number_of_points, Float)  #Easting
+    y = zeros(number_of_points, Float)  #Northing
+
+    if verbose: print 'Making triangular grid'
+    #Get zone of 1st point.
+    refzone, _, _ = redfearn(latitudes[0],longitudes[0])
+    
+    vertices = {}
+    i = 0
+    for k, lat in enumerate(latitudes):       
+        for l, lon in enumerate(longitudes):  
+
+            vertices[l,k] = i
+
+            zone, easting, northing = redfearn(lat,lon)
+
+            msg = 'Zone boundary crossed at longitude =', lon
+            #assert zone == refzone, msg
+            #print '%7.2f %7.2f %8.2f %8.2f' %(lon, lat, easting, northing)
+            x[i] = easting
+            y[i] = northing
+            i += 1
+
+
+    #Construct 2 triangles per 'rectangular' element
+    volumes = []
+    for l in range(number_of_longitudes-1):    #X direction
+        for k in range(number_of_latitudes-1): #Y direction
+            v1 = vertices[l,k+1]
+            v2 = vertices[l,k]
+            v3 = vertices[l+1,k+1]
+            v4 = vertices[l+1,k]
+
+            volumes.append([v1,v2,v3]) #Upper element
+            volumes.append([v4,v3,v2]) #Lower element
+
+    volumes = array(volumes)
+
+    # FIXME - use coords
+    zone = refzone
+    xllcorner = min(x)
+    yllcorner = min(y)
+
+    outfile.xllcorner = xllcorner
+    outfile.yllcorner = yllcorner
+    outfile.zone = zone
+
+
+    z = resize(bath_grid,outfile.variables['z'][:].shape)
+    outfile.variables['x'][:] = x - xllcorner
+    outfile.variables['y'][:] = y - yllcorner
+    outfile.variables['z'][:] = z
+    #outfile.variables['elevation'][:] = z  #FIXME HACK
+    #outfile.variables['time'][:] = times   #Store time relative
+    #outfile.variables['time'][:] = [0]   #Store time relative
+    outfile.variables['volumes'][:] = volumes.astype(Int32) #On Opteron 64
+
+    # do this to create an ok sww file?
+    #outfile.variables['time'][:] = [0]   #Store time relative
+    #outfile.variables['stage'] = z
+    #outfile.variables['xmomentum'] = z
+    #outfile.variables['ymomentum'] = z
+    
+    outfile.close()
+    
 def _read_asc(filename, verbose=False):
     """Read esri file from the following ASCII format (.asc)
@@ -3071 +3254 @@
     NODATA_value  -9999
     138.3698 137.4194 136.5062 135.5558 ..........
+
     """
 
@@ -3086 +3270 @@
     yllcorner = float(lines.pop(0).split()[1].strip())
     cellsize = float(lines.pop(0).split()[1].strip())
-    NODATA_value = int(lines.pop(0).split()[1].strip())
+    NODATA_value = float(lines.pop(0).split()[1].strip())
     
     assert len(lines) == nrows 
@@ -3100 +3284 @@
     grid = array(grid)
 
-    return xllcorner, yllcorner, cellsize, NODATA_value, grid
-    
+    return {'xllcorner':xllcorner,
+            'yllcorner':yllcorner,
+            'cellsize':cellsize,
+            'NODATA_value':NODATA_value}, grid
+    

inundation/pyvolution/test_data_manager.py

@@ -7 +7 @@
 from util import mean
 import tempfile
+import os
+from Scientific.IO.NetCDF import NetCDFFile
 
 from data_manager import *
@@ -2015 +2017 @@
 
         #Cleanup
-        import os
         os.remove('small.sww')
 
@@ -2678 +2679 @@
 """)
         fid.close()
-        xllcorner, yllcorner, cellsize, NODATA_value, grid = \
+        bath_metadata, grid = \
                    data_manager._read_asc(filename, verbose=False)
-        self.failUnless(xllcorner  == 2000.5,  'Failed')
-        self.failUnless(yllcorner  == 3000.5,  'Failed')
-        self.failUnless(cellsize  == 25,  'Failed')
-        self.failUnless(NODATA_value  == -9999,  'Failed')
+        self.failUnless(bath_metadata['xllcorner']  == 2000.5,  'Failed')
+        self.failUnless(bath_metadata['yllcorner']  == 3000.5,  'Failed')
+        self.failUnless(bath_metadata['cellsize']  == 25,  'Failed')
+        self.failUnless(bath_metadata['NODATA_value']  == -9999,  'Failed')
         self.failUnless(grid[0][0]  == 97.921,  'Failed')
         self.failUnless(grid[3][6]  == 514.980,  'Failed')
+
+        os.remove(filename)
         
-
+    def test_asc_csiro2sww(self):
+        import tempfile
+
+        bath_dir = tempfile.mkdtemp()
+        bath_dir_filename = bath_dir + os.sep +'ba19940524.000' 
+        #bath_dir = 'bath_data_manager_test'
+        #print "os.getcwd( )",os.getcwd( ) 
+        elevation_dir =  tempfile.mkdtemp()
+        #elevation_dir = 'elev_expanded'
+        elevation_dir_filename1 = elevation_dir + os.sep +'el19940524.000' 
+        elevation_dir_filename2 = elevation_dir + os.sep +'el19940524.001' 
+        
+        fid = open(bath_dir_filename, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+   -10.000   -20.000   -30.000 
+""")
+        fid.close()
+        
+        fid = open(elevation_dir_filename1, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+     0.000     0.000     0.000 
+""")
+        fid.close()
+
+        fid = open(elevation_dir_filename2, 'w')
+        fid.write(""" ncols             3
+ nrows             2
+ xllcorner    148.00000
+ yllcorner    -38.00000
+ cellsize       0.25
+ nodata_value   -9999.0
+    90.000    60.000    30.0
+    10.000    10.000    10.000 
+""")
+        fid.close()
+        
+        sww_file = 'test.sww'
+        asc_csiro2sww(bath_dir,elevation_dir, sww_file)
+
+        # check the sww file
+        
+        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        geo_ref = Geo_reference(NetCDFObject=fid)
+        
+        fid.close()
+    
+        #tidy up
+        os.remove(bath_dir_filename)
+        os.rmdir(bath_dir)
+
+        os.remove(elevation_dir_filename1)
+        os.remove(elevation_dir_filename2)
+        os.rmdir(elevation_dir)
+
+
+        # remove sww file
+        #co-ords
+        #Zone:   55    
+        #Easting:  587798.418  Northing: 5793713.242 
+        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  0 ' 0.00000 '' 
+        #Grid Convergence:  0  36 ' 56.52 ''  Point Scale: 0.99969494 
+
+        #Zone:   55    
+        #Easting:  587798.418  Northing: 5793713.242 
+        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  0 ' 0.00000 '' 
+        #Grid Convergence:  0  36 ' 56.52 ''  Point Scale: 0.99969494 
+
+        
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Data_Manager,'test')
+    suite = unittest.makeSuite(Test_Data_Manager,'test_asc_')
+    #suite = unittest.makeSuite(Test_Data_Manager,'test')
     #suite = unittest.makeSuite(Test_Data_Manager,'xxxtest')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_boundingbox')