Changeset 2022

Timestamp:

Nov 14, 2005, 10:15:25 AM (19 years ago)

Author:

duncan

Message:

continuing work on reading gippsland files.

Location:

inundation

Files:

• examples/run_gippsland.py (modified) (1 diff)
• pyvolution/combine_pts.py (modified) (1 diff)
• pyvolution/data_manager.py (modified) (10 diffs)
• pyvolution/test_data_manager.py (modified) (6 diffs)

inundation/examples/run_gippsland.py

@@ -1 +1 @@
 from pyvolution.data_manager import asc_csiro2sww
 
-asc_csiro2sww('bathymetry_expanded','elev_expanded', 'test.sww')
+asc_csiro2sww('bathymetry_expanded','elev_expanded','ucur_expanded',
+              'vcur_expanded', 'test.sww',zscale=1000, mean_stage = 100,
+              fail_on_NaN = False,
+              elevation_NaN_filler = 0)

inundation/pyvolution/combine_pts.py

@@ -148 +148 @@
     #print "extent",extent 
 
-    #get a list of in point indices
-    #FIXME closed doesn't seems to work here.
-    #FIXME (Ole): This functionality has a unittest - so we need a bit 
-    #more information than that :-)
-    #FIXME (dsg): Ha ha. I think changes that we've done have fixed this.
-    # Did the unit tests I did for combine_pts show that it didn't
-    #seem to work? Or that it works now? It's hard to tell.
-    # Let's close this fixme, due to everybogy thinking that there is
-    # no problem.
     inside_indices = inside_polygon(points['pointlist'],
                                     extent, closed=True)

inundation/pyvolution/data_manager.py

@@ -50 +50 @@
 import os
 
-from Numeric import concatenate, array, Float, resize
+from Numeric import concatenate, array, Float, resize, sometrue
 
 from coordinate_transforms.geo_reference import Geo_reference
@@ -1379 +1379 @@
         ptsname = basename_out + '.pts'
 
+    #FIXME (DSG-ON): use loadASCII export_points_file
     if verbose: print 'Store to NetCDF file %s' %ptsname
     # NetCDF file definition
@@ -3266 +3267 @@
     fid.close()
 
-def asc_csiro2sww(bath_dir, elevation_dir, sww_file, verbose=False):
-    """
+def asc_csiro2sww(bath_dir,
+                  elevation_dir,
+                  ucur_dir,
+                  vcur_dir,
+                  sww_file,
+                  zscale=1,
+                  mean_stage = 0,
+                  fail_on_NaN = True,
+                  elevation_NaN_filler = 0,
+                  bath_prefix='ba',
+                  elevation_prefix='el',
+                  verbose=False):
+    """
+
+    Also convert latitude and longitude to UTM. All coordinates are
+    assumed to be given in the GDA94 datum.
     
     assume:
@@ -3278 +3293 @@
     v velocity
 
