Changeset 832

Timestamp:

Feb 2, 2005, 5:53:35 PM (20 years ago)

Author:

ole

Message:

Work towards conversions for dems

Location:

inundation/ga/storm_surge/pyvolution

Files:

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

inundation/ga/storm_surge/pyvolution/data_manager.py

@@ -793 +793 @@
 
 def dem2xya(filename, verbose=False):
-    """Read Digitial Elevation model from the following ASCII format (.asc)
+    """Read Digitial Elevation model from the following NetCDF format (.asc)
 
     Example:
@@ -862 +862 @@
     #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
+    #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()
@@ -881 +881 @@
 
                                       
-
-def dem2netcdf(filename, verbose=False):
+def asciidem2netcdf(filename, verbose=False):
     """Read Digitial Elevation model from the following ASCII format (.asc)
 

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -161 +161 @@
 
     if verbose: print 'Read xya'
-    points, attributes = util.read_xya(xya_name)
+    points, attributes, _ = util.read_xya(xya_name)
 
     #Reduce number of points a bit

inundation/ga/storm_surge/pyvolution/test_data_manager.py

@@ -536 +536 @@
 
         #Convert to NetCDF xya
+        asciidem2netcdf(filename)
         dem2xya(filename)
 
         #Check contents
         #Get NetCDF
-        fid = NetCDFFile(root+'.xya.', 'r')
+        fid = NetCDFFile(root+'.xya', 'r')
       
         # Get the variables

inundation/ga/storm_surge/pyvolution/test_util.py

@@ -337 +337 @@
         
 
-    def test_xya(self):
+    def test_xya_ascii(self):
         import time, os
         FN = 'xyatest' + str(time.time()) + '.xya'
@@ -346 +346 @@
         fid.close()
         
-        points, attributes = read_xya(FN)
+        points, attributes, _ = read_xya(FN)
         
         assert allclose(points, [ [0,1], [1,0], [1,1] ])
@@ -354 +354 @@
         os.remove(FN)
         
+    def test_xya_ascii_w_names(self):
+        import time, os
+        FN = 'xyatest' + str(time.time()) + '.xya'
+        fid = open(FN, 'w')
+        fid.write('      %s %s %s\n' %('a1', 'a2', 'a3') )        
+        fid.write('%f %f %f %f %f\n' %(0,1,10,20,30) )
+        fid.write('%f %f %f %f %f\n' %(1,0,30,20,10) )
+        fid.write('%f %f %f %f %f\n' %(1,1,40.2,40.3,40.4) )            
+        fid.close()
+        
+        points, attributes, attribute_names = read_xya(FN)
+
+        assert attribute_names[0] == 'a1'
+        assert attribute_names[1] == 'a2'
+        assert attribute_names[2] == 'a3'        
+        assert allclose(points, [ [0,1], [1,0], [1,1] ])
+        assert allclose(attributes, [ [10,20,30], [30,20,10], 
+                [40.2,40.3,40.4] ])
+        
+        os.remove(FN)
         
     def test_polygon_function_constants(self):

inundation/ga/storm_surge/pyvolution/util.py

@@ -51 +51 @@
 
 
-
-#def point_on_line(point, line):
 def point_on_line(x, y, x0, y0, x1, y1): 
     """Determine whether a point is on a line segment 
@@ -61 +59 @@
         
     """
-    
-    #FIXME: Could do with some C-optimisation         
         
     from Numeric import array, dot, allclose
@@ -291 +287 @@
 
 def read_xya(filename):
-    """Read simple xya file, no header, just
+    """Read simple xya file, possibly with a header in the first line, just
     x y [attributes]
     separated by whitespace
@@ -297 +293 @@
     If xya is a NetCDF file it will be read otherwise attemot to read as ASCII
     
-    Return list of points and list of attributes
-    """
-
-
+    Return list of points, list of attributes and
+    attribute names if present otherwise None
+    """
+    
     from Scientific.IO.NetCDF import NetCDFFile
-
-
     try:
-        #Get NetCDF        
-        fid = NetCDFFile(filename, 'r')
+        #Get NetCDF
+        
+        #FIXME: How can we suppress error message if this fails?        
+        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
+        attribute_names = fid.variables['attribute_names']
+        #Don't close - arrays are needed outside this function,
+        #alternatively take a copy here
     except:
-        #Read as ASCII
+        #Read as ASCII file assuming that it is separated by whitespace
         fid = open(filename)
         lines = fid.readlines()
+        fid.close()
+
+        #Check if there is a header line
+        fields = lines[0].strip().split()
+        try:
+            float(fields[0])
+        except:
+            #This must be a header line
+            attribute_names = fields
+            lines = lines[1:]
+        else:
+            attribute_names = None
     
         points = []     
@@ -323 +331 @@
         for line in lines:
             fields = line.strip().split()
-            points.append( (float(fields[0]),  float(fields[1])) )
+            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]
+        
+    return points, attributes, attribute_names
     
 
@@ -407 +414 @@
     py = polygon[:,1]   
 
-    
     #Begin algorithm
     indices = []
@@ -420 +426 @@
             
             #Check for case where point is contained in line segment    
-            ##if point_on_line( (x,y), [ [px[i], py[i]], [px[j], py[j]] ]):
             if point_on_line(x, y, px[i], py[i], px[j], py[j]):
                 if closed: