Changeset 1137 for inundation/ga

Timestamp:

Mar 23, 2005, 6:06:51 PM (20 years ago)

Author:

ole

Message:

Work on sww2asc export, ensure_numeric and other minor stuff

Location:

inundation/ga/storm_surge

Files:

• pyvolution/data_manager.py (modified) (5 diffs)
• pyvolution/interpolate_sww.py (modified) (2 diffs)
• pyvolution/least_squares.py (modified) (7 diffs)
• pyvolution/test_data_manager.py (modified) (4 diffs)
• pyvolution/test_util.py (modified) (1 diff)
• pyvolution/util.py (modified) (12 diffs)
• pyvolution/util_ext.h (modified) (1 diff)
• wiki/issues.txt (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -16 +16 @@
 
 .asc: ASCII foramt of regular DEMs as output from ArcView
+.prj: Associated ArcView file giving more meta data for asc format
+
 .dem: NetCDF representation of regular DEM data
 
 .tsh: ASCII format for storing meshes and associated boundary and region info
+.msh: NetCDF format for storing meshes and associated boundary and region info
 
 .nc: Native ferret NetCDF format
@@ -40 +43 @@
 TSH -> SWW:  Conversion of TSH to sww viewable using Swollen
 
-TSH + Boundary SWW -> SWW: SImluation using pyvolution
+TSH + Boundary SWW -> SWW: Simluation using pyvolution
 
 
@@ -209 +212 @@
 
     def __init__(self, domain, mode = 'w',\
-                 max_size = 2000000000,recursion=False):
+                 max_size = 2000000000,
+                 recursion=False):
         from Scientific.IO.NetCDF import NetCDFFile
         from Numeric import Int, Float, Float32
@@ -978 +982 @@
     outfile.yllcorner = yllcorner #Northing of lower left corner
     outfile.false_easting = false_easting
-    outfile.false_northing =false_northing
+    outfile.false_northing = false_northing
 
     outfile.projection = projection
@@ -1019 +1023 @@
     outfile.close()
 
