Changeset 826

Timestamp:

Feb 1, 2005, 5:18:44 PM (20 years ago)

Author:

ole

Message:

Played with DEM's and NetCDF

Location:

inundation/ga/storm_surge/pyvolution

Files:

• data_manager.py (modified) (1 diff)
• least_squares.py (modified) (13 diffs)
• quantity.py (modified) (3 diffs)
• test_data_manager.py (modified) (1 diff)
• util.py (modified) (2 diffs)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -792 +792 @@
 
 
+def dem2xya(filename, verbose=False):
+    """Read Digitial Elevation model 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 ..........
+
+    Convert to NetCDF xyz format (.xya)
+    """ 
+
+    import os
+    from Scientific.IO.NetCDF import NetCDFFile
+    from Numeric import Float, arrayrange, concatenate    
+
+    root, ext = os.path.splitext(filename)
+    fid = open(filename)
+
+    if verbose: print 'Reading DEM from %s' %filename
+    lines = fid.readlines()
+    fid.close()
+
+    if verbose: print 'Got', len(lines), ' lines'
+    
+    ncols = int(lines[0].split()[1].strip())
+    nrows = int(lines[1].split()[1].strip())
+    xllcorner = float(lines[2].split()[1].strip())
+    yllcorner = float(lines[3].split()[1].strip())
+    cellsize = float(lines[4].split()[1].strip())
+    NODATA_value = int(lines[5].split()[1].strip())
+
+    assert len(lines) == nrows + 6
+
+    netcdfname = root + '.xya'
+    if verbose: print 'Store to NetCDF file %s' %netcdfname
+    # NetCDF file definition
+    fid = NetCDFFile(netcdfname, 'w')
+        
+    #Create new file
+    fid.institution = 'Geoscience Australia'
+    fid.description = 'NetCDF xya format for compact and portable storage ' +\
+                      'of spatial point data'
+
+    fid.ncols = ncols
+    fid.nrows = nrows    
+
+
+    # dimension definitions
+    fid.createDimension('number_of_points', nrows*ncols)    
+    fid.createDimension('number_of_attributes', 1) #Always 1 with the dem fmt
+    fid.createDimension('number_of_dimensions', 2) #This is 2d data
+    
+
+    # variable definitions
+    fid.createVariable('points', Float, ('number_of_points',
+                                         'number_of_dimensions'))
+    fid.createVariable('attributes', Float, ('number_of_points',
+                                             'number_of_attributes'))
+
+    # Get handles to the variables
+    points = fid.variables['points']
+    attributes = fid.variables['attributes']
+
+    #Store data
+    #FIXME: Could be faster using array operations
+    #FIXME: I think I swapped x and y here -
+    #       but this function is probably obsolete anyway
+    for i, line in enumerate(lines[6:]):
+        fields = line.split()
+        if verbose: print 'Processing row %d of %d' %(i, nrows)
+        
+        x = i*cellsize            
+        for j, field in enumerate(fields):
+            index = i*ncols + j
+            
+            y = j*cellsize
+            points[index, :] = [x,y]
+
+            z = float(field)                        
+            attributes[index, 0] = z
+
+    fid.close()
+
+                                      
+
+def dem2netcdf(filename, verbose=False):
+    """Read Digitial Elevation model 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 ..........
+
+    Convert to NetCDF format (.dem) mimcing the ASCII format closely.
+    """ 
+
+    import os
+    from Scientific.IO.NetCDF import NetCDFFile
+    from Numeric import Float, array
+
+    root, ext = os.path.splitext(filename)
+    fid = open(filename)
+
+    if verbose: print 'Reading DEM from %s' %filename
+    lines = fid.readlines()
+    fid.close()
+
+    if verbose: print 'Got', len(lines), ' lines'
+    
+    ncols = int(lines[0].split()[1].strip())
+    nrows = int(lines[1].split()[1].strip())
+    xllcorner = float(lines[2].split()[1].strip())
+    yllcorner = float(lines[3].split()[1].strip())
+    cellsize = float(lines[4].split()[1].strip())
+    NODATA_value = int(lines[5].split()[1].strip())
+
+    assert len(lines) == nrows + 6
+
+    netcdfname = root + '.dem'
+    if verbose: print 'Store to NetCDF file %s' %netcdfname
+    # NetCDF file definition
+    fid = NetCDFFile(netcdfname, 'w')
+        
+    #Create new file
+    fid.institution = 'Geoscience Australia'
+    fid.description = 'NetCDF DEM format for compact and portable storage ' +\
+                      'of spatial point data'
+
+    fid.ncols = ncols
+    fid.nrows = nrows
+    fid.xllcorner = xllcorner
+    fid.yllcorner = yllcorner
+    fid.cellsize = cellsize
+    fid.NODATA_value = NODATA_value
+
+    # dimension definitions
+    fid.createDimension('number_of_rows', nrows)
+    fid.createDimension('number_of_columns', ncols)            
+
+    # variable definitions
+    fid.createVariable('elevation', Float, ('number_of_rows',
+                                            'number_of_columns'))
+
+    # Get handles to the variables
+    elevation = fid.variables['elevation']
+
+    #Store data
+    for i, line in enumerate(lines[6:]):
+        fields = line.split()
+        if verbose: print 'Processing row %d of %d' %(i, nrows)
+
+        elevation[i, :] = array([float(x) for x in fields])
+
+    fid.close()
+
+    
+
+
 
 #OBSOLETE STUFF

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -114 +114 @@
                 point_coordinates,
                 point_attributes,
