Changeset 3738


Ignore:
Timestamp:
Oct 11, 2006, 10:22:45 AM (17 years ago)
Author:
ole
Message:

Added attributes into clip functions and tested

Location:
anuga_core/source/anuga/geospatial_data
Files:
2 edited

Legend:

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

    r3735 r3738  
    55
    66from os import access, F_OK, R_OK
     7from types import DictType
     8
    79from Numeric import concatenate, array, Float, shape, reshape, ravel, take
    810
     
    124126
    125127        else:
    126 ###                data = {}   
    127                 # watch for case where file name and points, attributes etc are provided!!
    128                 # if file name then all provided info will be removed!
    129                 self.import_points_file(file_name, delimiter, verbose)
     128            # watch for case where file name and points, attributes etc are provided!!
     129            # if file name then all provided info will be removed!
     130            self.import_points_file(file_name, delimiter, verbose)
    130131               
    131                 self.check_data_points(self.data_points)
    132                 self.set_attributes(self.attributes)
    133                 self.set_geo_reference(self.geo_reference)
    134                 self.set_default_attribute_name(default_attribute_name)
     132            self.check_data_points(self.data_points)
     133            self.set_attributes(self.attributes)
     134            self.set_geo_reference(self.geo_reference)
     135            self.set_default_attribute_name(default_attribute_name)
     136
     137
     138        assert self.attributes is None or isinstance(self.attributes, DictType)
     139       
    135140
    136141    def __len__(self):
     
    159164        """
    160165       
    161         from types import DictType
    162        
    163166        if attributes is None:
    164167            self.attributes = None
     
    235238        points = self.get_data_points()   
    236239        inside_indices = inside_polygon(points, polygon, closed)
    237         return Geospatial_data(take(points, inside_indices))
     240
     241        clipped_points = take(points, inside_indices)
     242
     243        # Clip all attributes
     244        attributes = self.get_all_attributes()
     245
     246        clipped_attributes = {}
     247        if attributes is not None:
     248            for key, att in attributes.items():
     249                clipped_attributes[key] = take(att, inside_indices)
     250
     251        return Geospatial_data(clipped_points, clipped_attributes)
    238252       
    239253       
     
    261275        points = self.get_data_points()   
    262276        outside_indices = outside_polygon(points, polygon, closed)
    263         return Geospatial_data(take(points, outside_indices))
     277
     278        clipped_points = take(points, outside_indices)
     279
     280        # Clip all attributes
     281        attributes = self.get_all_attributes()
     282
     283        clipped_attributes = {}
     284        if attributes is not None:
     285            for key, att in attributes.items():
     286                clipped_attributes[key] = take(att, outside_indices)
     287
     288        return Geospatial_data(clipped_points, clipped_attributes)
     289
    264290   
    265291    def _set_using_lat_long(self,
     
    350376        """
    351377        Return values for all attributes.
     378        The return value is either None or a dictionary (possibly empty).
    352379        """
    353380
     
    445472        Note: will throw an IOError if it can't load the file.
    446473        Catch these!
     474
     475        Post condition: self.attributes dictionary has been set
    447476        """
    448477       
  • anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

    r3736 r3738  
    395395                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
    396396
    397     def no_test_clip1(self):
     397    def test_clip0_with_attributes(self):
     398        """test_clip0_with_attributes(self):
     399       
     400        Test that point sets with attributes can be clipped by a polygon
     401        """
     402       
     403        from anuga.coordinate_transforms.geo_reference import Geo_reference
     404       
     405        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
     406                  [0, 0], [2.4, 3.3]]
     407
     408        attributes = [2, -4, 5, 76, -2, 0.1, 3]
     409        att_dict = {'att1': attributes,
     410                    'att2': array(attributes)+1}
     411       
     412        G = Geospatial_data(points, att_dict)
     413
     414        # First try the unit square   
     415        U = [[0,0], [1,0], [1,1], [0,1]]
     416        assert allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]])
     417        assert allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
     418        assert allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])               
     419
     420        # Then a more complex polygon
     421        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
     422        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
     423
     424        # This time just one attribute
     425        attributes = [2, -4, 5, 76, -2, 0.1]
     426        G = Geospatial_data(points, attributes)
     427
     428        assert allclose(G.clip(polygon).get_data_points(),
     429                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
     430        assert allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
     431       
     432
     433    def test_clip1(self):
    398434        """test_clip1(self):
    399435       
     
    408444        attributes = [2, -4, 5, 76, -2, 0.1, 3]
    409445        att_dict = {'att1': attributes,
    410                     'att2': array(attributes) +1}
     446                    'att2': array(attributes)+1}
    411447        G = Geospatial_data(points, att_dict)
    412448       
     
    417453
    418454        assert allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1])
    419         #assert allclose(G.clip(U).get_attributes('att2'), ???)
     455        assert allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1])                       
    420456       
    421457        # Then a more complex polygon
    422458        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
    423         G = Geospatial_data(points)
     459        attributes = [2, -4, 5, 76, -2, 0.1]       
     460        G = Geospatial_data(points, attributes)
    424461        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
    425462       
     
    427464        assert allclose(G.clip(polygon).get_data_points(),
    428465                        [[0.5, 0.5], [1, -0.5], [1.5, 0]])
     466        assert allclose(G.clip(polygon).get_attributes(), [-4, 5, 76])
     467       
    429468
    430469    def test_clip0_outside(self):
     
    438477        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
    439478                  [0, 0], [2.4, 3.3]]
    440         G = Geospatial_data(points)
     479        attributes = [2, -4, 5, 76, -2, 0.1, 3]       
     480        G = Geospatial_data(points, attributes)
    441481
    442482        # First try the unit square   
     
    444484        assert allclose(G.clip_outside(U).get_data_points(),
    445485                        [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
     486        #print G.clip_outside(U).get_attributes()
     487        assert allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])       
     488       
    446489
    447490        # Then a more complex polygon
    448491        polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]
    449492        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
    450         G = Geospatial_data(points)
     493        attributes = [2, -4, 5, 76, -2, 0.1]       
     494        G = Geospatial_data(points, attributes)
    451495
    452496        assert allclose(G.clip_outside(polygon).get_data_points(),
    453                         [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])   
     497                        [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
     498        assert allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])               
    454499
    455500
     
    465510        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
    466511                  [0, 0], [2.4, 3.3]]
    467         G = Geospatial_data(points)
     512        attributes = [2, -4, 5, 76, -2, 0.1, 3]       
     513        G = Geospatial_data(points, attributes)       
    468514
    469515        # First try the unit square   
     
    471517        assert allclose(G.clip_outside(U).get_data_points(),
    472518                        [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]])
     519        assert allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])       
    473520
    474521        # Then a more complex polygon
    475522        points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]]
    476         G = Geospatial_data(points)
     523        attributes = [2, -4, 5, 76, -2, 0.1]       
     524        G = Geospatial_data(points, attributes)
     525
    477526        polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])
    478527       
     
    480529        assert allclose(G.clip_outside(polygon).get_data_points(),
    481530                        [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]])
     531        assert allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])
     532       
    482533
    483534
     
    493544        points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3],
    494545                  [0, 0], [2.4, 3.3]]
    495 
    496         G = Geospatial_data(points)
     546        attributes = [2, -4, 5, 76, -2, 0.1, 3]       
     547        G = Geospatial_data(points, attributes)
    497548
    498549        # First try the unit square   
     
    500551        G1 = G.clip(U)
    501552        assert allclose(G1.get_data_points(),[[0.2, 0.5], [0.4, 0.3], [0, 0]])
     553        assert allclose(G.clip(U).get_attributes(), [-4, 76, 0.1])
    502554       
    503555        G2 = G.clip_outside(U)
    504556        assert allclose(G2.get_data_points(),[[-1, 4], [1.0, 2.1],
    505557                                              [3.0, 5.3], [2.4, 3.3]])
     558        assert allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3])               
    506559
    507560       
     
    510563                     [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]
    511564
     565        new_attributes = [-4, 76, 0.1, 2, 5, -2, 3]                 
     566       
    512567        assert allclose((G1+G2).get_data_points(), new_points)
     568        assert allclose((G1+G2).get_attributes(), new_attributes)
    513569
    514570        G = G1+G2
     
    523579        # Check result
    524580        assert allclose(G3.get_data_points(), new_points)       
    525        
     581        assert allclose(G3.get_attributes(), new_attributes)       
    526582       
    527583        os.remove(FN)
Note: See TracChangeset for help on using the changeset viewer.