-    Assume each type is in a seperate directory
-    assume one bath file with extention .000
+    Assumptions
+    Each type is in a seperate directory
+    One bath file with extention .000
+    The time period is less than 24hrs and uniform.
     """
     from Scientific.IO.NetCDF import NetCDFFile
@@ -3289 +3306 @@
     precision = Float # So if we want to change the precision its done here
     
-    # go in to the bath dir and load the 000 file,
+    # go in to the bath dir and load the only file,
     bath_files = os.listdir(bath_dir)
     #print "bath_files",bath_files 
@@ -3303 +3320 @@
     #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_dir_file = elevation_dir  + os.sep + elevation_prefix \
+                         + base_start
     elevation_metadata,elevation_grid =  _read_asc(elevation_dir_file)
     #print "elevation_dir_file",elevation_dir_file 
     #print "elevation_grid", elevation_grid 
+    
+    elevation_files = os.listdir(elevation_dir)
+    ucur_files = os.listdir(ucur_dir)
+    vcur_files = os.listdir(vcur_dir)
+
+    # the first elevation file should be the
+    # file with the same base name as the bath data
+    #print "elevation_files[0]",elevation_files[0]
+    #print "'el' + base_start",'el' + base_start 
+    assert elevation_files[0] == 'el' + base_start
     
     assert bath_metadata == elevation_metadata
@@ -3329 +3356 @@
     # reverse order of lat, so the fist lat represents the first grid row
     latitudes.reverse()
-    
+
+    #Work out the times
+    if len(elevation_files) > 1:
+        # Assume: The time period is less than 24hrs.
+        time_period = (int(elevation_files[1][-3:]) - \
+                      int(elevation_files[0][-3:]))*60*60
+        times = [x*time_period for x in range(len(elevation_files))]
+    else:
+        times = [0.0]
+    #print "times", times
     #print "number_of_latitudes", number_of_latitudes
     #print "number_of_longitudes", number_of_longitudes 
@@ -3335 +3371 @@
     #print "latitudes", latitudes
     #print "longitudes", longitudes
+
     
     ######### WRITE THE SWW FILE #############
@@ -3349 +3386 @@
     outfile.order = 1
     
-    #FIXME(DSG) - hack to get a good sww file
-    number_of_times = 1
     
     # dimension definitions
@@ -3426 +3461 @@
     volumes = array(volumes)
 
-    # FIXME - use coords
-    zone = refzone
-    xllcorner = min(x)
-    yllcorner = min(y)
-
-    outfile.xllcorner = xllcorner
-    outfile.yllcorner = yllcorner
-    outfile.zone = zone
-
+    geo_ref = Geo_reference(refzone,min(x),min(y))
+    geo_ref.write_NetCDF(outfile)
 
     z = resize(bath_grid,outfile.variables['z'][:].shape)
-    outfile.variables['x'][:] = x - xllcorner
-    outfile.variables['y'][:] = y - yllcorner
+    outfile.variables['x'][:] = x - geo_ref.get_xllcorner()
+    outfile.variables['y'][:] = y - geo_ref.get_yllcorner()
     outfile.variables['z'][:] = z
-    #outfile.variables['elevation'][:] = z  #FIXME HACK
-    #outfile.variables['time'][:] = times   #Store time relative
-    #outfile.variables['time'][:] = [0.0]   #Store time relative
+    outfile.variables['elevation'][:] = z  #FIXME HACK
     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
-    
+    # do this to create an ok sww file.
+    #outfile.variables['time'][:] = [0]   #Store time relative
+    #outfile.variables['stage'] = z
+    # put the next line up in the code after outfile.order = 1
+    #number_of_times = 1
+
+    stage = outfile.variables['stage']
+    xmomentum = outfile.variables['xmomentum']
+    ymomentum = outfile.variables['ymomentum']
+    
+    outfile.variables['time'][:] = times   #Store time relative
+    
+    if verbose: print 'Converting quantities'
+    for j in range(number_of_times):
+        # load in files
+        elevation_meta, elevation_grid = \
+           _read_asc(elevation_dir + os.sep + elevation_files[j])
+                         
+        u_meta, u_momentum_grid =  _read_asc(ucur_dir + os.sep + ucur_files[j])
+        v_meta, v_momentum_grid =  _read_asc(vcur_dir + os.sep + vcur_files[j])
+
+        
+        # handle missing values
+        missing = (elevation_grid == elevation_meta['NODATA_value'])
+        if sometrue (missing):
+            if fail_on_NaN:
+                msg = 'File %s contains missing values'\
+                      %(elevation_files[j])
+                raise msg
+            else:
+                elevation_grid = elevation_grid*(missing==0) + \
+                                 missing*elevation_NaN_filler
+
+
+        if verbose and j%((n+10)/10)==0: print '  Doing %d of %d' %(j, n)
+        i = 0
+        for k in range(number_of_latitudes):      #Y direction
+            for l in range(number_of_longitudes): #X direction
+                #w = zscale*amplitudes[j,k,l]/100 + mean_stage
+                w = zscale*elevation_grid[k,l] + mean_stage
+                stage[j,i] = w
+                h = w - z[i]
+                xmomentum[j,i] = u_momentum_grid[k,l]*h
+                ymomentum[j,i] = v_momentum_grid[k,l]*h
+                i += 1
+
     outfile.close()
     

inundation/pyvolution/test_data_manager.py

@@ -2840 +2840 @@
  cellsize       0.25
  nodata_value   -9999.0
-    9000.000    60.000    3000.0
-   -10.000   -20.000   -30.000 
+    9000.000    -1000.000    3000.0
+   -1000.000    9000.000  -1000.000 
 """)
         fid.close()
@@ -2852 +2852 @@
  cellsize       0.25
  nodata_value   -9999.0
-    90.000    60.000    30.0
-     0.000     0.000     0.000 
+    9000.000    0.000    3000.0
+     0.000     9000.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
+    9000.000    4000.000    4000.0
+    4000.000    9000.000    4000.000 
+""")
+        fid.close()
+         
+        ucur_dir =  tempfile.mkdtemp()
+        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
+        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
+        
+        fid = open(ucur_dir_filename1, 'w')
         fid.write(""" ncols             3
  nrows             2
@@ -2868 +2884 @@
 """)
         fid.close()
-        
-        sww_file = 'test.sww'
-        asc_csiro2sww(bath_dir,elevation_dir, sww_file)
+        fid = open(ucur_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()
+        
+        vcur_dir =  tempfile.mkdtemp()
+        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
+        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
+        
+        fid = open(vcur_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
+    10.000    10.000    10.000 
+""")
+        fid.close()
+        fid = open(vcur_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 = 'a_test.sww'
+        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir, sww_file)
 
         # check the sww file
@@ -2877 +2931 @@
         x = fid.variables['x'][:]
         y = fid.variables['y'][:]
+        z = fid.variables['z'][:]
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
         geo_ref = Geo_reference(NetCDFObject=fid)
         #print "geo_ref",geo_ref 
