Changeset 834

Timestamp:

Feb 4, 2005, 1:58:09 PM (20 years ago)

Author:

ole

Message:

Fixed up dem conversions

Location:

inundation/ga/storm_surge/pyvolution

Files:

• data_manager.py (modified) (3 diffs)
• least_squares.py (modified) (2 diffs)
• test_data_manager.py (modified) (2 diffs)
• test_least_squares.py (modified) (1 diff)
• test_util.py (modified) (2 diffs)
• util.py (modified) (3 diffs)

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -6 +6 @@
 
 All data is assumed to reside at vertex locations.  
+
+
+Formats used within GANUGA:
+
+.sww: Netcdf format for storing model output
+.xya: Either ASCII or NetCDF format for storing arbitrary points and associated attributed
+.dem: NetCDF representation of regular DEM data
+.tsh: ASCII format for storing meshes and associated boundary and region info
+FIXME: What else
 
 """
@@ -790 +799 @@
     return cls(domain, mode)
 
-
-
 def dem2xya(filename, verbose=False):
-    """Read Digitial Elevation model from the following NetCDF format (.asc)
+    """Read Digitial Elevation model from the following NetCDF format (.dem)
 
     Example:
@@ -813 +820 @@
 
     root, ext = os.path.splitext(filename)
-    fid = open(filename)
+
+    #Get NetCDF
+    infile = NetCDFFile(filename, 'r')  #Open existing netcdf file for read
 
     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
+    
+    ncols = infile.ncols[0]
+    nrows = infile.nrows[0]
+    xllcorner = infile.xllcorner[0]
+    yllcorner = infile.yllcorner[0]
+    cellsize = infile.cellsize[0]
+    NODATA_value = infile.NODATA_value[0]
+    elevation = infile.variables['elevation']
+
+    #Get output file                                
+    xyaname = root + '.xya'
+    if verbose: print 'Store to NetCDF file %s' %xyaname
     # NetCDF file definition
-    fid = NetCDFFile(netcdfname, 'w')
+    outfile = NetCDFFile(xyaname, 'w')
         
     #Create new file
-    fid.institution = 'Geoscience Australia'
-    fid.description = 'NetCDF xya format for compact and portable storage ' +\
+    outfile.institution = 'Geoscience Australia'
+    outfile.description = 'NetCDF xya format for compact and portable storage ' +\
                       'of spatial point data'
 
-    fid.ncols = ncols
-    fid.nrows = nrows    
+    outfile.ncols = ncols
+    outfile.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
+    outfile.createDimension('number_of_points', nrows*ncols)    
+    outfile.createDimension('number_of_attributes', 1) #Always 1 with the dem fmt
+    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
     
 
     # variable definitions
-    fid.createVariable('points', Float, ('number_of_points',
+    outfile.createVariable('points', Float, ('number_of_points',
                                          'number_of_dimensions'))
-    fid.createVariable('attributes', Float, ('number_of_points',
+    outfile.createVariable('attributes', Float, ('number_of_points',
                                              'number_of_attributes'))
 
     # Get handles to the variables
-    points = fid.variables['points']
-    attributes = fid.variables['attributes']
+    points = outfile.variables['points']
+    attributes = outfile.variables['attributes']
 
     #Store data
     #FIXME: Could be faster using array operations
-    #FIXME: Perhaps the y dimension needs to be reversed
-    #FIXME: x and y definitely needs to be swapped
-    for i, line in enumerate(lines[6:]):
-        fields = line.split()
+    for i in range(nrows):
         if verbose: print 'Processing row %d of %d' %(i, nrows)
         
-        x = i*cellsize            
-        for j, field in enumerate(fields):
+        y = (nrows-i)*cellsize            
+        for j in range(ncols):
             index = i*ncols + j
             
-            y = j*cellsize
+            x = j*cellsize
             points[index, :] = [x,y]
-
-            z = float(field)                        
-            attributes[index, 0] = z
-
-    fid.close()
+            attributes[index, 0] = elevation[i, j]            
+
+    infile.close()            
+    outfile.close()
+
+
+# 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: Perhaps the y dimension needs to be reversed
+#     #FIXME: x and y definitely needs to be swapped
+#     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 asciidem2netcdf(filename, verbose=False):
+def convert_dem_from_ascii2netcdf(filename, verbose=False):
     """Read Digitial Elevation model from the following ASCII format (.asc)
 

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -148 +148 @@
     
 def xya2rectangular(xya_name, M, N, alpha = DEFAULT_ALPHA,
-                    verbose = False, reduction = 1):
+                    verbose = False, reduction = 1, format = 'netcdf'):
     """Fits attributes from xya file to MxN rectangular mesh 
     
@@ -161 +161 @@
 
     if verbose: print 'Read xya'
-    points, attributes, _ = util.read_xya(xya_name)
+    points, attributes, _ = util.read_xya(xya_name, format)
 
     #Reduce number of points a bit

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -523 +523 @@
         ref_attributes = []
         for i in range(6):
-            x = i*25.0            
+            y = (6-i)*25.0            
             for j in range(5):
-                y = j*25.0
+                x = j*25.0
                 z = x+2*y
 
@@ -536 +536 @@
 
         #Convert to NetCDF xya
-        asciidem2netcdf(filename)
-        dem2xya(filename)
+        convert_dem_from_ascii2netcdf(filename)
+        dem2xya(root + '.dem')
 
         #Check contents

inundation/ga/storm_surge/pyvolution/test_least_squares.py

@@ -499 +499 @@
         fid.close()
         
-        points, triangles, boundary, attributes = xya2rectangular(FN, 4, 4)
+        points, triangles, boundary, attributes =\
+                xya2rectangular(FN, 4, 4, format = 'asc')
         
-        
-        
-
         z1 = [2, 4]
         z2 = [4, 8]

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -346 +346 @@
         fid.close()
         
-        points, attributes, _ = read_xya(FN)
+        points, attributes, _ = read_xya(FN, format = 'asc')
         
         assert allclose(points, [ [0,1], [1,0], [1,1] ])
@@ -364 +364 @@
         fid.close()
         
-        points, attributes, attribute_names = read_xya(FN)
+        points, attributes, attribute_names = read_xya(FN, format = 'asc')
 
         assert attribute_names[0] == 'a1'

inundation/ga/storm_surge/pyvolution/util.py

@@ -286 +286 @@
             return res
 
-def read_xya(filename):
+def read_xya(filename, format = 'netcdf'):
     """Read simple xya file, possibly with a header in the first line, just
     x y [attributes]
     separated by whitespace
 
-    If xya is a NetCDF file it will be read otherwise attemot to read as ASCII
-    
+    Format can be either ASCII or NetCDF
+   
     Return list of points, list of attributes and
     attribute names if present otherwise None
@@ -298 +298 @@
     
     from Scientific.IO.NetCDF import NetCDFFile
-    try:
+
+    if format.lower() == 'netcdf': 
         #Get NetCDF
         
-        #FIXME: How can we suppress error message if this fails?        
         fid = NetCDFFile(filename, 'r') 
       
@@ -310 +310 @@
         #Don't close - arrays are needed outside this function,
         #alternatively take a copy here
-    except:
+    else:    
         #Read as ASCII file assuming that it is separated by whitespace
         fid = open(filename)