Changeset 3714


Ignore:
Timestamp:
Oct 9, 2006, 3:53:01 PM (17 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:
2 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/geospatial_data/geospatial_data.py

    r3712 r3714  
    238238       
    239239       
    240    
     240    def clip_outside(self, polygon, closed=True):
     241        """Clip geospatial date by a polygon, keeping data OUTSIDE of polygon
     242
     243        Input
     244          polygon - Either a list of points, an Nx2 array or
     245                    a Geospatial data object.
     246          closed - (optional) determine whether points on boundary should be
     247          regarded as belonging to the polygon (closed = True)
     248          or not (closed = False). Default is True.
     249         
     250        Output
     251          Geospatial data object representing point OUTSIDE specified polygon
     252       
     253        """
     254
     255        from anuga.utilities.polygon import outside_polygon
     256
     257        if isinstance(polygon, Geospatial_data):
     258            # Polygon is an object - extract points
     259            polygon = polygon.get_data_points()
     260
     261        points = self.get_data_points()   
     262        outside_indices = outside_polygon(points, polygon, closed)
     263        return Geospatial_data(take(points, outside_indices))
     264   
    241265    def _set_using_lat_long(self,
    242266                            latitudes,
  • anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

    r3702 r3714  
    423423
    424424
    425        
     425    def test_clip0_outside(self):
     426        """test_clip0(self):
     427       
     428        Test that point sets can be clipped outside of a polygon
     429        """
     430       
     431        from anuga.coordinate_transforms.geo_reference import Geo_reference
     432       
     433        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
     434                  [0, 0], [2.4, 3.3]]
     435        G = Geospatial_data(points)
     436
     437        # First try the unit square   
     438        U = [[0,0], [1,0], [1,1], [0,1]]
     439        assert allclose(G.clip_outside(U).get_data_points(),
     440                        [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
     441
     442        # Then a more complex polygon
     443        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
     444        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
     445        G = Geospatial_data(points)
     446
     447        assert allclose(G.clip_outside(polygon).get_data_points(),
     448                        [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])   
     449
     450
     451    def test_clip1_outside(self):
     452            """test_clip1(self):
     453           
     454            Test that point sets can be clipped outside of a polygon given as
     455            another Geospatial dataset
     456            """
     457           
     458            from anuga.coordinate_transforms.geo_reference import Geo_reference
     459           
     460            points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
     461                      [0, 0], [2.4, 3.3]]
     462            G = Geospatial_data(points)
     463
     464            # First try the unit square   
     465            U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
     466            assert allclose(G.clip_outside(U).get_data_points(),
     467                            [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
     468
     469            # Then a more complex polygon
     470            points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
     471            G = Geospatial_data(points)
     472            polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
     473           
     474
     475            assert allclose(G.clip_outside(polygon).get_data_points(),
     476                            [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
     477
    426478
    427479
Note: See TracChangeset for help on using the changeset viewer.