Changeset 356

Timestamp:

Oct 6, 2004, 10:50:22 AM (20 years ago)

Author:

duncan

Message:

minor work on interpolate sww

Location:

inundation/ga/storm_surge/pyvolution

Files:

• interpolate_sww.py (modified) (5 diffs)
• least_squares.py (modified) (2 diffs)
• test_interpolate_sww.py (modified) (3 diffs)

inundation/ga/storm_surge/pyvolution/interpolate_sww.py

@@ -1 +1 @@
 """ Interpolation of a sww file.
+Used to interpolate height information from sww files.
 
+When using as a stand alone application,
+Input;
+ - The sww file with stage and bed elevation information
+ - An xya file specifying the points (x,y) where the height values
+   (w.r.t. time) are needed.
+
+Ouput;
+An xya file with x,y and height values w.r.t. time.
+
+The first row of the output xya file has the time value for each height
+Each following row has the format x,y,[height,]
+
+NOTE: stage = bed elevation + height
 
    Ole Nielsen, Stephen Roberts, Duncan Gray, Christopher Zoppou
@@ -6 +20 @@
 """
 
-#from general_mesh import General_Mesh
-#from Numeric import zeros, array, Float, Int, dot, transpose
-#from LinearAlgebra import solve_linear_equations
+from Numeric import transpose
+from least_squares import Interpolation 
 
-
-from Numeric import transpose
-
-# def read_sww(file_name):
-#         import time, os
-#         from Numeric import array, zeros, allclose, Float, concatenate
-#         from Scientific.IO.NetCDF import NetCDFFile           
-
-#         #Check contents
-#         #Get NetCDF
-#         fid = NetCDFFile(file_name, 'r') 
-      
-#         # Get the variables
-#         x = fid.variables['x'][:]
-#         y = fid.variables['y'][:]
-#         # z = fid.variables['z'][:] # We don't care about the bed elevation
-#         time = fid.variables['time'][:]
-#         stage = fid.variables['stage'][:,:]   # 2D   
-        
-        
-#         return x,y,time,stage
 
 class Interpolate_sww:
     def __init__(self,file_name):
-        from least_squares import Interpolation 
 
         
-        x, y, volumes, time, stage = self.read_sww(file_name)       
+        x, y,z, volumes, time, stage = self.read_sww(file_name)       
 
         vertex_coordinates = transpose([x,y])
         
-        if True:
+        if False:
             print "****************************"
-            print "x ",x
+            print "x "
+            print x
             print "****************************"
-            print "Y ",y
+            print "Y "
+            print y
             print "****************************"
-            print "V ",volumes
+            print "Z "
+            print z
             print "****************************"
-            print "Time ", time
+            print "V "
+            print volumes
             print "****************************"
-            print "Stage ",stage
+            print "Time "
+            print time
             print "****************************"
-            print "vertex_coordinates",vertex_coordinates
+            print "Stage "
+            print stage
+            print "****************************"
+            print "vertex_coordinates"
+            print vertex_coordinates
             print "****************************"
 
-        interp = Interpolation(vertex_coordinates, volumes )
+        self.interp = Interpolation(vertex_coordinates, volumes )
+        self.time = time
+        self.stage = transpose(stage)
+        self.bed_elevation = z
+        #print "self.stage",self.stage
         
+    def interpolate_xya(self,file_name):
+        """
+        Given a point file, interpolate the height w.r.t. time at the points
+        specified in the point file.
 
+        Input;
+        file_name - the xya file
+        """
+        #
+        self.point_coordinates = self.read_xya(file_name)
+        self.interp.build_interpolation_matrix_A(self.point_coordinates)
+        #print "self.interp.A",self.interp.A
+
+        #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.point_heights = transpose(transpose(point_stage) - point_bed_elevation)
+        
+        #print "self.point_heights", self.point_heights
+        
     def read_sww(self,file_name):
+        """
+        Read in an sww file.
+
+        Input;
+        file_name - the sww file
+
+        Output;
+        x - Vector of x values
+        y - Vector of y values
+        z - Vector of bed elevation
+        volumes - Array.  Each row has 3 values, representing
+                  the vertices that define the volume
+        time - Vector of the times where there is stage information
+        stage - array with respect to time and vertices (x,y)
+        """
         
-    #FIXME Have this reader as part of data_manager?
+        #FIXME Have this reader as part of data_manager?
 
         #import time, os
@@ -74 +113 @@
         x = fid.variables['x'][:]
         y = fid.variables['y'][:]
-        # z = fid.variables['z'][:] # We don't care about the bed elevation
+        z = fid.variables['z'][:] 
         volumes = fid.variables['volumes'][:] 
         time = fid.variables['time'][:]
@@ -81 +120 @@
         fid.close()
         
-        return x, y, volumes, time, stage
+        return x, y, z, volumes, time, stage
+
+    def read_xya(self,point_file):
+        """
+        Read in an sww file.
+
+        Input;
+        point_file - the xya file
+
+        Output;
+        point_coordinates - the points (x,y) where the height values
+                             (w.r.t. time) are needed
+        """
+        from load_mesh.loadASCII import  load_xya_file
+        
+        point_dict = load_xya_file(point_file)
+        return point_dict['pointlist']
 
 
+    def write_point_heights_xya(self,point_file, delimiter = ','):
+        """
+        pre condition:
+        point attributes have been determined
+        (interpolate_xya has been called)
+        The time list is defined
+        """
+        from load_mesh.loadASCII import  export_xya_file
 
