Changeset 7693

Timestamp:

Apr 23, 2010, 4:46:14 PM (14 years ago)

Author:

hudson

Message:

ticket 292: sww2csv_gauges needs to ignore gauges that are not in the domain.
Now, a warning message is output, instead of generating a maningless file.

Location:

anuga_core/source/anuga/abstract_2d_finite_volumes

Files:

• gauge.py (modified) (3 diffs)
• test_gauge.py (modified) (1 diff)
• util.py (modified) (2 diffs)

anuga_core/source/anuga/abstract_2d_finite_volumes/gauge.py

@@ -23 +23 @@
 from math import sqrt
 
+def _quantities2csv(quantities, point_quantities, centroids, point_i):
+    points_list = []
+    
+    for quantity in quantities:
+        #core quantities that are exported from the interpolator     
+        if quantity == 'stage':
+            points_list.append(point_quantities[0])
+            
+        if quantity == 'elevation':
+            points_list.append(point_quantities[1])
+            
+        if quantity == 'xmomentum':
+            points_list.append(point_quantities[2])
+            
+        if quantity == 'ymomentum':
+            points_list.append(point_quantities[3])
+
+        #derived quantities that are calculated from the core ones
+        if quantity == 'depth':
+            points_list.append(point_quantities[0] 
+                               - point_quantities[1])
+
+        if quantity == 'momentum':
+            momentum = sqrt(point_quantities[2]**2 
+                            + point_quantities[3]**2)
+            points_list.append(momentum)
+            
+        if quantity == 'speed':
+            #if depth is less than 0.001 then speed = 0.0
+            if point_quantities[0] - point_quantities[1] < 0.001:
+                vel = 0.0
+            else:
+                if point_quantities[2] < 1.0e6:
+                    momentum = sqrt(point_quantities[2]**2
+                                    + point_quantities[3]**2)
+                    vel = momentum / (point_quantities[0] 
+                                      - point_quantities[1])
+                else:
+                    momentum = 0
+                    vel = 0
+                
+            points_list.append(vel)
+            
+        if quantity == 'bearing':
+            points_list.append(calc_bearing(point_quantities[2],
+                                            point_quantities[3]))
+        if quantity == 'xcentroid':
+            points_list.append(centroids[point_i][0])
+
+        if quantity == 'ycentroid':
+            points_list.append(centroids[point_i][1])
+
+    return points_list
+    
+    
 ##
 # @brief Take gauge readings, given a gauge file and a sww mesh
@@ -163 +218 @@
     heading.insert(0,'time')
     heading.insert(1,'hours')
-
-    #create a list of csv writers for all the points and write header
-    points_writer = []
-    for point_i,point in enumerate(points):
-        points_writer.append(writer(file(dir_name + sep + gauge_file
-                                         + point_name[point_i] + '.csv', "wb")))
-        points_writer[point_i].writerow(heading)
     
     if verbose: log.critical('Writing csv files')
@@ -187 +235 @@
             quake_offset_time = callable_sww.starttime
 
+    for point_i, point in enumerate(points_array):
+        is_opened = False       
         for time in callable_sww.get_time():
