Changeset 3411

Timestamp:

Jul 25, 2006, 12:14:42 PM (19 years ago)

Author:

duncan

Message:

adding sww to csv function

Location:

inundation/fit_interpolate

Files:

• interpolate.py (modified) (4 diffs)
• test_interpolate.py (modified) (2 diffs)

inundation/fit_interpolate/interpolate.py

@@ -21 +21 @@
 import os
 from warnings import warn
+from math import sqrt
+from csv import writer, DictWriter
 
 from Numeric import zeros, array, Float, Int, dot, transpose, concatenate, \
@@ -33 +35 @@
 from utilities.numerical_tools import ensure_numeric, mean, INF
 from utilities.polygon import in_and_outside_polygon
-from geospatial_data.geospatial_data import Geospatial_data
+from geospatial_data.geospatial_data import Geospatial_data, ensure_absolute
 from fit_interpolate.search_functions import search_tree_of_vertices
 from fit_interpolate.general_fit_interpolate import FitInterpolate
-
+from pyvolution.util import file_function
 
 
@@ -276 +278 @@
         return A
 
+def interpolate_sww2csv(sww_file, points,
+                        depth_file,
+                        velocity_file,
+                                #quantities = ['depth', 'velocity'],
+                                verbose=True,
+                                 use_cache = True):
+    """
+    Interpolate the quantities at a given set of locations, given
+    an sww file.
+    The results are written to a csv file.
+
+    In the future let points be a csv or xya file.
+    And the user choose the quantities.
+    """
+    quantities =  ['stage', 'elevation', 'xmomentum', 'ymomentum']
+    #print "points",points 
+    points = ensure_absolute(points)
+    point_count = len(points)
+    callable_sww = file_function(sww_file,
+                                 quantities=quantities,
+                                 interpolation_points=points,
+                                 verbose=verbose,
+                                 use_cache=use_cache)
+    
+    depth_writer = writer(file(depth_file, "wb"))
+    velocity_writer = writer(file(velocity_file, "wb"))
+    # Write heading
+    depth_writer.writerow([str(x[0])+ ':' + str(x[1]) for x in points])
+    velocity_writer.writerow([str(x[0])+ ':' + str(x[1]) for x in points])
+    #for row in someiterable:
+        #csvwriter.writerow(row)
+    for time in callable_sww.get_time():
+        depths = []
+        velocitys = []
+        for point_i, point in enumerate(points):
+            quantities = callable_sww(time,point_i)
+            #print "quantities", quantities
+            
+            w = quantities[0]
+            z = quantities[1]
+            uh = quantities[2]
+            vh = quantities[3]
+            depth = w - z 
+              
+            if w == INF or z == INF or uh == INF or vh == INF:
+                velocity = INF
+            else:
+                momentum = sqrt(uh*uh + vh*vh)
+                if depth > 1.e-30: # use epsilon
+                    velocity = momentum / depth  #Absolute velocity
+                else:
+                    velocity = 0
+                    
+            print "depth",depth
+            print "velocity",velocity 
+            depths.append(depth)
+            velocitys.append(velocity)
+        depth_writer.writerow(depths)
+        velocity_writer.writerow(velocitys)
+
 class Interpolation_function:
     """Interpolation_interface - creates callable object f(t, id) or f(t,x,y)
@@ -715 +777 @@
 #-------------------------------------------------------------
 if __name__ == "__main__":
-    a = [0.0, 0.0]
-    b = [0.0, 2.0]
-    c = [2.0,0.0]
-    points = [a, b, c]
-    vertices = [ [1,0,2] ]   #bac
-    
-    data = [ [2.0/3, 2.0/3] ] #Use centroid as one data point
-    
-    interp = Interpolate(points, vertices) #, data)
-    A  = interp._build_interpolation_matrix_A(data, verbose=True)
-    A = A.todense()
-    print "A",A
-    assert allclose(A, [[1./3, 1./3, 1./3]])
-    print "finished"
+    names = ["x","y"]
+    someiterable = [[1,2],[3,4]]
+    csvwriter = writer(file("some.csv", "wb"))
+    csvwriter.writerow(names)
+    for row in someiterable:
+        csvwriter.writerow(row)

inundation/fit_interpolate/test_interpolate.py

@@ -23 +23 @@
 from data_manager import get_dataobject
 from geospatial_data.geospatial_data import Geospatial_data
+from pmesh.mesh import Mesh
 
 def distance(x, y):
@@ -1351 +1352 @@
             self.failUnless( z[i,1] == answer[11,1], 'Fail!')
             self.failUnless( z[i,0] == answer[11,0], 'Fail!')
-        
+
+    def ztest_interpolate_sww2csv(self):
+
+        def elevation_function(x, y):
+            return -x
+        
+        # create mesh
+        mesh_file = tempfile.mktemp(".tsh")    
+        points = [[0.0,0.0],[6.0,0.0],[6.0,6.0],[0.0,6.0]]
+        m = Mesh()
+        m.add_vertices(points)
+        m.autoSegment()
+        m.generate_mesh(verbose=False)
+        m.export_mesh_file(mesh_file)
+        
+        #Create shallow water domain
+        domain = Domain(mesh_file)
+        os.remove(mesh_file)
+        
+        domain.default_order=2
+        domain.beta_h = 0
+
+        #Set some field values
+        domain.set_quantity('elevation', elevation_function)
+        domain.set_quantity('friction', 0.03)
+        domain.set_quantity('xmomentum', 10.0)
+        domain.set_quantity('ymomentum', 100.0)
+
+        ######################
+        # Boundary conditions
+        B = Transmissive_boundary(domain)
+        domain.set_boundary( {'exterior': B})
+
+        # This call mangles the stage values.
+        domain.distribute_to_vertices_and_edges()
+        domain.set_quantity('stage', 1.0)
+
+        #sww_file = tempfile.mktemp("")
+        domain.filename = 'datatest' + str(time.time())
+        #domain.filename = sww_file
+        #print "domain.filename",domain.filename 
+        domain.format = 'sww'
+        domain.smooth = True
+        domain.reduction = mean
+
+        sww = get_dataobject(domain)
+        sww.store_connectivity()
+        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
+        domain.set_quantity('stage', -1.0)
+        domain.time = 2.
+        sww.store_timestep(['stage', 'xmomentum', 'ymomentum'])
+
+        # test the function
+        points = [[5.0,1.],[0.5,2.]]
+        interpolate_sww2csv(sww.filename, points, "aa.csv", "bb.csv")
+
+        print "sww.filename",sww.filename 
+        #os.remove(sww.filename)
 #-------------------------------------------------------------
 if __name__ == "__main__":
 
-    suite = unittest.makeSuite(Test_Interpolate,'test')
-    #suite = unittest.makeSuite(Test_Interpolate,'test_interpolation_precompute_points')
+    #suite = unittest.makeSuite(Test_Interpolate,'test')
+    suite = unittest.makeSuite(Test_Interpolate,'test_interpolate_sww2csv')
     runner = unittest.TextTestRunner(verbosity=1)
     runner.run(suite)