Changeset 2594


Ignore:
Timestamp:
Mar 24, 2006, 6:07:58 PM (18 years ago)
Author:
ole
Message:

Duncan and Ole refactoring and dealing with georeferencing

Location:
inundation
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • inundation/coordinate_transforms/geo_reference.py

    r2593 r2594  
    77
    88import types, sys
    9 from Numeric import array,Float,ArrayType
     9from Numeric import array,Float,ArrayType,reshape
     10from utilities.numerical_tools import ensure_numeric
    1011
    1112DEFAULT_ZONE = sys.maxint
     
    131132        If the points do not have a geo ref, assume 'absolute' values
    132133        """
     134
     135        # FIXME: Use ensure_numeric
    133136       
    134137        is_list = False
    135138        if type(points) == types.ListType:
    136139            is_list = True
    137             if len(points)>0 and type(points[0]) not in [types.ListType,types.TupleType]:
    138                 #a single point is being passed.  make it a list of lists
    139                 points = [points]
    140         elif type(points) == ArrayType:
    141             if len(points.shape) == 1:
    142                 points = [points]
    143        
    144         # FIXME: currently untested
    145         # Always returns a list of lists
    146         if points_geo_ref is self:
    147             return points       
    148            
    149         #convert into array
    150         points = array(points).astype(Float)
    151         #add point geo ref to points
    152         if not points_geo_ref is None:
    153             points[:,0] += points_geo_ref.xllcorner
    154             points[:,1] += points_geo_ref.yllcorner
    155        
    156         #subtract primary geo ref from points
    157         points[:,0] -= self.xllcorner
    158         points[:,1] -= self.yllcorner
     140
     141        points = ensure_numeric(points)
     142        if len(points.shape) == 1:
     143            #One point has been passed
     144            msg = 'Single point must have two elements'
     145            assert len(points) == 2, msg
     146            points = reshape(points, (1,2))
     147
     148           
     149        if points_geo_ref is not self:
     150            #add point geo ref to points
     151            if not points_geo_ref is None:
     152                points[:,0] += points_geo_ref.xllcorner
     153                points[:,1] += points_geo_ref.yllcorner
     154       
     155            #subtract primary geo ref from points
     156            points[:,0] -= self.xllcorner
     157            points[:,1] -= self.yllcorner
     158
     159           
    159160        if is_list:
    160161            points = points.tolist()
     162           
    161163        return points
     164   
    162165   
    163166    def get_absolute(self, points):
     
    188191        return points
    189192
     193
     194    def reconcile_zones(self, other):
     195
     196        if self.zone == other.zone or self.zone == DEFAULT_ZONE and other.zone == DEFAULT_ZONE:
     197            pass
     198        elif self.zone == DEFAULT_ZONE:
     199            self.zone = other.zone
     200        elif other.zone == DEFAULT_ZONE:
     201            other.zone = self.zone           
     202        else:   
     203            msg = 'Both geospatial_data objects must be in the same \
     204            ZONE to allow reconciliation. I got zone %d and %d' %(self.zone, other.zone)
     205            raise msg
    190206   
    191207    #def easting_northing2geo_reffed_point(self, x, y):
  • inundation/coordinate_transforms/test_geo_reference.py

    r2261 r2594  
    159159        g = Geo_reference(56,x,y)
    160160        lofl = array([[3.0,323.0]])
    161         new_lofl = g.change_points_geo_ref(lofl)
     161
     162       
     163        #print "5 lofl before",lofl         
     164        new_lofl = g.change_points_geo_ref(lofl.copy())
    162165        #print "5 lofl",lofl
    163166        #print "5 new_lofl",new_lofl
     
    165168        self.failUnless(type(new_lofl) == ArrayType, ' failed')
    166169        self.failUnless(type(new_lofl) == type(lofl), ' failed')
     170
     171
    167172        for point,new_point in map(None,lofl,new_lofl):
    168173            self.failUnless(point[0]-x==new_point[0], ' failed')
     
    174179        g = Geo_reference(56,x,y)
    175180        lofl = array([355.0,3.0])
    176         new_lofl = g.change_points_geo_ref(lofl)
     181        new_lofl = g.change_points_geo_ref(lofl.copy())       
    177182        #print "lofl",lofl
    178183        #print "new_lofl",new_lofl
     
    223228     
    224229        self.failUnless(g == new_g, 'test___cmp__ failed')   
     230
     231
     232    def test_reconcile(self):
     233        g1 = Geo_reference(56,2,5)
     234        g2 = Geo_reference(50,4,5)
     235        g3 = Geo_reference(50,66,6)       
     236        g_default = Geo_reference()               
     237
     238
     239        g2.reconcile_zones(g3)
     240        assert g2.get_zone() == g3.get_zone()
     241
     242        g_default.reconcile_zones(g3)
     243        assert g_default.get_zone() == g3.get_zone()
     244
     245        g_default = Geo_reference()               
     246        g3.reconcile_zones(g_default)
     247        assert g_default.get_zone() == g3.get_zone()               
     248
     249        try:
     250            g1.reconcile_zones(g2)
     251        except:
     252            pass
     253        else:
     254            msg = 'Should have raised an exception'
     255            raise msg
     256       
     257       
     258       
    225259       
    226260       
  • inundation/geospatial_data/geospatial_data.py

    r2592 r2594  
    176176        Default: absolute is True.
    177177        """
     178
    178179        if absolute is True:
    179             xll = self.geo_reference.xllcorner
    180             yll = self.geo_reference.yllcorner
    181             return self.data_points + [xll, yll]
     180            return self.geo_reference.get_absolute(self.data_points)
    182181        else:
    183             return self.data_points
     182            return self.data_points       
     183       
    184184   
    185185    def get_attributes(self, attribute_name = None):
     
    202202        return self.attributes[attribute_name]
    203203
     204
    204205    def __add__(self, other):
     206        """
     207        returns the addition of 2 geospatical objects,
     208        objects are concatencated to the end of each other
     209           
     210        NOTE: doesn't add if objects contain different
     211        attributes 
     212        """
     213
     214
     215        # find objects zone and checks if the same
     216        geo_ref1 = self.get_geo_reference()
     217        zone1 = geo_ref1.get_zone()
     218       
     219        geo_ref2 = other.get_geo_reference()
     220        zone2 = geo_ref2.get_zone()
     221
     222        geo_ref1.reconcile_zones(geo_ref2)
     223
     224
     225        # sets xll and yll as the smallest from self and other
     226        # FIXME (Duncan and Ole): use lower left corner derived from absolute coordinates
     227        if self.geo_reference.xllcorner <= other.geo_reference.xllcorner:
     228            xll = self.geo_reference.xllcorner
     229        else:
     230            xll = other.geo_reference.xllcorner
     231
     232        if self.geo_reference.yllcorner <= other.geo_reference.yllcorner:
     233            yll = self.geo_reference.yllcorner
     234        else:
     235            yll = other.geo_reference.yllcorner
     236        new_geo_ref = Geo_reference(geo_ref1.get_zone(), xll, yll)
     237
     238        xll = yll = 0.
     239        relative_points1 = self.get_data_points(absolute=False)
     240        relative_points2 = other.get_data_points(absolute=False)
     241           
     242        #print 'R1', relative_points1
     243        #print 'R2', relative_points2       
     244
     245       
     246        new_relative_points1 = new_geo_ref.change_points_geo_ref(relative_points1, geo_ref1)
     247        new_relative_points2 = new_geo_ref.change_points_geo_ref(relative_points2, geo_ref2)
     248
     249        #print 'NR1', new_relative_points1
     250        #print 'NR2', new_relative_points2       
     251       
     252        #Now both point sets are relative to new_geo_ref and zones have been reconciled
     253
     254
     255        # Concatenate points
     256        new_points = concatenate((new_relative_points1,
     257                                  new_relative_points2),
     258                                 axis = 0)
     259
     260     
     261        # Concatenate attributes
     262        new_attributes = {}
     263        for x in self.attributes.keys():
     264            if other.attributes.has_key(x):
     265
     266                attrib1 = self.attributes[x]
     267                attrib2 = other.attributes[x]
     268                new_attributes[x] = concatenate((attrib1, attrib2))
     269
     270            else:
     271                msg = 'Both geospatial_data objects must have the same \n'
     272                msg += 'attributes to allow addition.'
     273                raise msg
     274
     275       
     276        # Instantiate new data object and return   
     277        return Geospatial_data(new_points,
     278                               new_attributes,
     279                               new_geo_ref)
     280
     281
     282
     283    def xxx__add__(self, other):
    205284        """
    206285        returns the addition of 2 geospatical objects,
  • inundation/geospatial_data/test_geospatial_data.py

    r2591 r2594  
    234234
    235235        P = G.get_data_points(absolute=True)
     236
     237        P_relative = G.get_data_points(absolute=False)
     238       
     239        assert allclose(P_relative, P - [0.1, 2.0])
     240
    236241        assert allclose(P, concatenate( (P1,P2) ))
    237 
    238         P_relative = G.get_data_points(absolute=False)
    239         assert allclose(P_relative, P - [0.1, 2.0])                     
    240242        assert allclose(P, [[2.0, 4.1], [4.0, 7.3],
    241243                            [5.1, 9.1], [6.1, 6.3]])
    242 
     244       
     245
     246
     247       
     248
     249    def test_add_with_geo_absolute (self):
     250        """
     251        Difference in Geo_reference resolved
     252        """
     253        points1 = array([[2.0, 4.1], [4.0, 7.3]])
     254        points2 = array([[5.1, 9.1], [6.1, 6.3]])       
     255        attributes1 = [2, 4]
     256        attributes2 = [5, 76]
     257        geo_ref1= Geo_reference(55, 1.0, 2.0)
     258        geo_ref2 = Geo_reference(55, 2.0, 3.0)
     259
     260       
     261                               
     262        G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()],
     263                             attributes1, geo_ref1)
     264       
     265        G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()],
     266                             attributes2, geo_ref2)
     267
     268        #Check that absolute values are as expected
     269        P1 = G1.get_data_points(absolute=True)
     270        assert allclose(P1, points1)
     271
     272        P1 = G1.get_data_points(absolute=False)
     273        assert allclose(P1, points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()])       
     274
     275        P2 = G2.get_data_points(absolute=True)
     276        assert allclose(P2, points2)
     277
     278        P2 = G2.get_data_points(absolute=False)
     279        assert allclose(P2, points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()])               
     280       
     281        G = G1 + G2
     282       
     283        assert allclose(G.get_geo_reference().get_xllcorner(), 1.0)
     284        assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)
     285
     286        P = G.get_data_points(absolute=True)
     287
     288        P_relative = G.get_data_points(absolute=False)
     289       
     290        assert allclose(P_relative, [[1.0, 2.1], [3.0, 5.3], [4.1, 7.1], [5.1, 4.3]])
     291
     292        assert allclose(P, concatenate( (points1,points2) ))
     293                           
     294       
    243295
    244296
Note: See TracChangeset for help on using the changeset viewer.