+        #create the first line, giving the times 
+        if len(self.time) != 0:
+            for i in range(len(self.time)):
+                if i == 0:
+                    title = str(self.time[0])
+                else:
+                    title = title + delimiter + str(self.time[i])
+        else:
+            title  = ""
+        
+        #print ">" +  str(self.time) + "<"
+        #print ">" + title + "<"
+        
+        xya_dict = {}
+        xya_dict['pointlist'] = self.point_coordinates
+        xya_dict['pointattributelist'] = self.point_heights
+        
+        export_xya_file( point_file, xya_dict, title, delimiter = delimiter)
         
 #-------------------------------------------------------------
 if __name__ == "__main__":
     """
-    Load in an sww file and an xya file and return ?? 
+    Load in an sww file and an xya file and return an xya file 
     """
     import os, sys
-    usage = "usage: %s pyvolution_results.sww point.xya " %         os.path.basename(sys.argv[0])
+    usage = "usage: %s pyvolution_results.sww point.xya point_attributes_out.xya" %         os.path.basename(sys.argv[0])
 
-    if len(sys.argv) < 2:
+    if len(sys.argv) < 4:
         print usage
     else:
         sww_file = sys.argv[1]
-        point_file = sys.argv[2]
-        #mesh_output_file = sys.argv[3]
-        #fit_to_mesh_file(mesh_file, point_file, mesh_output_file, alpha)
-        c = Interpolate_sww(sww_file)
+        point_file_in = sys.argv[2]
+        point_file_out = sys.argv[3]
+        
+        interp = Interpolate_sww(sww_file)
+        interp.interpolate_xya(point_file_in)
+        interp.write_point_heights_xya(point_file_out)
 
 
@@ -111 +194 @@
 
 
-
-

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -25 +25 @@
 from LinearAlgebra import solve_linear_equations
 
+try:
+    from util import gradient
+
+except ImportError, e: 
+    #FIXME reduce the dependency of modules in pyvolution
+    def gradient(x0, y0, x1, y1, x2, y2, q0, q1, q2):
+        """
+        """
+    
+        det = (y2-y0)*(x1-x0) - (y1-y0)*(x2-x0)            
+        a = (y2-y0)*(q1-q0) - (y1-y0)*(q2-q0)
+        a /= det
+
+        b = (x1-x0)*(q2-q0) - (x2-x0)*(q1-q0)
+        b /= det            
+
+        return a, b
+
 ALPHA_INITIAL = 0.001
 
@@ -238 +256 @@
         #"element stiffness matrices:
 
-        from util import gradient
         
         m = self.mesh.coordinates.shape[0] #Nbr of basis functions (1/vertex)

inundation/ga/storm_surge/pyvolution/test_interpolate_sww.py

@@ -14 +14 @@
 from data_manager import *
 #from config import epsilon
+
 
 class dataTestCase(unittest.TestCase):
@@ -136 +137 @@
  
     def test_interpolate_sww(self):
-        """Not a test, rather a look at the sww format
+        """Not reaa unit test, rather a system test for 
         """
 
         import time, os
-        from Numeric import array, zeros, allclose, Float, concatenate
-        from Scientific.IO.NetCDF import NetCDFFile        
+        from Numeric import array, zeros, allclose, Float, concatenate, \
+             transpose
+        from Scientific.IO.NetCDF import NetCDFFile    
+        import tempfile
+        from load_mesh.loadASCII import  load_xya_file     
 
         self.domain.filename = 'datatest' + str(time.time())
@@ -154 +158 @@
         sww.store_timestep('level')
 
-        print "self.domain.filename",self.domain.filename
-        c = Interpolate_sww(sww.filename)
+        #print "self.domain.filename",self.domain.filename
+        interp = Interpolate_sww(sww.filename)
+        
+        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]
+        #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])
+
+        # create an .xya file
+        point_file = tempfile.mktemp(".xya")
+        fd = open(point_file,'w')
+        fd.write("# demo \n 0.0, 0.6,2.,4 \n 0.0, 0.9,4,8 \n 0.0,0.1,4.,8 \n 0.4,1.0,4.,8 \n")
+        fd.close()
+
+        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.point_heights,answer)
+
+        # create an output .xya file
+        point_file_out = tempfile.mktemp(".xya")
+        interp.write_point_heights_xya(point_file_out)
+        
+        #check the output file
+        xya_dict = load_xya_file(point_file_out)
+
+        assert allclose(interp.point_coordinates, xya_dict['pointlist'])
+        assert allclose(interp.point_heights,
+                        xya_dict['pointattributelist'] )
+        
+        title = xya_dict['title'].split('\n')
+        #print "title",title
+        string_list = title[0].split(',') # assume a title has only one line
+        time_list = []
+        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)
+        
         
         #Cleanup
         os.remove(sww.filename)
+        os.remove(point_file_out)
  
 #-------------------------------------------------------------