-            for point_i, point in enumerate(points_array):
-                #add domain starttime to relative time.
-                quake_time = time + quake_offset_time
-                points_list = [quake_time, quake_time/3600.]# fudge around SWW time bug
-                point_quantities = callable_sww(time, point_i) # __call__ is overridden
-                            
-                for quantity in quantities:
-                    if quantity == NAN:
-                        log.critical('quantity does not exist in %s'
-                                     % callable_sww.get_name)
-                    else:
-                        #core quantities that are exported from the interpolator     
-                        if quantity == 'stage':
-                            points_list.append(point_quantities[0])
-                            
-                        if quantity == 'elevation':
-                            points_list.append(point_quantities[1])
-                            
-                        if quantity == 'xmomentum':
-                            points_list.append(point_quantities[2])
-                            
-                        if quantity == 'ymomentum':
-                            points_list.append(point_quantities[3])
-
-                        #derived quantities that are calculated from the core ones
-                        if quantity == 'depth':
-                            points_list.append(point_quantities[0] 
-                                               - point_quantities[1])
-
-                        if quantity == 'momentum':
-                            momentum = sqrt(point_quantities[2]**2 
-                                            + point_quantities[3]**2)
-                            points_list.append(momentum)
-                            
-                        if quantity == 'speed':
-                            #if depth is less than 0.001 then speed = 0.0
-                            if point_quantities[0] - point_quantities[1] < 0.001:
-                                vel = 0.0
-                            else:
-                                if point_quantities[2] < 1.0e6:
-                                    momentum = sqrt(point_quantities[2]**2
-                                                    + point_quantities[3]**2)
-                                    vel = momentum / (point_quantities[0] 
-                                                      - point_quantities[1])
-                                else:
-                                    momentum = 0
-                                    vel = 0
-                                
-                            points_list.append(vel)
-                            
-                        if quantity == 'bearing':
-                            points_list.append(calc_bearing(point_quantities[2],
-                                                            point_quantities[3]))
-                        if quantity == 'xcentroid':
-                            points_list.append(callable_sww.centroids[point_i][0])
-
-                        if quantity == 'ycentroid':
-                            points_list.append(callable_sww.centroids[point_i][1])
-                            
-                points_writer[point_i].writerow(points_list)
-
+            #add domain starttime to relative time.
+            quake_time = time + quake_offset_time
+            point_quantities = callable_sww(time, point_i) # __call__ is overridden
+
+            if point_quantities[0] != NAN:
+                if is_opened == False:
+                    points_writer = writer(file(dir_name + sep + gauge_file
+                                                        + point_name[point_i] + '.csv', "wb"))
+                    points_writer.writerow(heading)
+                    is_opened = True
+                points_list = [quake_time, quake_time/3600.] +  _quantities2csv(quantities, point_quantities, callable_sww.centroids, point_i)
+                points_writer.writerow(points_list)
+            else:
+                msg = 'gauge' + point_name[point_i] + 'falls off the mesh in file ' + sww_file + '.'
+                log.warning(msg)
 ##
 # @brief Read a .sww file and plot the time series.

anuga_core/source/anuga/abstract_2d_finite_volumes/test_gauge.py

@@ -301 +301 @@
 
        
-    # def test_sww2csv_gauge_point_off_mesh(self):
-        # from anuga.pmesh.mesh import Mesh
-        # from anuga.shallow_water import Domain, Transmissive_boundary
-        # from csv import reader,writer
-        # import time
-        # import string
-        
-        # """Most of this test was copied from test_interpolate
-        # test_interpole_sww2csv
-        
-        # This is testing the sww2csv_gauges function with one gauge off the mesh, by creating a sww file and
-        # then exporting the gauges and checking the results.
-        
-        # This tests the correct values for when a gauge is off the mesh, which is important for parallel.
-        # """
-        
-        # domain = self.domain
-        # self._create_sww()     
+    def test_sww2csv_gauge_point_off_mesh(self):
+        from anuga.pmesh.mesh import Mesh
+        from anuga.shallow_water import Domain, Transmissive_boundary
+        from csv import reader,writer
+        import time
+        import string
+        
+        """Most of this test was copied from test_interpolate
+        test_interpole_sww2csv
+        
+        This is testing the sww2csv_gauges function with one gauge off the mesh, by creating a sww file and
+        then exporting the gauges and checking the results.
+        
+        This tests the correct values for when a gauge is off the mesh, which is important for parallel.
+        """
+
+        domain = self.domain
+        sww = self._create_sww()
      