@@ -2898 +2955 @@
         #Easting:  609748.788  Northing: 5793447.860 
         #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 '' 
-        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.860 - y_ref))
-
+        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
+
+        assert allclose(z[0],9000.0 )
+        assert allclose(stage[0][1],0.0 )
+
+        #(4000+1000)*60
+        assert allclose(xmomentum[1][1],300000.0 )
+        
+        
         fid.close()
     
@@ -2914 +2978 @@
         os.remove(sww_file)
         
+    def test_asc_csiro2sww2(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
+    9000.000    -1000.000    3000.0
+   -1000.000    9000.000  -1000.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
+    9000.000    0.000    3000.0
+     0.000     -9999.000     -9999.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
+    9000.000    4000.000    4000.0
+    4000.000    9000.000    4000.000 
+""")
+        fid.close()
+         
+        ucur_dir =  tempfile.mkdtemp()
+        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
+        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
+        
+        fid = open(ucur_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
+    10.000    10.000    10.000 
+""")
+        fid.close()
+        fid = open(ucur_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()
+        
+        vcur_dir =  tempfile.mkdtemp()
+        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
+        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
+        
+        fid = open(vcur_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
+    10.000    10.000    10.000 
+""")
+        fid.close()
+        fid = open(vcur_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()
+        
+        try:
+            asc_csiro2sww(bath_dir,elevation_dir, ucur_dir,
+                          vcur_dir, sww_file)
+        except:
+            #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)
+        else:
+            #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)
+            raise 'Should raise exception'
+        
+     
+    def test_asc_csiro2sww3(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
+    9000.000    -1000.000    3000.0
+   -1000.000    9000.000  -1000.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
+    9000.000    0.000    3000.0
+     0.000     -9999.000     -9999.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
+    9000.000    4000.000    4000.0
+    4000.000    9000.000    4000.000 
+""")
+        fid.close()
+         
+        ucur_dir =  tempfile.mkdtemp()
+        ucur_dir_filename1 = ucur_dir + os.sep +'uc19940524.000'
+        ucur_dir_filename2 = ucur_dir + os.sep +'uc19940524.001'
+        
+        fid = open(ucur_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
+    10.000    10.000    10.000 
+""")
+        fid.close()
+        fid = open(ucur_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()
+        
+        vcur_dir =  tempfile.mkdtemp()
+        vcur_dir_filename1 = vcur_dir + os.sep +'vc19940524.000'
+        vcur_dir_filename2 = vcur_dir + os.sep +'vc19940524.001'
+        
+        fid = open(vcur_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
+    10.000    10.000    10.000 
+""")
+        fid.close()
+        fid = open(vcur_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 = 'a_test.sww'
+        asc_csiro2sww(bath_dir,elevation_dir, ucur_dir, vcur_dir,
+                      sww_file, fail_on_NaN = False, elevation_NaN_filler = 0,
+                      mean_stage = 100)
+
+        # check the sww file
+        
+        fid = NetCDFFile(sww_file, 'r')    #Open existing file for read
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['z'][:]
+        stage = fid.variables['stage'][:]
+        xmomentum = fid.variables['xmomentum'][:]
+        geo_ref = Geo_reference(NetCDFObject=fid)
+        #print "geo_ref",geo_ref 
+        x_ref = geo_ref.get_xllcorner()
+        y_ref = geo_ref.get_yllcorner()
+        self.failUnless(geo_ref.get_zone() == 55,  'Failed')
+        assert allclose(x_ref, 587798.418) # (-38, 148)
+        assert allclose(y_ref, 5793123.477)# (-38, 148.5)
+        
+        #Zone:   55    
+        #Easting:  588095.674  Northing: 5821451.722 
+        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148 0 ' 0.00000 ''
+        assert allclose((x[0],y[0]), (588095.674 - x_ref, 5821451.722 - y_ref))
+
+        #Zone:   55    
+        #Easting:  632145.632  Northing: 5820863.269 
+        #Latitude:   -37  45 ' 0.00000 ''  Longitude: 148  30 ' 0.00000 '' 
+        assert allclose((x[2],y[2]), (632145.632 - x_ref, 5820863.269 - y_ref))
+
+        #Zone:   55    
+        #Easting:  609748.788  Northing: 5793447.860 
+        #Latitude:   -38  0 ' 0.00000 ''  Longitude: 148  15 ' 0.00000 '' 
+        assert allclose((x[4],y[4]), (609748.788  - x_ref, 5793447.86 - y_ref))
+
+        assert allclose(z[0],9000.0 )
+        assert allclose(stage[0][4],100.0 )
+        assert allclose(stage[0][5],100.0 )
+
+        #(100.0 - 9000)*10
+        assert allclose(xmomentum[0][4], -89000.0 )
+        
+        #(100.0 - -1000.000)*10
+        assert allclose(xmomentum[0][5], 11000.0 )
+        
+        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
+        os.remove(sww_file)
+   
+        
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Data_Manager,'test_asc_')
-    #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')