Changeset 2594
- Timestamp:
- Mar 24, 2006, 6:07:58 PM (19 years ago)
- Location:
- inundation
- Files:
-
- 4 edited
Legend:
- Unmodified
- Added
- Removed
-
inundation/coordinate_transforms/geo_reference.py
r2593 r2594 7 7 8 8 import types, sys 9 from Numeric import array,Float,ArrayType 9 from Numeric import array,Float,ArrayType,reshape 10 from utilities.numerical_tools import ensure_numeric 10 11 11 12 DEFAULT_ZONE = sys.maxint … … 131 132 If the points do not have a geo ref, assume 'absolute' values 132 133 """ 134 135 # FIXME: Use ensure_numeric 133 136 134 137 is_list = False 135 138 if type(points) == types.ListType: 136 139 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 159 160 if is_list: 160 161 points = points.tolist() 162 161 163 return points 164 162 165 163 166 def get_absolute(self, points): … … 188 191 return points 189 192 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 190 206 191 207 #def easting_northing2geo_reffed_point(self, x, y): -
inundation/coordinate_transforms/test_geo_reference.py
r2261 r2594 159 159 g = Geo_reference(56,x,y) 160 160 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()) 162 165 #print "5 lofl",lofl 163 166 #print "5 new_lofl",new_lofl … … 165 168 self.failUnless(type(new_lofl) == ArrayType, ' failed') 166 169 self.failUnless(type(new_lofl) == type(lofl), ' failed') 170 171 167 172 for point,new_point in map(None,lofl,new_lofl): 168 173 self.failUnless(point[0]-x==new_point[0], ' failed') … … 174 179 g = Geo_reference(56,x,y) 175 180 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()) 177 182 #print "lofl",lofl 178 183 #print "new_lofl",new_lofl … … 223 228 224 229 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 225 259 226 260 -
inundation/geospatial_data/geospatial_data.py
r2592 r2594 176 176 Default: absolute is True. 177 177 """ 178 178 179 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) 182 181 else: 183 return self.data_points 182 return self.data_points 183 184 184 185 185 def get_attributes(self, attribute_name = None): … … 202 202 return self.attributes[attribute_name] 203 203 204 204 205 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): 205 284 """ 206 285 returns the addition of 2 geospatical objects, -
inundation/geospatial_data/test_geospatial_data.py
r2591 r2594 234 234 235 235 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 236 241 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])240 242 assert allclose(P, [[2.0, 4.1], [4.0, 7.3], 241 243 [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 243 295 244 296
Note: See TracChangeset
for help on using the changeset viewer.