-                alpha = DEFAULT_ALPHA):
+                alpha = DEFAULT_ALPHA,
+                verbose = False):
     """
     Fit a smooth surface to a trianglulation,
@@ -138 +139 @@
                            triangles,
                            point_coordinates,
-                           alpha = alpha)
-    
-    vertex_attributes = interp.fit_points(point_attributes)
+                           alpha = alpha,
+                           verbose = verbose)
+    
+    vertex_attributes = interp.fit_points(point_attributes, verbose = verbose)
     return vertex_attributes
 
     
     
-def xya2rectangular(xya_name, M, N, alpha = DEFAULT_ALPHA):
+def xya2rectangular(xya_name, M, N, alpha = DEFAULT_ALPHA,
+                    verbose = False, reduction = 1):
     """Fits attributes from xya file to MxN rectangular mesh 
     
@@ -156 +159 @@
     
     import util, mesh_factory
-    
-    points, attributes = util.read_xya(xya_name) 
-    
+
+    if verbose: print 'Read xya'
+    points, attributes = util.read_xya(xya_name)
+
+    #Reduce number of points a bit
+    points = points[::reduction]
+    attributes = attributes[::reduction]
+    
+    if verbose: print 'Got %d data points' %len(points)
+
+    if verbose: print 'Create mesh'
     #Find extent
     max_x = min_x = points[0][0]
@@ -179 +190 @@
                         triangles,
                         points,
-                        attributes, alpha=alpha)
+                        attributes, alpha=alpha, verbose=verbose)
 
 
@@ -193 +204 @@
                  triangles,
                  point_coordinates = None,
-                 alpha = DEFAULT_ALPHA):
+                 alpha = DEFAULT_ALPHA,
+                 verbose = False):
 
         
@@ -216 +228 @@
         """
 
-
         #Convert input to Numeric arrays
         vertex_coordinates = array(vertex_coordinates).astype(Float)
@@ -222 +233 @@
         
         #Build underlying mesh
+        if verbose: print 'Building mesh' 
         self.mesh = General_mesh(vertex_coordinates, triangles)
 
@@ -228 +240 @@
 
         #Build coefficient matrices
-        self.build_coefficient_matrix_B(point_coordinates)    
+        self.build_coefficient_matrix_B(point_coordinates, verbose = verbose)
+        
 
     def set_point_coordinates(self, point_coordinates):
@@ -236 +249 @@
         self.build_coefficient_matrix_B(point_coordinates)
         
-    def build_coefficient_matrix_B(self, point_coordinates=None):
+    def build_coefficient_matrix_B(self, point_coordinates=None,
+                                   verbose = False):
         """Build final coefficient matrix"""
         
 
         if self.alpha <> 0:
+            if verbose: print 'Building smoothing matrix'         
             self.build_smoothing_matrix_D()
         
         if point_coordinates:
-
+            if verbose: print 'Building interpolation matrix'         
             self.build_interpolation_matrix_A(point_coordinates)
 
@@ -357 +372 @@
                     #print "PDSG - xi1", xi1
                     #print "PDSG - xi2", xi2
-                   #print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))" \
+                    #print "PDSG element %i verts((%f, %f),(%f, %f),(%f, %f))"\
                     #   % (k, xi0[0], xi0[1], xi1[0], xi1[1], xi2[0], xi2[1])
                     
@@ -366 +381 @@
 
                     
-
                     #Compute interpolation
                     sigma2 = dot((x-xi0), n2)/dot((xi2-xi0), n2)
@@ -603 +617 @@
 
             
-    def fit_points(self, z):
+    def fit_points(self, z, verbose=False):
         """Like fit, but more robust when each point has two or more attributes
         FIXME (Ole): The name fit_points doesn't carry any meaning 
