Changeset 2009 for inundation/pyvolution/data_manager.py

Timestamp:

Nov 7, 2005, 2:21:12 PM (19 years ago)

Author:

ole

Message:

Added conversion from HEC-RAS to NetCDF pts + test

File:

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

inundation/pyvolution/data_manager.py

@@ -1203 +1203 @@
 
     infile.close()
+    outfile.close()
+
+
+
+def _read_hecras_cross_sections(lines):
+    """Return block of surface lines for each cross section
+    Starts with SURFACE LINE,
+    Ends with END CROSS-SECTION  
+    """
+
+    points = []
+
+    reading_surface = False
+    for i, line in enumerate(lines): 
+        
+        if len(line.strip()) == 0:    #Ignore blanks
+            continue                       
+
+        if lines[i].strip().startswith('SURFACE LINE'):
+            reading_surface = True
+            continue
+
+        if lines[i].strip().startswith('END') and reading_surface:
+            yield points
+            reading_surface = False
+            points = []
+
+        if reading_surface:
+            fields = line.strip().split(',')
+            easting = float(fields[0])
+            northing = float(fields[1])
+            elevation = float(fields[2])                        
+            points.append([easting, northing, elevation])                
+
+
+
+
+def hecras_cross_sections2pts(basename_in,
+                              basename_out=None,
+                              verbose=False):
+    """Read HEC-RAS Elevation datal from the following ASCII format (.sdf)
+    
+    Example:
+    
+
+# RAS export file created on Mon 15Aug2005 11:42
+# by HEC-RAS Version 3.1.1
+
+BEGIN HEADER:
+  UNITS: METRIC
+  DTM TYPE: TIN
+  DTM: v:\1\cit\perth_topo\river_tin
+  STREAM LAYER: c:\local\hecras\21_02_03\up_canning_cent3d.shp
+  CROSS-SECTION LAYER: c:\local\hecras\21_02_03\up_can_xs3d.shp
+  MAP PROJECTION: UTM
+  PROJECTION ZONE: 50
+  DATUM: AGD66
+  VERTICAL DATUM: 
+  NUMBER OF REACHES:  19 
+  NUMBER OF CROSS-SECTIONS:  14206 
+END HEADER:
+
+
+Only the SURFACE LINE data of the following form will be utilised
+
+  CROSS-SECTION:
+    STREAM ID:Southern-Wungong
+    REACH ID:Southern-Wungong
+    STATION:19040.* 
+    CUT LINE:
+      405548.671603161 , 6438142.7594925 
+      405734.536092045 , 6438326.10404912 
+      405745.130459356 , 6438331.48627354 
+      405813.89633823 , 6438368.6272789 
+    SURFACE LINE:
+     405548.67,   6438142.76,   35.37
+     405552.24,   6438146.28,   35.41
+     405554.78,   6438148.78,   35.44
+     405555.80,   6438149.79,   35.44
+     405559.37,   6438153.31,   35.45
+     405560.88,   6438154.81,   35.44
+     405562.93,   6438156.83,   35.42
+     405566.50,   6438160.35,   35.38
+     405566.99,   6438160.83,   35.37
+     ...
+   END CROSS-SECTION  
+
+    Convert to NetCDF pts format which is
+
+    points:  (Nx2) Float array
+    elevation: N Float array
+    """
+
+    #FIXME: Can this be written feasibly using write_pts?
+
+    import os
+    from Scientific.IO.NetCDF import NetCDFFile
+    from Numeric import Float, zeros, reshape
+
+    root = basename_in
+
+    #Get ASCII file
+    infile = open(root + '.sdf', 'r')  #Open SDF file for read
+
+    if verbose: print 'Reading DEM from %s' %(root + '.sdf')
+
+    lines = infile.readlines()
+    infile.close()
+
+    if verbose: print 'Converting to pts format'
+
+    i = 0
+    while lines[i].strip() == '' or lines[i].strip().startswith('#'):
+        i += 1
+
+    assert lines[i].strip().upper() == 'BEGIN HEADER:'
+    i += 1
+    
+    assert lines[i].strip().upper().startswith('UNITS:')
+    units = lines[i].strip().split()[1]
+    i += 1    
+
+    assert lines[i].strip().upper().startswith('DTM TYPE:')    
+    i += 1
+
+    assert lines[i].strip().upper().startswith('DTM:')    
+    i += 1            
+
+    assert lines[i].strip().upper().startswith('STREAM')    
+    i += 1
+
+    assert lines[i].strip().upper().startswith('CROSS')    
+    i += 1
+
+    assert lines[i].strip().upper().startswith('MAP PROJECTION:')
+    projection = lines[i].strip().split(':')[1]
+    i += 1
+
+    assert lines[i].strip().upper().startswith('PROJECTION ZONE:')
+    zone = int(lines[i].strip().split(':')[1])
+    i += 1
+
+    assert lines[i].strip().upper().startswith('DATUM:')
+    datum = lines[i].strip().split(':')[1]
+    i += 1
+
+    assert lines[i].strip().upper().startswith('VERTICAL DATUM:')
+    i += 1
+
+    assert lines[i].strip().upper().startswith('NUMBER OF REACHES:') 
+    i += 1
+
+    assert lines[i].strip().upper().startswith('NUMBER OF CROSS-SECTIONS:')
+    number_of_cross_sections = int(lines[i].strip().split(':')[1])
+    i += 1                        
+
+
+    #Now read all points
+    points = []
+    elevation = []
+    for j, entries in enumerate(_read_hecras_cross_sections(lines[i:])):
+        for k, entry in enumerate(entries):
+            points.append(entry[:2])
+            elevation.append(entry[2])            
+
+
+    msg = 'Actual #number_of_cross_sections == %d, Reported as %d'\
+          %(j+1, number_of_cross_sections)          
+    assert j+1 == number_of_cross_sections, msg
+    
+    #Get output file
+    if basename_out == None:
+        ptsname = root + '.pts'
+    else:
+        ptsname = basename_out + '.pts'
+
+    if verbose: print 'Store to NetCDF file %s' %ptsname
+    # NetCDF file definition
+    outfile = NetCDFFile(ptsname, 'w')
+
+    #Create new file
+    outfile.institution = 'Geoscience Australia'
+    outfile.description = 'NetCDF pts format for compact and portable ' +\
+                          'storage of spatial point data derived from HEC-RAS'
+    
+    #Georeferencing
+    outfile.zone = zone
+    outfile.xllcorner = 0.0
+    outfile.yllcorner = 0.0
+    outfile.false_easting = 500000    #FIXME: Use defaults from e.g. config.py
+    outfile.false_northing = 1000000
+
+    outfile.projection = projection
+    outfile.datum = datum
+    outfile.units = units
+
+
+    outfile.createDimension('number_of_points', len(points))
+    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
+
+    # variable definitions
+    outfile.createVariable('points', Float, ('number_of_points',
+                                             'number_of_dimensions'))
+    outfile.createVariable('elevation', Float, ('number_of_points',))
+
+    # Get handles to the variables
+    outfile.variables['points'][:] = points
+    outfile.variables['elevation'][:] = elevation
+
     outfile.close()