Changeset 2009

Timestamp:

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

Author:

ole

Message:

Added conversion from HEC-RAS to NetCDF pts + test

Location:

inundation/pyvolution

Files:

• data_manager.py (modified) (1 diff)
• test_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()
 

inundation/pyvolution/test_data_manager.py

@@ -739 +739 @@
         os.remove(root + '.asc')
         os.remove(root + '.prj')
+
+
+
+    def test_hecras_cross_sections2pts(self):
+        """Test conversion from HECRAS cross sections in ascii format
+        to native NetCDF pts format
+        """
+
+        import time, os
+        from Numeric import array, zeros, allclose, Float, concatenate
+        from Scientific.IO.NetCDF import NetCDFFile
+
+        #Write test asc file
+        root = 'hecrastest'
+
+        filename = root+'.sdf'
+        fid = open(filename, 'w')
+        fid.write("""
+# 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:\\x_local\hecras\21_02_03\up_canning_cent3d.shp
+  CROSS-SECTION LAYER: c:\\x_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:  2 
+END HEADER:
+
+
+BEGIN CROSS-SECTIONS:
+
+  CROSS-SECTION:
+    STREAM ID:Southern-Wungong
+    REACH ID:Southern-Wungong
+    STATION:21410   
+    CUT LINE:
+      407546.08 , 6437277.542 
+      407329.32 , 6437489.482 
+      407283.11 , 6437541.232 
+    SURFACE LINE:
+     407546.08,   6437277.54,   52.14
+     407538.88,   6437284.58,   51.07
+     407531.68,   6437291.62,   50.56
+     407524.48,   6437298.66,   49.58
+     407517.28,   6437305.70,   49.09
+     407510.08,   6437312.74,   48.76
+  END:
+
+  CROSS-SECTION:
+    STREAM ID:Swan River      
+    REACH ID:Swan Mouth      
+    STATION:840.*   
+    CUT LINE:
+      381178.0855 , 6452559.0685 
+      380485.4755 , 6453169.272 
+    SURFACE LINE:
+     381178.09,   6452559.07,   4.17
+     381169.49,   6452566.64,   4.26
+     381157.78,   6452576.96,   4.34
+     381155.97,   6452578.56,   4.35
+     381143.72,   6452589.35,   4.43
+     381136.69,   6452595.54,   4.58
+     381114.74,   6452614.88,   4.41
+     381075.53,   6452649.43,   4.17
+     381071.47,   6452653.00,   3.99
+     381063.46,   6452660.06,   3.67
+     381054.41,   6452668.03,   3.67
+  END:  
+END CROSS-SECTIONS:
+""")
+
+        fid.close()
+
+
+        #Convert to NetCDF pts
+        hecras_cross_sections2pts(root)
+
+        #Check contents
+        #Get NetCDF
+        fid = NetCDFFile(root+'.pts', 'r')
+
+        # Get the variables
+        #print fid.variables.keys()
+        points = fid.variables['points']
+        elevation = fid.variables['elevation']
+
+        #Check values
+        ref_points = [[407546.08, 6437277.54],
+                      [407538.88, 6437284.58],
+                      [407531.68, 6437291.62],
+                      [407524.48, 6437298.66],
+                      [407517.28, 6437305.70],
+                      [407510.08, 6437312.74]]
+
+        ref_points += [[381178.09, 6452559.07], 
+                       [381169.49, 6452566.64], 
+                       [381157.78, 6452576.96], 
+                       [381155.97, 6452578.56], 
+                       [381143.72, 6452589.35], 
+                       [381136.69, 6452595.54], 
+                       [381114.74, 6452614.88], 
+                       [381075.53, 6452649.43], 
+                       [381071.47, 6452653.00], 
+                       [381063.46, 6452660.06], 
+                       [381054.41, 6452668.03]] 
+        
+                      
+        ref_elevation = [52.14, 51.07, 50.56, 49.58, 49.09, 48.76]
+        ref_elevation += [4.17, 4.26, 4.34, 4.35, 4.43, 4.58, 4.41, 4.17, 3.99, 3.67, 3.67]
+
+        #print points[:]
+        #print ref_points
+        assert allclose(points, ref_points)
+
+        #print attributes[:]
+        #print ref_elevation
+        assert allclose(elevation, ref_elevation)
+
+        #Cleanup
+        fid.close()
+
+
+        os.remove(root + '.sdf')
+        os.remove(root + '.pts')