Changeset 1992

Timestamp:

Nov 2, 2005, 10:20:13 AM (17 years ago)

Author:

duncan

Message:

adding code to convert gippsland to sww

Location:

inundation/pyvolution

Files:

• data_manager.py (modified) (2 diffs)
• test_data_manager.py (modified) (2 diffs)

inundation/pyvolution/data_manager.py

@@ -49 +49 @@
 """
 
-from Numeric import concatenate
+from Numeric import concatenate, array, Float
 
 from coordinate_transforms.geo_reference import Geo_reference
@@ -2839 +2839 @@
     outfile.close()
 
+
+    
+def sww2asc_obsolete(basename_in, basename_out = None,
+            quantity = None,
+            timestep = None,
+            reduction = None,
+            cellsize = 10,
+            verbose = False,
+            origin = None):
+    """Read SWW file and convert to Digitial Elevation model format (.asc)
+
+    Example:
+
+    ncols         3121
+    nrows         1800
+    xllcorner     722000
+    yllcorner     5893000
+    cellsize      25
+    NODATA_value  -9999
+    138.3698 137.4194 136.5062 135.5558 ..........
+
+    Also write accompanying file with same basename_in but extension .prj
+    used to fix the UTM zone, datum, false northings and eastings.
+
+    The prj format is assumed to be as
+
+    Projection    UTM
+    Zone          56
+    Datum         WGS84
+    Zunits        NO
+    Units         METERS
+    Spheroid      WGS84
+    Xshift        0.0000000000
+    Yshift        10000000.0000000000
+    Parameters
+
+
+    if quantity is given, out values from quantity otherwise default to
+    elevation
+
+    if timestep (an index) is given, output quantity at that timestep
+
+    if reduction is given use that to reduce quantity over all timesteps.
+
+    """
+    from Numeric import array, Float, concatenate, NewAxis, zeros,\
+         sometrue
+    from utilities.polygon import inside_polygon
+
+    #FIXME: Should be variable
+    datum = 'WGS84'
+    false_easting = 500000
+    false_northing = 10000000
+
+    if quantity is None:
+        quantity = 'elevation'
+
+    if reduction is None:
+        reduction = max
+
+    if basename_out is None:
+        basename_out = basename_in + '_%s' %quantity
+
+    swwfile = basename_in + '.sww'
+    ascfile = basename_out + '.asc'
+    prjfile = basename_out + '.prj'
+
+
+    if verbose: print 'Reading from %s' %swwfile
+    #Read sww file
+    from Scientific.IO.NetCDF import NetCDFFile
+    fid = NetCDFFile(swwfile)
+
+    #Get extent and reference
+    x = fid.variables['x'][:]
+    y = fid.variables['y'][:]
+    volumes = fid.variables['volumes'][:]
+
+    ymin = min(y); ymax = max(y)
+    xmin = min(x); xmax = max(x)
+
+    number_of_timesteps = fid.dimensions['number_of_timesteps']
+    number_of_points = fid.dimensions['number_of_points']
+    if origin is None:
+
+        #Get geo_reference
+        #sww files don't have to have a geo_ref
+        try:
+            geo_reference = Geo_reference(NetCDFObject=fid)
+        except AttributeError, e:
+            geo_reference = Geo_reference() #Default georef object
+
+        xllcorner = geo_reference.get_xllcorner()
+        yllcorner = geo_reference.get_yllcorner()
+        zone = geo_reference.get_zone()
+    else:
+        zone = origin[0]
+        xllcorner = origin[1]
+        yllcorner = origin[2]
+
+
+    #Get quantity and reduce if applicable
+    if verbose: print 'Reading quantity %s' %quantity
+
+    if quantity.lower() == 'depth':
+        q = fid.variables['stage'][:] - fid.variables['elevation'][:]
+    else:
+        q = fid.variables[quantity][:]
+
+
+    if len(q.shape) == 2:
+        if verbose: print 'Reducing quantity %s' %quantity
+        q_reduced = zeros( number_of_points, Float )
+
+        for k in range(number_of_points):
+            q_reduced[k] = reduction( q[:,k] )
+
+        q = q_reduced
+
+    #Now q has dimension: number_of_points
+
+    #Create grid and update xll/yll corner and x,y
+    if verbose: print 'Creating grid'
+    ncols = int((xmax-xmin)/cellsize)+1
+    nrows = int((ymax-ymin)/cellsize)+1
+
+    newxllcorner = xmin+xllcorner
+    newyllcorner = ymin+yllcorner
+
+    x = x+xllcorner-newxllcorner
+    y = y+yllcorner-newyllcorner
+
+    vertex_points = concatenate ((x[:, NewAxis] ,y[:, NewAxis]), axis = 1)
+    assert len(vertex_points.shape) == 2
+
+
+    from Numeric import zeros, Float
+    grid_points = zeros ( (ncols*nrows, 2), Float )
+
+
+    for i in xrange(nrows):
+        yg = i*cellsize
+        for j in xrange(ncols):
+            xg = j*cellsize
+            k = i*ncols + j
+
+            grid_points[k,0] = xg
+            grid_points[k,1] = yg
+
+    #Interpolate
+    from least_squares import Interpolation
+    
+
+    #FIXME: This should be done with precrop = True, otherwise it'll
+    #take forever. With expand_search  set to False, some grid points might
+    #miss out....
+
+    interp = Interpolation(vertex_points, volumes, grid_points, alpha=0.0,
+                           precrop = False, expand_search = True,
+                           verbose = verbose)
+
+    #Interpolate using quantity values
+    if verbose: print 'Interpolating'
+    grid_values = interp.interpolate(q).flat
+
+    #Write
+    #Write prj file
+    if verbose: print 'Writing %s' %prjfile
+    prjid = open(prjfile, 'w')
+    prjid.write('Projection    %s\n' %'UTM')
+    prjid.write('Zone          %d\n' %zone)
+    prjid.write('Datum         %s\n' %datum)
+    prjid.write('Zunits        NO\n')
+    prjid.write('Units         METERS\n')
+    prjid.write('Spheroid      %s\n' %datum)
+    prjid.write('Xshift        %d\n' %false_easting)
+    prjid.write('Yshift        %d\n' %false_northing)
+    prjid.write('Parameters\n')
+    prjid.close()
+
+
+    
+    if verbose: print 'Writing %s' %ascfile
+    NODATA_value = -9999
+ 
+    ascid = open(ascfile, 'w')
+
+    ascid.write('ncols         %d\n' %ncols)
+    ascid.write('nrows         %d\n' %nrows)
+    ascid.write('xllcorner     %d\n' %newxllcorner)
+    ascid.write('yllcorner     %d\n' %newyllcorner)
+    ascid.write('cellsize      %f\n' %cellsize)
+    ascid.write('NODATA_value  %d\n' %NODATA_value)
+
+
+    #Get bounding polygon from mesh
+    P = interp.mesh.get_boundary_polygon()
+    inside_indices = inside_polygon(grid_points, P)
+
+    for i in range(nrows):
+        if verbose and i%((nrows+10)/10)==0:
+            print 'Doing row %d of %d' %(i, nrows)
+
+        for j in range(ncols):
+            index = (nrows-i-1)*ncols+j
+
+            if sometrue(inside_indices == index):
+                ascid.write('%f ' %grid_values[index])
+            else:
+                ascid.write('%d ' %NODATA_value)
+
+        ascid.write('\n')
+
+    #Close
+    ascid.close()
+    fid.close()
+
+def asc_csiro2sww(bath_file_in, verbose=False):
+    pass
+
+def _read_asc(filename, verbose=False):
+    """Read esri file from the following ASCII format (.asc)
+
+    Example:
+
+    ncols         3121
+    nrows         1800
+    xllcorner     722000
+    yllcorner     5893000
+    cellsize      25
+    NODATA_value  -9999
+    138.3698 137.4194 136.5062 135.5558 ..........
+    """
+
+    datafile = open(filename)
+
+    if verbose: print 'Reading DEM from %s' %(filename)
+    lines = datafile.readlines()
+    datafile.close()
+
+    if verbose: print 'Got', len(lines), ' lines'
+
+    ncols = int(lines.pop(0).split()[1].strip())
+    nrows = int(lines.pop(0).split()[1].strip())
+    xllcorner = float(lines.pop(0).split()[1].strip())
+    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())
+    
+    assert len(lines) == nrows 
+    
+    #Store data
+    grid = []
+    
+    n = len(lines)
+    for i, line in enumerate(lines):
+        cells = line.split()
+        assert len(cells) == ncols
+        grid.append(array([float(x) for x in cells]))
+    grid = array(grid)
+
+    return xllcorner, yllcorner, cellsize, NODATA_value, grid
+    

inundation/pyvolution/test_data_manager.py

@@ -6 +6 @@
 from Numeric import zeros, array, allclose, Float
 from util import mean
+import tempfile
 
 from data_manager import *
@@ -2653 +2654 @@
                 verbose = True)
 
-
-
+    def test_read_asc(self):
+        """Test conversion from dem in ascii format to native NetCDF xya format
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        import data_manager
+        #Write test asc file
+        filename = tempfile.mktemp(".000")
+        fid = open(filename, 'w')
+        fid.write("""ncols         7
+nrows         4
+xllcorner     2000.5
+yllcorner     3000.5
+cellsize      25
+NODATA_value  -9999
+    97.921    99.285   125.588   180.830   258.645   342.872   415.836
+   473.157   514.391   553.893   607.120   678.125   777.283   883.038
+   984.494  1040.349  1008.161   900.738   730.882   581.430   514.980
+   502.645   516.230   504.739   450.604   388.500   338.097   514.980
+""")
+        fid.close()
+        xllcorner, yllcorner, cellsize, NODATA_value, 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(grid[0][0]  == 97.921,  'Failed')
+        self.failUnless(grid[3][6]  == 514.980,  'Failed')
+        
 
 #-------------------------------------------------------------