Changeset 631 for inundation/ga/storm_surge/pyvolution/interpolate_sww.py

Timestamp:

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

Author:

duncan

Message:

user can choose the quantity to interpolate

File:

• inundation/ga/storm_surge/pyvolution/interpolate_sww.py (modified) (10 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)