@@ -610 +624 @@
         
         try:
+            if verbose: print 'Solving penalised least_squares problem'
             return self.fit(z)
         except VectorShapeError, e:

inundation/ga/storm_surge/pyvolution/quantity.py

@@ -308 +308 @@
         """
 
-        from Numeric import array, Float, Int
+        from Numeric import array, Float, Int, allclose
 
         values = array(values).astype(Float)
@@ -347 +347 @@
             
         elif location == 'unique vertices':
-            assert len(values.shape) == 1, 'Values array must be 1d'
-            self.set_vertex_values(values, indexes = indexes)
+            assert len(values.shape) == 1 or allclose(values.shape[1:], 1),\
+                   'Values array must be 1d'
+
+            self.set_vertex_values(values.flat, indexes = indexes)
         else:
             if len(values.shape) == 1:
@@ -395 +397 @@
         A = array(A, Float)
 
+        #print 'SHAPE A', A.shape
         assert len(A.shape) == 1
 

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -495 +495 @@
         #Cleanup
         os.remove(sww.filename)
+
+
+
+    def test_dem2xya(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
+
+        #Write test dem file
+        root = 'demtest'+str(time.time())
+
+        filename = root+'.asc'
+        fid = open(filename, 'w')
+        fid.write("""ncols         5
+nrows         6
+xllcorner     2000.5
+yllcorner     3000.5
+cellsize      25
+NODATA_value  -9999
+""")
+        #Create linear function
+
+        ref_points = []
+        ref_attributes = []
+        for i in range(6):
+            x = i*25.0            
+            for j in range(5):
+                y = j*25.0
+                z = x+2*y
+
+                ref_points.append( [x,y] )
+                ref_attributes.append([z])
+                fid.write('%f ' %z)
+            fid.write('\n')    
+
+        fid.close()
+
+        #Convert to NetCDF xya
+        dem2xya(filename)
+
+        #Check contents
+        #Get NetCDF
+        fid = NetCDFFile(root+'.xya.', 'r')
+      
+        # Get the variables
+        points = fid.variables['points']
+        attributes = fid.variables['attributes']
+
+        #Check values
+
+        #print points[:]
+        #print ref_points
+        assert allclose(points, ref_points)
+
+        #print attributes[:]
+        #print ref_attributes
+        assert allclose(attributes, ref_attributes)        
+      
+        #Cleanup
+        fid.close()
+        
+        os.remove(filename)
+        os.remove(root + '.xya')        
+
+        
         
 #-------------------------------------------------------------

inundation/ga/storm_surge/pyvolution/util.py

@@ -1 @@
-
+"""This modul contains various auxiliary function used by pyvolution.
+
+It is also a clearing house for function tat may later earn a module
+of their own.
+"""
 
 def angle(v):
@@ -286 +290 @@
             return res
 
-def read_xya(file_name):
+def read_xya(filename):
     """Read simple xya file, no header, just
     x y [attributes]
     separated by whitespace
+
+    If xya is a NetCDF file it will be read otherwise attemot to read as ASCII
     
     Return list of points and list of attributes
     """
-    
-    fid = open(file_name)
-    lines = fid.readlines()
-    
-    points = [] 
-    attributes = []
-    for line in lines:
-        fields = line.strip().split()
-        
-        points.append( (float(fields[0]),  float(fields[1])) )
-        attributes.append( [float(z) for z in fields[2:] ] )
-        
-    return points, attributes
+
+
+    from Scientific.IO.NetCDF import NetCDFFile
+
+
+    try:
+        #Get NetCDF        
+        fid = NetCDFFile(filename, 'r')
+      
+        # Get the variables
+        points = fid.variables['points']
+        attributes = fid.variables['attributes']
+        ncols = fid.ncols
+        nrows = fid.nrows
+        #fid.close()  #Don't close - arrays are needed outside this function
+    except:
+        #Read as ASCII
+        fid = open(filename)
+        lines = fid.readlines()
+    
+        points = []     
+        attributes = []
+        for line in lines:
+            fields = line.strip().split()
+            points.append( (float(fields[0]),  float(fields[1])) )
+            attributes.append( [float(z) for z in fields[2:] ] )
+        nrows = ncols = None #FIXME: HACK    
+        
+    return points, attributes #, nrows[0], ncols[0]
     
 
 #####################################
-#POLYFON STUFF
+#POLYGON STUFF
 #
 #FIXME: ALl these should be put into new module polygon.py