-        # # test the function
-        # points = [[50.0,1.],[50.5,-20.25]]
-
-# #        points_file = tempfile.mktemp(".csv")
-        # points_file = 'test_point.csv'
-        # file_id = open(points_file,"w")
-        # file_id.write("name,easting,northing \n\
-# point1, 50.0, 1.0\n\
-# point2, 50.5, 20.25\n")
-        # file_id.close()
-
-        # sww2csv_gauges(sww.filename, 
-                            # points_file,
-                            # quantities=['stage', 'elevation'],
-                            # use_cache=False,
-                            # verbose=False)
-
-        # point1_answers_array = [[0.0,1.0,-5.0], [2.0,10.0,-5.0]]
-        # point1_filename = 'gauge_point1.csv'
-        # point1_handle = file(point1_filename)
-        # point1_reader = reader(point1_handle)
-        # point1_reader.next()
-
-        # line=[]
-        # for i,row in enumerate(point1_reader):
-# #            print 'i',i,'row',row
-            # # note the 'hole' (element 1) below - skip the new 'hours' field
-            # line.append([float(row[0]),float(row[2]),float(row[3])])
-            # #print 'line',line[i],'point1',point1_answers_array[i]
-            # assert num.allclose(line[i], point1_answers_array[i])
-
-        # point2_answers_array = [[0.0,1.0,-0.5], [2.0,10.0,-0.5]]
-        # point2_filename = 'gauge_point2.csv' 
-        # point2_handle = file(point2_filename)
-        # point2_reader = reader(point2_handle)
-        # point2_reader.next()
-                        
-        # line=[]
-        # for i,row in enumerate(point2_reader):
-# #            print 'i',i,'row',row
-            # # note the 'hole' (element 1) below - skip the new 'hours' field
-            # line.append([float(row[0]),float(row[2]),float(row[3])])
-# #            print 'line',line[i],'point1',point1_answers_array[i]
-            # assert num.allclose(line[i], point2_answers_array[i])
-                         
-        # # clean up
-        # point1_handle.close()
-        # point2_handle.close()
-        # os.remove(points_file)
-# #        os.remove(point1_filename)
-# #        os.remove(point2_filename)
+        # test the function
+        points = [[50.0,1.],[50.5,-20.25]]
+
+#        points_file = tempfile.mktemp(".csv")
+        points_file = 'test_point.csv'
+        file_id = open(points_file,"w")
+        file_id.write("name,easting,northing \n\
+offmesh1, 50.0, 1.0\n\
+offmesh2, 50.5, 20.25\n")
+        file_id.close()
+
+        points_files = ['offmesh1.csv', 'offmesh2.csv']        
+        
+        for point_filename in points_files:
+            if os.path.exists(point_filename): os.remove(point_filename)         
+        
+        sww2csv_gauges(self.sww.filename, 
+                            points_file,
+                            quantities=['stage', 'elevation', 'bearing'],
+                            use_cache=False,
+                            verbose=False)
+
+        for point_filename in points_files: 
+            assert not os.path.exists(point_filename)
+            
+        os.remove(points_file)
         
         

anuga_core/source/anuga/abstract_2d_finite_volumes/util.py

@@ -14 +14 @@
 from shutil import copy
 
-from anuga.utilities.numerical_tools import ensure_numeric
+from anuga.utilities.numerical_tools import ensure_numeric, angle
 
 from math import sqrt, atan, degrees
@@ -300 +300 @@
     # * moving the reference direction from [1,0] to North
     # * changing from counter clockwise to clocwise.
-        
-    angle = degrees(atan(vh/(uh+1.e-15)))
-
-    if (0 < angle < 90.0):
-        if vh > 0:
-            bearing = 90.0 - abs(angle)
-        if vh < 0:
-            bearing = 270.0 - abs(angle)
-    
-    if (-90 < angle < 0):
-        if vh < 0:
-            bearing = 90.0 - (angle)
-        if vh > 0:
-            bearing = 270.0 - (angle)
-    if angle == 0: bearing = 0.0
-
-    return bearing
+    
+    return degrees(angle([uh, vh], [0, -1]))