Changeset 1833

Timestamp:

Sep 15, 2005, 8:53:36 AM (20 years ago)

Author:

tdhu

Message:

Testing is as yet uncomplete. Test_data_manager is currently running data_manager to import an sww to a ermapper grid.

Location:

inundation/pyvolution

Files:

• data_manager.py (modified) (1 diff)
• ermapper_grids.py (modified) (4 diffs)
• test_data_manager.py (modified) (1 diff)

inundation/pyvolution/data_manager.py

@@ -2709 +2709 @@
     outfile.close()
 
+def sww2ers(basename_in, basename_out = None,
+            quantity = None,
+            timestep = None,
+            reduction = None,
+            cellsize = 10,
+            verbose = False,
+            origin = None,
+            datum = 'WGS84'):
+    
+    """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, reshape, sometrue 
+    import ermapper_grids
+
+    header = {}
+    false_easting = 500000
+    false_northing = 10000000
+    NODATA_value = 0
+
+    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'
+    # Note the use of a .ers extension is optional (write_ermapper_grid will
+    # deal with either option
+    ersfile = basename_out
+##    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
+        yg = (nrows-i)*cellsize # this will flip the order of the y values
+        for j in xrange(ncols):
+            xg = j*cellsize
+            k = i*ncols + j
+
+##            print k, xg, yg
+            grid_points[k,0] = xg
+            grid_points[k,1] = yg
+
+    #Interpolate
+    from least_squares import Interpolation
+    from util import inside_polygon
+
+    #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
+    grid_values = reshape(grid_values,(nrows, ncols))
+
+
+    # setup header information
+    header['datum'] = '"' + datum + '"'
+    # FIXME The use of hardwired UTM and zone number needs to be made optional
+    # FIXME Also need an automatic test for coordinate type (i.e. EN or LL)
+    header['projection'] = '"UTM-' + str(zone) + '"'
+    header['coordinatetype'] = 'EN'
+    if header['coordinatetype'] == 'LL':
+        header['longitude'] = str(newxllcorner)
+        header['latitude'] = str(newyllcorner)
+    elif header['coordinatetype'] == 'EN':
+        header['eastings'] = str(newxllcorner)
+        header['northings'] = str(newyllcorner)
+    header['nullcellvalue'] = str(NODATA_value)
+    header['xdimension'] = str(cellsize)
+    header['ydimension'] = str(cellsize)
+    header['value'] = '"' + quantity + '"'
+
+
+    #Write
+    if verbose: print 'Writing %s' %ersfile
+    ermapper_grids.write_ermapper_grid(ersfile,grid_values, header)
+
+##    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()

inundation/pyvolution/ermapper_grids.py

@@ -5 +5 @@
 
 def write_ermapper_grid(ofile, data, header = {}):
-    # function write_ermapper_grid(ofile, data, header = {}):
-    #
-    # Function to write a 2D Numeric array to an ERMapper grid
-    #
-    # Input Parameters:
-    # ofile
-
+    """
+    write_ermapper_grid(ofile, data, header = {}):
+    
+    Function to write a 2D Numeric array to an ERMapper grid.  There are a series of conventions adopted within
+    this code, specifically:
+    1)  The registration coordinate for the data is the SW (or lower-left) corner of the data
+    2)  The registration coordinates refer to cell centres
+    3)  The data is a 2D Numeric array with the NW-most data in element (0,0) and the SE-most data in element (N,M)
+        where N is the last line and M is the last column
+    4)  There has been no testng of the use of a rotated grid.  Best to keep data in an NS orientation
+    
+    Input Parameters:
+    ofile:      string - filename for output (note the output will consist of two files
+                ofile and ofile.ers.  Either of these can be entered into this function
+    data:       Numeric.array - 2D array containing the data to be output to the grid
+    header:     dictionary - contains spatial information about the grid, in particular:
+                    header['datum'] datum for the data ('"GDA94"')
+                    header['projection'] - either '"GEOGRAPHIC"' or '"PROJECTED"' 
+                    header['coordinatetype'] - either 'EN' (for eastings/northings) or
+                                                      'LL' (for lat/long)
+                    header['rotation'] - rotation of grid ('0:0:0.0')
+                    header['celltype'] - data type for writing data ('IEEE4ByteReal')
+                    header['nullcellvalue'] - value for null cells ('-99999')
+                    header['xdimension'] - cell size in x-dir in units dictated by 'coordinatetype' ('100')
+                    header['registrationcellx'] == '0'
+                    header['ydimension'] - cell size in y-dir in units dictated by 'coordinatetype' ('100')
+                    header['longitude'] - co-ordinate of registration cell ('0:0:0')
+                    header['latitude'] - co-ordinate of registration line ('0:0:0')
+                    header['nrofbands'] - number of bands ('1')
+                    header['value'] - name of grid ('"Default_Band"')
+
+                    Some entries are determined automatically from the data
+                    header['nroflines'] - number of lines in data
+                    header['nrofcellsperline'] - number of columns in data
+                    header['registrationcelly'] == last line of data
+    """
     # extract filenames for header and data files from ofile
     ers_index = ofile.find('.ers')
@@ -79 +108 @@
     fid.write('\t\tCoordinateType\t = ' + header['coordinatetype'] + '\n')
     fid.write('\t\tRotation\t\t = ' + header['rotation'] + '\n')
+    fid.write('\t\tUnits\t\t = ' + header['units'] + '\n')
     fid.write('\tCoordinateSpace End\n')
 
@@ -97 +127 @@
     # Write the registration coordinate information
     fid.write('\t\tRegistrationCoord Begin\n')
+    print X_Class
     fid.write('\t\t\t' + X_Class + '\t\t\t = ' + header[X_Class.lower()] + '\n')
     fid.write('\t\t\t' + Y_Class + '\t\t\t = ' + header[Y_Class.lower()] + '\n')
@@ -163 +194 @@
     if not header.has_key('rotation'):
         header['rotation'] = '0:0:0.0'
+    if not header.has_key('units'):
+        header['units'] = '"METERS"'
     if not header.has_key('celltype'):
         header['celltype'] = 'IEEE4ByteReal'

inundation/pyvolution/test_data_manager.py

@@ -2052 +2052 @@
         os.remove(root + '_100.dem')
 
+    def test_sww2ers_simple(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Setup
+        self.domain.filename = 'datatest'
+
+        headerfile = self.domain.filename + '.ers'
+        swwfile = self.domain.filename + '.sww'
+
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+
+        self.domain.geo_reference = Geo_reference(56,308500,6189000)
+
+        sww = get_dataobject(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep('stage')
+
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep('stage')
+
+        cellsize = 0.25
+        #Check contents
+        #Get NetCDF
+
+        fid = NetCDFFile(sww.filename, 'r')
+
+        # Get the variables
+        x = fid.variables['x'][:]
+        y = fid.variables['y'][:]
+        z = fid.variables['elevation'][:]
+        time = fid.variables['time'][:]
+        stage = fid.variables['stage'][:]
+
+
+        #Export to ascii/prj files
+        sww2ers(self.domain.filename,
+                quantity = 'elevation',
+                cellsize = cellsize,
+                verbose = False)
+
+
+##        #Check prj (meta data)
+##        prjid = open(prjfile)
+##        lines = prjid.readlines()
+##        prjid.close()
+##
+##        L = lines[0].strip().split()
+##        assert L[0].strip().lower() == 'projection'
+##        assert L[1].strip().lower() == 'utm'
+##
+##        L = lines[1].strip().split()
+##        assert L[0].strip().lower() == 'zone'
+##        assert L[1].strip().lower() == '56'
+##
+##        L = lines[2].strip().split()
+##        assert L[0].strip().lower() == 'datum'
+##        assert L[1].strip().lower() == 'wgs84'
+##
+##        L = lines[3].strip().split()
+##        assert L[0].strip().lower() == 'zunits'
+##        assert L[1].strip().lower() == 'no'
+##
+##        L = lines[4].strip().split()
+##        assert L[0].strip().lower() == 'units'
+##        assert L[1].strip().lower() == 'meters'
+##
+##        L = lines[5].strip().split()
+##        assert L[0].strip().lower() == 'spheroid'
+##        assert L[1].strip().lower() == 'wgs84'
+##
+##        L = lines[6].strip().split()
+##        assert L[0].strip().lower() == 'xshift'
+##        assert L[1].strip().lower() == '500000'
+##
+##        L = lines[7].strip().split()
+##        assert L[0].strip().lower() == 'yshift'
+##        assert L[1].strip().lower() == '10000000'
+##
+##        L = lines[8].strip().split()
+##        assert L[0].strip().lower() == 'parameters'
+##
+##
+##        #Check asc file
+##        ascid = open(ascfile)
+##        lines = ascid.readlines()
+##        ascid.close()
+##
+##        L = lines[0].strip().split()
+##        assert L[0].strip().lower() == 'ncols'
+##        assert L[1].strip().lower() == '5'
+##
+##        L = lines[1].strip().split()
+##        assert L[0].strip().lower() == 'nrows'
+##        assert L[1].strip().lower() == '5'
+##
+##        L = lines[2].strip().split()
+##        assert L[0].strip().lower() == 'xllcorner'
+##        assert allclose(float(L[1].strip().lower()), 308500)
+##
+##        L = lines[3].strip().split()
+##        assert L[0].strip().lower() == 'yllcorner'
+##        assert allclose(float(L[1].strip().lower()), 6189000)
+##
+##        L = lines[4].strip().split()
+##        assert L[0].strip().lower() == 'cellsize'
+##        assert allclose(float(L[1].strip().lower()), cellsize)
+##
+##        L = lines[5].strip().split()
+##        assert L[0].strip() == 'NODATA_value'
+##        assert L[1].strip().lower() == '-9999'
+##
+##        #Check grid values
+##        for j in range(5):
+##            L = lines[6+j].strip().split()
+##            y = (4-j) * cellsize
+##            for i in range(5):
+##                assert allclose(float(L[i]), -i*cellsize - y)
+##
+
+        fid.close()
+        print fid
+        #Cleanup
+        #FIXME the file clean-up doesn't work (eg Permission Denied Error)
+        #os.remove(sww.filename)
+
+    def testz_sww2ers_real(self):
+        """Test that sww information can be converted correctly to asc/prj
+        format readable by e.g. ArcView
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+
+
+
+        #Export to ascii/prj files
+        sww2ers('karratha_100m.sww',
+                quantity = 'depth',
+                cellsize = 5,
+                verbose = False)
+
+
+
 
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    suite = unittest.makeSuite(Test_Data_Manager,'test')
+    suite = unittest.makeSuite(Test_Data_Manager,'testz')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')