Changeset 826 for inundation/ga/storm_surge/pyvolution/data_manager.py

Timestamp:

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

Author:

ole

Message:

Played with DEM's and NetCDF

File:

• inundation/ga/storm_surge/pyvolution/data_manager.py (modified) (1 diff)

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