Changeset 2062

Timestamp:

Nov 24, 2005, 1:05:20 PM (19 years ago)

Author:

ole

Message:

Found, unit tested, and fixed bug in outside polygon.
Then refactored sww2dem, to a more efficient way of assigning missing values

Location:

inundation

Files:

• pyvolution/data_manager.py (modified) (2 diffs)
• pyvolution/test_data_manager.py (modified) (1 diff)
• utilities/polygon.py (modified) (1 diff)
• utilities/test_polygon.py (modified) (1 diff)

inundation/pyvolution/data_manager.py

@@ -1480 +1480 @@
     from Numeric import array2string
     
-    from utilities.polygon import inside_polygon
+    from utilities.polygon import inside_polygon, outside_polygon, separate_points_by_polygon
     from util import apply_expression_to_dictionary
     
@@ -1676 +1676 @@
 
 
-    #Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
-    P = interp.mesh.get_boundary_polygon()
-    inside_indices = inside_polygon(grid_points, P)
-
-    #FIXME: Remove when new interpolation takes care of this
-    for i in range(nrows):
-        if verbose and i%((nrows+10)/10)==0:
-            print 'Doing row %d of %d' %(i, nrows)
-
-        base_index = (nrows-i-1)*ncols     
-        for j in range(ncols):
-            index = base_index+j
-
-            if not sometrue(inside_indices == index):
-                grid_values[index] = NODATA_value
-
-    #############            
-
-    
-
     if verbose:
         print 'Interpolated values are in [%f, %f]' %(min(grid_values),
                                                       max(grid_values))
+
+    #Assign NODATA_value to all points outside bounding polygon (from interpolation mesh)
+    P = interp.mesh.get_boundary_polygon()
+    outside_indices = outside_polygon(grid_points, P, closed=True)
+
+    for i in outside_indices:
+        grid_values[i] = NODATA_value        
+
+
+
 
     if format.lower() == 'ers':

inundation/pyvolution/test_data_manager.py

@@ -3590 +3590 @@
     suite = unittest.makeSuite(Test_Data_Manager,'test')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_missing_points')
-    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_boundingbox')    
+    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2dem_asc_elevation')    
     #suite = unittest.makeSuite(Test_Data_Manager,'test_dem2pts_bounding_box')
     #suite = unittest.makeSuite(Test_Data_Manager,'test_decimate_dem')

inundation/utilities/polygon.py

@@ -116 +116 @@
         return count != 1
     else:
-        return indices[count:][::-1]  #return reversed
+        if count == len(indices):
+            # No points are outside
+            return []
+        else:    
+            return indices[count:][::-1]  #return reversed
 
 

inundation/utilities/test_polygon.py

@@ -247 +247 @@
         #evaluate to True as the point 0.5, 1.0 is outside the unit square
 
+    def test_all_outside_polygon(self):
+        """Test case where all points are outside poly
+        """
+        
+        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square    
+
+        points = [[2,2], [1,3], [-1,1], [0,2]] #All outside
+
+
+        indices, count = separate_points_by_polygon(points, U)
+        #print indices, count
+        assert count == 0 #None inside
+        assert allclose(indices, [3,2,1,0])
+
+        indices = outside_polygon(points, U, closed = True)
+        assert allclose(indices, [0,1,2,3])
+
+        indices = inside_polygon(points, U, closed = True)
+        assert allclose(indices, [])                
+
+
+    def test_all_inside_polygon(self):
+        """Test case where all points are inside poly
+        """
+        
+        U = [[0,0], [1,0], [1,1], [0,1]] #Unit square    
+
+        points = [[0.5,0.5], [0.2,0.3], [0,0.5]] #All inside (or on edge)
+
+
+        indices, count = separate_points_by_polygon(points, U)
+        assert count == 3 #All inside
+        assert allclose(indices, [0,1,2])
+
+        indices = outside_polygon(points, U, closed = True)
+        assert allclose(indices, [])
+
+        indices = inside_polygon(points, U, closed = True)
+        assert allclose(indices, [0,1,2])
+        
+
     def test_separate_points_by_polygon(self):
         U = [[0,0], [1,0], [1,1], [0,1]] #Unit square