Changeset 8541

Timestamp:

Aug 30, 2012, 10:35:14 AM (13 years ago)

Author:

steve

Message:

Dealing with holes

Location:

trunk/anuga_core/source/anuga

Files:

• file_conversion/sww2dem.py (modified) (2 diffs)
• fit_interpolate/interpolate.py (modified) (2 diffs)
• fit_interpolate/test_interpolate.py (modified) (1 diff)

trunk/anuga_core/source/anuga/file_conversion/sww2dem.py

@@ -330 +330 @@
     if verbose: log.critical('Interpolating')
     grid_values = interp.interpolate(result, grid_points).flatten()
+    outside_indices = interp.get_outside_poly_indices()
 
     if verbose:
@@ -336 +337 @@
 
     # 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)
+    
+#    P = interp.mesh.get_boundary_polygon()
+#    outside_indices = outside_polygon(grid_points, P, closed=True)
 
     for i in outside_indices:

trunk/anuga_core/source/anuga/fit_interpolate/interpolate.py

@@ -371 +371 @@
 
 
+    def get_outside_poly_indices(self):
+        """
+        Return index of those data points outside (and in holes)
+        the mesh
+
+        Precondition: interpolation or interpolation_block has been called
+        """
+        return self.outside_poly_indices
+
+
     def _get_point_data_z(self, f, verbose=False):
         """
@@ -413 +423 @@
         if verbose: log.critical('Getting indices inside mesh boundary')
 
-        # Quick test against boundary, but will not deal with holes in the mesh 
+        # Quick test against boundary, but will not deal with holes in the mesh,
+        # that is done below
         inside_boundary_indices, outside_poly_indices = \
             in_and_outside_polygon(point_coordinates,

trunk/anuga_core/source/anuga/fit_interpolate/test_interpolate.py

@@ -1715 +1715 @@
             self.failUnless( z[i,0] == answer[11,0], 'Fail!')
 
+
+
+
+    def test_points_in_hole(self):
+
+        v0  = [0.0,     0.0]
+        v1  = [1.0/3.0, 0.0]
+        v2  = [2.0/3.0, 0.0]
+        v3  = [1.0,     0.0]
+        v4  = [0.0,     1.0/3.0]
+        v5  = [1.0/3.0, 1.0/3.0]
+        v6  = [2.0/3.0, 1.0/3.0]
+        v7  = [1.0,     1.0/3.0]
+        v8  = [0.0,     2.0/3.0]
+        v9  = [1.0/3.0, 2.0/3.0]
+        v10 = [2.0/3.0, 2.0/3.0]
+        v11 = [1.0,     2.0/3.0]
+        v12 = [0.0,     1.0]
+        v13 = [1.0/3.0, 1.0]
+        v14 = [2.0/3.0, 1.0]
+        v15 = [1.0,     1.0]
+
+
+        vertices = [v0, v1, v2, v3, v4, v5, v6, v7, v8, v9, v10, v11, v12, v13, v14, v15]
+
+        triangles = [[1,4,0], [1,5,4],    [2,5,1], [2,6,5],       [3,6,2],  [3,7,6],
+                     [5,8,4], [5,9,8],                            [7,10,6], [7,11,10],
+                     [9,12,8], [9,13,12], [10,13,9], [10,14,13], [11,14,10], [11,15,14]]
+
+
+
+
+        point_coords = [[0.25, 0.25],
+                        [0.25, 0.2],
+                        [10.0, 12.0], # point Outside poly
+                        [0.5, 0.5], # point in hole
+                        [0.5, 0.4], # point in hole
+                        [0.0, 0.0],
+                        [1.0, 0.0]]
+
+        geo_data = Geospatial_data(data_points = point_coords)
+
+        interp = Interpolate(vertices, triangles)
+        f = num.array([linear_function(vertices),2*linear_function(vertices)])
+        f = num.transpose(f)
+        #print "f",f
+        z = interp.interpolate(f, geo_data)
+        #z = interp.interpolate(f, point_coords)
+        answer = [linear_function(point_coords),
+                  2*linear_function(point_coords) ]
+        answer = num.transpose(answer)
+        answer[2,:] = [NAN, NAN]
+        answer[3,:] = [NAN, NAN]
+        answer[4,:] = [NAN, NAN]
+        #print "z",z
+        #print "answer _ fixed",answer
+        assert num.allclose(z[0:1], answer[0:1])
+        assert num.allclose(z[5:6], answer[5:6])
+        for i in [2,3,4]:
+            self.failUnless( z[i,1] == answer[2,1], 'Fail!')
+            self.failUnless( z[i,0] == answer[2,0], 'Fail!')
+
+
     def test_interpolate_sww2csv(self):