+
+def sww2asc(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 
+
+    #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
+        
+    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'][:]
+
+    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:
+        xllcorner = fid.xllcorner[0]
+        yllcorner = fid.yllcorner[0]
+        zone = fid.zone[0]
+    else:
+        zone = origin[0]
+        xllcorner = origin[1]
+        yllcorner = origin[2]
+        
+
+    #Get quantity and reduce if applicable
+    q = fid.variables[quantity]
+
+    if len(q.shape) == 2:
+        q_reduced = array( 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
+    
+    #Write prj file
+    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        %f\n' %false_easting)
+    prjid.write('Yshift        %f\n' %false_northing)
+    prjid.write('Parameters\n')
+    prjid.close()
+
+    #Create grid
+    ncols = int((xmax-xmin)/cellsize)+1
+    nrows = int((ymax-ymin)/cellsize)+1
+
+    from Numeric import zeros, Float 
+    points = zeros ( (ncols*nrows, 2), Float )
+
+
+    for i in xrange(ncols):
+        xg = i*cellsize
+        for j in xrange(nrows):
+            yg = j*cellsize
+            k = i*nrows + j
+
+            #print k, xg, yg
+            points[k,0] = xg
+            points[k,1] = yg
+
+    #Interpolate
+    from least_squares import Interpolation
+    
+
+
+    
+    #Close
+    fid.close()
+    
 
 def convert_dem_from_ascii2netcdf(basename_in, basename_out = None,

inundation/ga/storm_surge/pyvolution/interpolate_sww.py

@@ -26 +26 @@
 
 class Interpolate_sww:
-    def __init__(self,file_name, quantity_name):
+    def __init__(self, file_name, quantity_name):
 
         #It's bad to have the quantity_name passed in here.
@@ -62 +62 @@
         self.quantity = quantity 
         
-    def interpolate_xya(self,file_name):
+    def interpolate_xya(self, file_name):
         """
         Given a point file, interpolate the height w.r.t. time at the points

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -99 +99 @@
         old_title_list = mesh_dict['vertex_attribute_titles']
         
-    if verbose:print "tsh file loaded"
+    if verbose: print 'tsh file %s loaded' %mesh_file
    
     # load in the .pts file
@@ -302 +302 @@
           
         """
-
+        from util import ensure_numeric
+        
         #Convert input to Numeric arrays
-        triangles = array(triangles).astype(Int)
-        vertex_coordinates = array(vertex_coordinates).astype(Float)
+        triangles = ensure_numeric(triangles, Int)
+        vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
 
         #Build underlying mesh
         if verbose: print 'Building mesh' 
         #self.mesh = General_mesh(vertex_coordinates, triangles,
-        #FIXME: Trying the normal mesh while testing precrop
+        #FIXME: Trying the normal mesh while testing precrop,
+        #       The functionality of boundary_polygon is needed for that 
         self.mesh = Mesh(vertex_coordinates, triangles,
                          origin = mesh_origin)
@@ -387 +389 @@
 
         from quad import build_quadtree
+        from util import ensure_numeric
 
         if data_origin is None:
@@ -393 +396 @@
         
         #Convert input to Numeric arrays just in case. 
-        point_coordinates = array(point_coordinates).astype(Float)
+        point_coordinates = ensure_numeric(point_coordinates, Float)        
 
 
@@ -595 +598 @@
         """
 
-
+        from util import ensure_numeric
         
         #Convert input to Numeric arrays
-        point_coordinates = array(point_coordinates).astype(Float)
+        point_coordinates = ensure_numeric(point_coordinates, Float)         
         
         #Build n x m interpolation matrix        
@@ -764 +767 @@
           z: Single 1d vector or array of data at the point_coordinates.
         """
-
+        
         #Convert input to Numeric arrays
-        z = array(z).astype(Float)
-
+        from util import ensure_numeric
+        z = ensure_numeric(z, Float)
+        
         if len(z.shape) > 1 :
             raise VectorShapeError, 'Can only deal with 1d data vector'
@@ -808 +812 @@
 
             #Convert input to Numeric arrays
-            z = array(z).astype(Float)
+            from util import ensure_numeric
+            z = ensure_numeric(z, Float)            
             
             #Build n x m interpolation matrix        

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -365 +365 @@
 
 
-
         #Check contents
         #Get NetCDF
@@ -378 +377 @@
         stage = fid.variables['stage']
 
-        #Check values
         #Check values
         Q = self.domain.quantities['stage']
@@ -551 +549 @@
         fid.close()
 
-        #Convert to NetCDF xya
+        #Convert to NetCDF pts
         convert_dem_from_ascii2netcdf(root)
         dem2pts(root)
@@ -582 +580 @@
         os.remove(root + '.asc')
         os.remove(root + '.prj')
+
+
+
+    def test_sww2asc(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' + str(id(self))
+        self.domain.set_datadir('.')
+        self.domain.format = 'sww'
+        self.domain.smooth = True
+        self.domain.set_quantity('elevation', lambda x,y: -x-y)
+        
+        sww = get_dataobject(self.domain)
+        sww.store_connectivity()
+        sww.store_timestep('stage')
+
+        self.domain.evolve_to_end(finaltime = 0.01)
+        sww.store_timestep('stage')
+
+
+        #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']
+        #print x[:]
+        #print y[:]
+        #print z[:]
+        #print stage[:]        
+
+        #Export to ascii/prj files
+        sww2asc(self.domain.filename, 
+                quantity = 'elevation',                         
+                cellsize = 0.25,
+                origin = (56,0,0))
+
+
+        raise 'Under construction (Ole)'
+
+        #Check values
+        #Q = self.domain.quantities['stage']
+        #Q0 = Q.vertex_values[:,0]
+        #Q1 = Q.vertex_values[:,1]
+        #Q2 = Q.vertex_values[:,2]
+        # 
+        #A = stage[1,:]
+        #assert allclose(A[0], min(Q2[0], Q1[1]))
+        #assert allclose(A[1], min(Q0[1], Q1[3], Q2[2]))
+        #assert allclose(A[2], Q0[3])
+        #assert allclose(A[3], min(Q0[0], Q1[5], Q2[4]))
+        #
+        ##Center point
+        #assert allclose(A[4], min(Q1[0], Q2[1], Q0[2],\
+        #                          Q0[5], Q2[6], Q1[7]))
+
+
+        fid.close()
+
+        #Cleanup
+        os.remove(sww.filename)
 
 

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -195 +195 @@
 
 
-
+    def test_ensure_numeric(self):
+        from util import ensure_numeric
+        from Numeric import ArrayType, Float, array
+
+        A = [1,2,3,4]
+        B = ensure_numeric(A)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'l'
+        assert B[0] == 1 and B[1] == 2 and B[2] == 3 and B[3] == 4
+
+
+        A = [1,2,3.14,4]
+        B = ensure_numeric(A)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'd'
+        assert B[0] == 1 and B[1] == 2 and B[2] == 3.14 and B[3] == 4
+
+
+        A = [1,2,3,4]
+        B = ensure_numeric(A, Float)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'd'
+        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
+
+
+        A = [1,2,3,4]
+        B = ensure_numeric(A, Float)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'd'
+        assert B[0] == 1.0 and B[1] == 2.0 and B[2] == 3.0 and B[3] == 4.0
+
+
+        A = array([1,2,3,4])
+        B = ensure_numeric(A)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'l'        
+        assert A == B    
+        assert A is B   #Same object
+
+
+        A = array([1,2,3,4])
+        B = ensure_numeric(A, Float)
+        assert type(B) == ArrayType
+        assert B.typecode() == 'd'        
+        assert A == B    
+        assert A is not B   #Not the same object        
+
+
+        
+        
 
     def test_spatio_temporal_file_function_time(self):

inundation/ga/storm_surge/pyvolution/util.py

@@ -94 +94 @@
       return False  
     
+def ensure_numeric(A, typecode = None):
+    """Ensure that sequence is a Numeric array.
+    Inputs:
+        A: Sequance. If A is already a Numeric array it will be returned
+                     unaltered
+                     If not, an attempt is made to convert it to a Numeric
+                     array
+        typecode: Numeric type. If specified, use this in the conversion.
+                                If not, let Numeric decide
+
+    This function is necessary as array(A) can cause memory overflow.
+    """
+
+    from Numeric import ArrayType, array
+
+    if typecode is None:
+        if type(A) == ArrayType:
+            return A
+        else:
+            return array(A)
+    else:
+        if type(A) == ArrayType:
+            if A.typecode == typecode:
+                return array(A)            
+            else:
+                return A.astype(typecode)
+        else:
+            return array(A).astype(typecode)
+        
+    
 
 def file_function(filename, domain=None, quantities = None, interpolation_points = None, verbose = False):
@@ -253 +283 @@
         
             try:
-                P = array(interpolation_points) 
+                P = ensure_numeric(interpolation_points) 
             except:
                 msg = 'Interpolation points must be an N x 2 Numeric array or a list of points\n'
@@ -498 +528 @@
             values.append(float(value))
 
-        q = array(values)
+        q = ensure_numeric(values)
         
         msg = 'ERROR: File must contain at least one independent value'
@@ -724 +754 @@
             values.append(float(value))
 
-        q = array(values)
+        q = ensure_numeric(values)
         
         msg = 'ERROR: File must contain at least one independent value'
@@ -932 +962 @@
     if verbose: print 'Checking input to inside_polygon'    
     try:
-        points = array(points).astype(Float)       
+        points = ensure_numeric(points, Float)             
     except:
         msg = 'Points could not be converted to Numeric array'
@@ -938 +968 @@
         
     try:
-        polygon = array(polygon).astype(Float)             
+        polygon = ensure_numeric(polygon, Float)           
     except:
         msg = 'Polygon could not be converted to Numeric array'
@@ -967 +997 @@
     if verbose: print 'Checking input to outside_polygon'    
     try:
-        points = array(points).astype(Float)       
+        points = ensure_numeric(points, Float)             
     except:
         msg = 'Points could not be converted to Numeric array'
@@ -973 +1003 @@
         
     try:
-        polygon = array(polygon).astype(Float)             
+        polygon = ensure_numeric(polygon, Float)           
     except:
         msg = 'Polygon could not be converted to Numeric array'
@@ -1041 +1071 @@
     #Input checks
     try:
-        points = array(points).astype(Float)       
+        points = ensure_numeric(points, Float)             
     except:
         msg = 'Points could not be converted to Numeric array'
@@ -1047 +1077 @@
         
     try:
-        polygon = array(polygon).astype(Float)             
+        polygon = ensure_numeric(polygon, Float)           
     except:
         msg = 'Polygon could not be converted to Numeric array'
@@ -1134 +1164 @@
     #Input checks
     try:
-        #FIXME: This is where it can die due to lack of memeory
-        #Maybe comment out and test for Numeric type
-        points = array(points).astype(Float)
+        points = ensure_numeric(points, Float)
     except:
         msg = 'Points could not be converted to Numeric array'
@@ -1143 +1171 @@
     #if verbose: print 'Checking input to separate_points_by_polygon 2' 
     try:
-        #FIXME: SAme thing
-        polygon = array(polygon).astype(Float)
+        polygon = ensure_numeric(polygon, Float)
     except:
         msg = 'Polygon could not be converted to Numeric array'

inundation/ga/storm_surge/pyvolution/util_ext.h

@@ -125 +125 @@
   
   //Convert to consecutive array
-  A = (PyArrayObject*) PyArray_ContiguousFromObject((PyObject*) B, B -> descr -> type, 0, 0);
+  A = (PyArrayObject*) PyArray_ContiguousFromObject((PyObject*) B, 
+                                                    B -> descr -> type, 0, 0);
   
   Py_DECREF(B); //FIXME: Is this really needed??

inundation/ga/storm_surge/wiki/issues.txt

@@ -2 +2 @@
 OPEN ISSUES:
 ------------
+
+
+Issue: Pmesh do not write in msh format
+Importance: Mid-High
+Action: See under pmesh/documentation
+
+Issue: geo_reference is not passed into domain and further on to sww file. Hence we can't export results properly to e.g GIS 
+Importance: High
+Action: Make it happen!
+