Changeset 3702


Ignore:
Timestamp:
Oct 5, 2006, 5:32:43 PM (17 years ago)
Author:
ole
Message:

Added clipping method to Geospatial data.

Location:
anuga_core/source/anuga
Files:
3 edited

Legend:

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

    r3514 r3702  
    55
    66from os import access, F_OK, R_OK
    7 from Numeric import concatenate, array, Float, shape, reshape, ravel
    8 
     7from Numeric import concatenate, array, Float, shape, reshape, ravel, take
    98
    109from anuga.utilities.numerical_tools import ensure_numeric
     
    137136    def __len__(self):
    138137        return len(self.data_points)
     138
     139    def __repr__(self):
     140        return str(self.get_data_points(absolute=True))
     141   
    139142   
    140143    def check_data_points(self, data_points):
     
    207210        else:
    208211            verbose = False
     212
     213    def clip(self, polygon, closed=True):
     214        """Clip geospatial date by a polygon
     215
     216        Input
     217          polygon - Either a list of points, an Nx2 array or
     218                    a Geospatial data object.
     219          closed - (optional) determine whether points on boundary should be
     220          regarded as belonging to the polygon (closed = True)
     221          or not (closed = False). Default is True.
     222         
     223        Output
     224          Geospatial data object representing point inside specified polygon
     225       
     226        """
     227
     228        from anuga.utilities.polygon import inside_polygon
     229
     230        if isinstance(polygon, Geospatial_data):
     231            # Polygon is an object - extract points
     232            polygon = polygon.get_data_points()
     233
     234        points = self.get_data_points()   
     235        inside_indices = inside_polygon(points, polygon, closed)
     236        return Geospatial_data(take(points, inside_indices))
     237       
     238       
    209239   
    210240    def _set_using_lat_long(self,
  • anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

    r3663 r3702  
    2222
    2323    def test_0(self):
    24         """Basic points
     24        """test_0(self):
     25       
     26        Basic points
    2527        """
    2628        from anuga.coordinate_transforms.geo_reference import Geo_reference
     
    3133        assert allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]])
    3234
     35        #Check __repr__
     36
     37        rep = `G`
     38        ref = '[[ 1.   2.1]\n [ 3.   5.3]]'
     39        assert rep == ref
     40
     41
     42        #Check getter
     43        assert allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3]])
     44       
    3345        #Check defaults
    3446        assert G.attributes is None
     
    358370        assert allclose(P, concatenate( (points1,points2) ))
    359371                           
     372       
     373    def test_clip0(self):
     374        """test_clip0(self):
     375       
     376        Test that point sets can be clipped by a polygon
     377        """
     378       
     379        from anuga.coordinate_transforms.geo_reference import Geo_reference
     380       
     381        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
     382                  [0, 0], [2.4, 3.3]]
     383        G = Geospatial_data(points)
     384
     385        # First try the unit square   
     386        U = [[0,0], [1,0], [1,1], [0,1]]
     387        assert allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]])
     388
     389        # Then a more complex polygon
     390        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
     391        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
     392        G = Geospatial_data(points)
     393
     394        assert allclose(G.clip(polygon).get_data_points(),
     395                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
     396
     397    def test_clip1(self):
     398        """test_clip1(self):
     399       
     400        Test that point sets can be clipped by a polygon given as
     401        another Geospatial dataset
     402        """
     403       
     404        from anuga.coordinate_transforms.geo_reference import Geo_reference
     405       
     406        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
     407                  [0, 0], [2.4, 3.3]]
     408        G = Geospatial_data(points)
     409
     410        # First try the unit square   
     411        U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]])
     412        assert allclose(G.clip(U).get_data_points(),
     413                        [[0.2, 0.5], [0.4, 0.3], [0, 0]])
     414
     415        # Then a more complex polygon
     416        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
     417        G = Geospatial_data(points)
     418        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
     419       
     420
     421        assert allclose(G.clip(polygon).get_data_points(),
     422                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
     423
     424
    360425       
    361426
  • anuga_core/source/anuga/utilities/polygon.py

    r3633 r3702  
    8080        raise NameError, e
    8181    except:
    82         msg = 'Points could not be converted to Numeric array'
     82        # FIXME(Ole): This message is wrong. Shouldn't it be "Absolute"
     83        # rather than "Numeric"?
     84        msg = 'Points could not be converted to Numeric array'
    8385        raise msg
    8486
     
    8890        raise NameError, e
    8991    except:
     92        # FIXME(Ole): This message is wrong. Shouldn't it be "Absolute"
     93        # rather than "Numeric"?
    9094        msg = 'Polygon %s could not be converted to Numeric array' %(str(polygon))
    9195        raise msg
Note: See TracChangeset for help on using the changeset viewer.