Changeset 631

Timestamp:

Nov 26, 2004, 12:23:25 PM (19 years ago)

Author:

duncan

Message:

user can choose the quantity to interpolate

Location:

inundation/ga/storm_surge/pyvolution

Files:

• interpolate_sww.py (modified) (10 diffs)
• test_interpolate_sww.py (modified) (4 diffs)

inundation/ga/storm_surge/pyvolution/interpolate_sww.py

@@ -23 +23 @@
 from least_squares import Interpolation 
 
+DEFAULT_QUANTITY = "height"
 
 class Interpolate_sww:
-    def __init__(self,file_name):
+    def __init__(self,file_name, quantity_name):
 
+        #It's bad to have the quantity_name passed in here.
+        # But it works with how the program is currently used.
+        # Refactor when necessary. - DSG
         
-        x, y,z, volumes, time, stage = self.read_sww(file_name)       
-
+        x, y, volumes, time, quantity = self.read_sww(file_name, quantity_name)
         vertex_coordinates = transpose([x,y])
         
@@ -40 +43 @@
             print y
             print "****************************"
-            print "Z "
-            print z
-            print "****************************"
             print "V "
             print volumes
@@ -49 +49 @@
             print time
             print "****************************"
-            print "Stage "
-            print stage
+            print "quantity "
+            print quantity
             print "****************************"
             print "vertex_coordinates"
@@ -58 +58 @@
         self.interp = Interpolation(vertex_coordinates, volumes )
         self.time = time
-        self.stage = transpose(stage)
-        self.bed_elevation = z
-        #print "self.stage",self.stage
+
+        self.quantity_name = quantity_name
+        self.quantity = quantity 
         
     def interpolate_xya(self,file_name):
@@ -70 +70 @@
         file_name - the xya file
         """
-        #
         self.point_coordinates = self.read_xya(file_name)
         self.interp.build_interpolation_matrix_A(self.point_coordinates)
-
-        #Change this if you want to generalise this application -
-        # This is height specific
-        point_stage =  self.interp.interpolate( self.stage)
-        point_bed_elevation = self.interp.interpolate( self.bed_elevation)
-        self.depth = transpose(transpose(point_stage) - point_bed_elevation)
+        self.interpolated_quantity = self.interp.interpolate( transpose(self.quantity))
         
-        #print "self.depth", self.depth
-        
-    def read_sww(self,file_name):
+    def read_sww(self,file_name, quantity_name):
         """
         Read in an sww file.
@@ -101 +93 @@
         #FIXME Have this reader as part of data_manager?
 
-        #import time, os
-        #from Numeric import array, zeros, allclose, Float, concatenate
         from Scientific.IO.NetCDF import NetCDFFile             
             
@@ -112 +102 @@
         x = fid.variables['x'][:]
         y = fid.variables['y'][:]
-        z = fid.variables['z'][:] 
         volumes = fid.variables['volumes'][:] 
         time = fid.variables['time'][:]
-        stage = fid.variables['stage'][:,:]   # 2D 
-        
+        try:
+            if quantity_name == 'height':
+                z = fid.variables['z'][:]
+                stage = fid.variables['stage'][:,:]
+                quantity = stage - z  # 2D, using broadcasting 
+            else:
+                quantity = fid.variables[quantity_name][:,:]   # 2D 
+        except (KeyError, IndexError),e:
+            fid.close()
+            raise KeyError
+            
         fid.close()
-        
-        return x, y, z, volumes, time, stage
+        return x, y, volumes, time, quantity
 
     def read_xya(self,point_file):
@@ -162 +159 @@
         xya_dict = {}
         xya_dict['pointlist'] = self.point_coordinates
-        xya_dict['pointattributelist'] = self.depth
+        xya_dict['pointattributelist'] = self.interpolated_quantity
         
         export_xya_file(point_file, xya_dict, title, delimiter = delimiter)
