Changeset 3702

Timestamp:

Oct 5, 2006, 5:32:43 PM (18 years ago)

Author:

ole

Message:

Added clipping method to Geospatial data.

Location:

anuga_core/source/anuga

Files:

• geospatial_data/geospatial_data.py (modified) (3 diffs)
• geospatial_data/test_geospatial_data.py (modified) (3 diffs)
• utilities/polygon.py (modified) (2 diffs)

anuga_core/source/anuga/geospatial_data/geospatial_data.py

@@ -5 +5 @@
 
 from os import access, F_OK, R_OK
-from Numeric import concatenate, array, Float, shape, reshape, ravel
-
+from Numeric import concatenate, array, Float, shape, reshape, ravel, take
 
 from anuga.utilities.numerical_tools import ensure_numeric
@@ -137 +136 @@
     def __len__(self):
         return len(self.data_points)
+
+    def __repr__(self):
+        return str(self.get_data_points(absolute=True))
+    
     
     def check_data_points(self, data_points):
@@ -207 +210 @@
         else:
             verbose = False
+
+    def clip(self, polygon, closed=True):
+        """Clip geospatial date by a 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 inside specified polygon
+        
+        """
+
+        from anuga.utilities.polygon import inside_polygon
+
+        if isinstance(polygon, Geospatial_data):
+            # Polygon is an object - extract points
+            polygon = polygon.get_data_points()
+
+        points = self.get_data_points()    
+        inside_indices = inside_polygon(points, polygon, closed)
+        return Geospatial_data(take(points, inside_indices))
+        
+        
    
     def _set_using_lat_long(self,

anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

@@ -22 +22 @@
 
     def test_0(self):
-        """Basic points
+        """test_0(self):
+        
+        Basic points
         """
         from anuga.coordinate_transforms.geo_reference import Geo_reference
@@ -31 +33 @@
         assert allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]])
 
+        #Check __repr__
+
+        rep = `G`
+        ref = '[[ 1.   2.1]\n [ 3.   5.3]]'
+        assert rep == ref
+
+
+        #Check getter
+        assert allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3]])
+        
         #Check defaults
         assert G.attributes is None
@@ -358 +370 @@
         assert allclose(P, concatenate( (points1,points2) ))
                            
+        
+    def test_clip0(self):
+        """test_clip0(self):
+        
+        Test that point sets can be clipped by 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(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]])
+
+        # 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(polygon).get_data_points(),
+                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
+
+    def test_clip1(self):
+        """test_clip1(self):
+        
+        Test that point sets can be clipped by 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(U).get_data_points(),
+                        [[0.2, 0.5], [0.4, 0.3], [0, 0]])
+
+        # 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(polygon).get_data_points(),
+                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
+
+
         
 

anuga_core/source/anuga/utilities/polygon.py

@@ -80 +80 @@
         raise NameError, e
     except:
-        msg = 'Points could not be converted to Numeric array'
+        # FIXME(Ole): This message is wrong. Shouldn't it be "Absolute"
+        # rather than "Numeric"?
+        msg = 'Points could not be converted to Numeric array' 
         raise msg
 
@@ -88 +90 @@
         raise NameError, e
     except:
+        # FIXME(Ole): This message is wrong. Shouldn't it be "Absolute"
+        # rather than "Numeric"?
         msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
         raise msg