Changeset 3714

Timestamp:

Oct 9, 2006, 3:53:01 PM (18 years ago)

Author:

sexton

Message:

incorporate clip_outside function to grab data which doesn't fall in a given polygon

Location:

anuga_core/source/anuga/geospatial_data

Files:

• geospatial_data.py (modified) (1 diff)
• test_geospatial_data.py (modified) (1 diff)

anuga_core/source/anuga/geospatial_data/geospatial_data.py

@@ -238 +238 @@
         
         
-   
+    def clip_outside(self, polygon, closed=True):
+        """Clip geospatial date by a polygon, keeping data OUTSIDE of polygon
+
+        Input
+          polygon - Either a list of points, an Nx2 array or
+                    a Geospatial data object.
+          closed - (optional) determine whether points on boundary should be
+          regarded as belonging to the polygon (closed = True)
+          or not (closed = False). Default is True.
+          
+        Output
+          Geospatial data object representing point OUTSIDE specified polygon
+        
+        """
+
+        from anuga.utilities.polygon import outside_polygon
+
+        if isinstance(polygon, Geospatial_data):
+            # Polygon is an object - extract points
+            polygon = polygon.get_data_points()
+
+        points = self.get_data_points()    
+        outside_indices = outside_polygon(points, polygon, closed)
+        return Geospatial_data(take(points, outside_indices))
+    
     def _set_using_lat_long(self,
                             latitudes,

anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

@@ -423 +423 @@
 
 
-        
+    def test_clip0_outside(self):
+        """test_clip0(self):
+        
+        Test that point sets can be clipped outside of a polygon
+        """
+        
+        from anuga.coordinate_transforms.geo_reference import Geo_reference
+        
+        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
+                  [0, 0], [2.4, 3.3]]
+        G = Geospatial_data(points)
+
+        # First try the unit square    
+        U = [[0,0], [1,0], [1,1], [0,1]]
+        assert allclose(G.clip_outside(U).get_data_points(),
+                        [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
+
+        # Then a more complex polygon
+        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
+        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+        G = Geospatial_data(points)
+
+        assert allclose(G.clip_outside(polygon).get_data_points(),
+                        [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])    
+
+
+    def test_clip1_outside(self):
+            """test_clip1(self):
+            
+            Test that point sets can be clipped outside of a polygon given as
+            another Geospatial dataset
+            """
+            
+            from anuga.coordinate_transforms.geo_reference import Geo_reference
+            
+            points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
+                      [0, 0], [2.4, 3.3]]
+            G = Geospatial_data(points)
+
+            # First try the unit square    
+            U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 
+            assert allclose(G.clip_outside(U).get_data_points(),
+                            [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
+
+            # Then a more complex polygon
+            points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
+            G = Geospatial_data(points)
+            polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
+            
+
+            assert allclose(G.clip_outside(polygon).get_data_points(),
+                            [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
+