@@ -172 +169 @@
     """
     import os, sys
-    usage = "usage: %s pyvolution_results.sww points.xya depth.xya" %         os.path.basename(sys.argv[0])
-
+    usage = "usage: %s pyvolution_results.sww points.xya depth.xya [height|stage|(other quantities)]" %         os.path.basename(sys.argv[0])
     if len(sys.argv) < 4:
         print usage
@@ -180 +176 @@
         point_file_in = sys.argv[2]
         point_file_out = sys.argv[3]
-        
-        interp = Interpolate_sww(sww_file)
+        if len(sys.argv) == 5:
+            quantity_name = sys.argv[4]
+        else:
+            quantity_name = DEFAULT_QUANTITY   
+        #print "quantity",quantity
+        try:
+            interp = Interpolate_sww(sww_file, quantity_name)
+        except KeyError:
+            print "Error: Unknown quantity"
+            sys.exit(1)
+
         interp.interpolate_xya(point_file_in)
         interp.write_depth_xya(point_file_out)

inundation/ga/storm_surge/pyvolution/test_interpolate_sww.py

@@ -159 +159 @@
 
         #print "self.domain.filename",self.domain.filename
-        interp = Interpolate_sww(sww.filename)
+        interp = Interpolate_sww(sww.filename, 'height')
         
         assert allclose(interp.time,[0.0,2.0])
-        answer = [ 0.15, 0.1, 0., -0.3, -0.35, -0.4, -0.7, -0.8, -0.850]
+
+        #answer = [ 0.15, 0.1, 0., -0.3, -0.35, -0.4, -0.7, -0.8, -0.850]
         #print "answer",answer
         #print interp.stage[0]
-        stage_t = transpose(interp.stage)
-        assert allclose(stage_t[0], answer)
-        assert allclose(stage_t[1],stage_t[0])
+        #stage_t = transpose(interp.stage)
+        #assert allclose(stage_t[0], answer)
+        #assert allclose(stage_t[1],stage_t[0])
 
         # create an .xya file
@@ -177 +178 @@
         interp.interpolate_xya(point_file)
 
-
         answer = [[0.08, 0.08], [0.02, 0.02], [0.14, 0.14], [.08,.08]]
         #print "answer",answer
-        assert allclose(interp.depth,answer)
+        assert allclose(interp.interpolated_quantity,answer)
 
         # create an output .xya file
@@ -190 +190 @@
 
         assert allclose(interp.point_coordinates, xya_dict['pointlist'])
-        assert allclose(interp.depth,
+        assert allclose(interp.interpolated_quantity,
                         xya_dict['pointattributelist'] )
         
@@ -198 +198 @@
         time_list = []
 
-        # this is if titles start with x,y
-        #answer = string_list.pop(0)
-        #self.failUnless( answer == 'x', 'Title is wrong!')
-        #self.failUnless( string_list.pop(0) == 'y', 'Title is wrong!')
-        
-        for time in string_list:
-            time_list.append(float(time))
-        #print "interp.time", interp.time
-        #print "time_list", time_list
-        assert allclose(interp.time,
-                        time_list)
-        
-        
+        # Try another quantity      
+        interp = Interpolate_sww(sww.filename, 'stage')
+        interp.interpolate_xya(point_file)
+
+        answer = [[0.08, 0.08], [0.02, 0.02], [0.14, 0.14], [-.32,-.32]]
+        #print "answer",answer
+        assert allclose(interp.interpolated_quantity,answer)
+
+        
+        # look at error catching
+        try:
+            interp = Interpolate_sww(sww.filename, 'funky!')
+        except KeyError:
+            pass
+        else:
+            self.failUnless(0==1,
+                        'bad key did not raise an error!')
+        
+        # look at error catching
+        try:
+            interp = Interpolate_sww(sww.filename, 'z')
+        except KeyError:
+            pass
+        else:
+            self.failUnless(0==1,
+                        'bad key did not raise an error!')
+                
         #Cleanup
         os.remove(sww.filename)