- Timestamp:
- Feb 10, 2009, 11:11:04 AM (15 years ago)
- Location:
- branches/numpy
- Files:
-
- 1 edited
- 1 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/numpy/anuga/geospatial_data/test_geospatial_data.py
r6153 r6304 1 1 #!/usr/bin/env python 2 3 2 4 3 import unittest 5 4 import os 5 import tempfile 6 6 from math import sqrt, pi 7 import tempfile8 7 from sets import ImmutableSet 9 8 10 import Numericas num9 import numpy as num 11 10 12 11 from anuga.geospatial_data.geospatial_data import * … … 16 15 from anuga.utilities.system_tools import get_host_name 17 16 from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a 17 from anuga.config import netcdf_float 18 18 19 19 20 class Test_Geospatial_data(unittest.TestCase): … … 24 25 pass 25 26 26 27 27 def test_0(self): 28 28 #Basic points 29 29 from anuga.coordinate_transforms.geo_reference import Geo_reference 30 30 31 31 points = [[1.0, 2.1], [3.0, 5.3]] 32 32 G = Geospatial_data(points) 33 34 33 assert num.allclose(G.data_points, [[1.0, 2.1], [3.0, 5.3]]) 35 34 … … 38 37 rep = `G` 39 38 ref = '[[ 1. 2.1]\n [ 3. 5.3]]' 40 41 msg = 'Representation %s is not equal to %s' %(rep, ref) 39 msg = 'Representation %s is not equal to %s' % (rep, ref) 42 40 assert rep == ref, msg 43 41 44 42 #Check getter 45 43 assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3]]) 46 44 47 45 #Check defaults 48 46 assert G.attributes is None 49 50 47 assert G.geo_reference.zone == Geo_reference().zone 51 48 assert G.geo_reference.xllcorner == Geo_reference().xllcorner 52 49 assert G.geo_reference.yllcorner == Geo_reference().yllcorner 53 54 50 55 51 def test_1(self): 56 52 points = [[1.0, 2.1], [3.0, 5.3]] 57 53 attributes = [2, 4] 58 G = Geospatial_data(points, attributes) 54 G = Geospatial_data(points, attributes) 59 55 assert G.attributes.keys()[0] == DEFAULT_ATTRIBUTE 60 56 assert num.allclose(G.attributes.values()[0], [2, 4]) 61 62 57 63 58 def test_2(self): 64 59 from anuga.coordinate_transforms.geo_reference import Geo_reference 60 65 61 points = [[1.0, 2.1], [3.0, 5.3]] 66 62 attributes = [2, 4] … … 72 68 assert G.geo_reference.yllcorner == 200 73 69 74 75 70 def test_get_attributes_1(self): 76 71 from anuga.coordinate_transforms.geo_reference import Geo_reference 72 77 73 points = [[1.0, 2.1], [3.0, 5.3]] 78 74 attributes = [2, 4] … … 80 76 geo_reference=Geo_reference(56, 100, 200)) 81 77 82 83 78 P = G.get_data_points(absolute=False) 84 assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]]) 79 assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]]) 85 80 86 81 P = G.get_data_points(absolute=True) 87 assert num.allclose(P, [[101.0, 202.1], [103.0, 205.3]]) 82 assert num.allclose(P, [[101.0, 202.1], [103.0, 205.3]]) 88 83 89 84 V = G.get_attributes() #Simply get them … … 95 90 def test_get_attributes_2(self): 96 91 #Multiple attributes 97 98 99 92 from anuga.coordinate_transforms.geo_reference import Geo_reference 93 100 94 points = [[1.0, 2.1], [3.0, 5.3]] 101 95 attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]} … … 104 98 default_attribute_name='a1') 105 99 106 107 100 P = G.get_data_points(absolute=False) 108 assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]]) 109 101 assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]]) 102 110 103 V = G.get_attributes() #Get default attribute 111 104 assert num.allclose(V, [2, 4]) … … 125 118 pass 126 119 else: 127 raise 'Should have raised exception'120 raise Exception, 'Should have raised exception' 128 121 129 122 def test_get_data_points(self): 130 points_ab = [[12.5, 34.7],[-4.5,-60.0]]123 points_ab = [[12.5, 34.7], [-4.5, -60.0]] 131 124 x_p = -10 132 125 y_p = -40 133 126 geo_ref = Geo_reference(56, x_p, y_p) 134 127 points_rel = geo_ref.change_points_geo_ref(points_ab) 135 128 136 129 spatial = Geospatial_data(points_rel, geo_reference=geo_ref) 137 138 130 results = spatial.get_data_points(absolute=False) 139 140 131 assert num.allclose(results, points_rel) 141 132 142 133 x_p = -1770 143 134 y_p = 4.01 144 135 geo_ref = Geo_reference(56, x_p, y_p) 145 136 points_rel = geo_ref.change_points_geo_ref(points_ab) 146 results = spatial.get_data_points \ 147 ( geo_reference=geo_ref) 148 137 results = spatial.get_data_points(geo_reference=geo_ref) 138 149 139 assert num.allclose(results, points_rel) 150 140 151 152 141 def test_get_data_points_lat_long(self): 153 # lat long [-30.],[130] 154 #Zone: 52 155 #Easting: 596450.153 Northing: 6680793.777 156 # lat long [-32.],[131] 157 #Zone: 52 158 #Easting: 688927.638 Northing: 6457816.509 159 160 points_Lat_long = [[-30.,130], [-32,131]] 161 162 spatial = Geospatial_data(latitudes=[-30, -32.], 163 longitudes=[130, 131]) 164 142 # lat long [-30.], [130] 143 # Zone: 52 144 # Easting: 596450.153 Northing: 6680793.777 145 # lat long [-32.], [131] 146 # Zone: 52 147 # Easting: 688927.638 Northing: 6457816.509 148 149 points_Lat_long = [[-30., 130], [-32, 131]] 150 151 spatial = Geospatial_data(latitudes=[-30, -32.], longitudes=[130, 131]) 165 152 results = spatial.get_data_points(as_lat_long=True) 166 #print "test_get_data_points_lat_long - results", results167 #print "points_Lat_long",points_Lat_long168 153 assert num.allclose(results, points_Lat_long) 169 154 170 155 def test_get_data_points_lat_longII(self): 171 156 # x,y North,east long,lat 172 157 boundary_polygon = [[ 250000, 7630000]] 173 158 zone = 50 174 159 175 160 geo_reference = Geo_reference(zone=zone) 176 geo = Geospatial_data(boundary_polygon ,geo_reference=geo_reference)161 geo = Geospatial_data(boundary_polygon ,geo_reference=geo_reference) 177 162 seg_lat_long = geo.get_data_points(as_lat_long=True) 178 lat_result = degminsec2decimal_degrees(-21,24,54) 179 long_result = degminsec2decimal_degrees(114,35,17.89) 180 #print "seg_lat_long", seg_lat_long [0][0] 181 #print "lat_result",lat_result 163 lat_result = degminsec2decimal_degrees(-21, 24, 54) 164 long_result = degminsec2decimal_degrees(114, 35, 17.89) 165 assert num.allclose(seg_lat_long[0][0], lat_result) #lat 166 assert num.allclose(seg_lat_long[0][1], long_result) #long 167 168 def test_get_data_points_lat_longIII(self): 169 # x,y North,east long,lat 170 # for northern hemisphere 171 boundary_polygon = [[419944.8, 918642.4]] 172 zone = 47 173 174 geo_reference = Geo_reference(zone=zone) 175 geo = Geospatial_data(boundary_polygon, geo_reference=geo_reference) 176 seg_lat_long = geo.get_data_points(as_lat_long=True, 177 isSouthHemisphere=False) 178 lat_result = degminsec2decimal_degrees(8.31, 0, 0) 179 long_result = degminsec2decimal_degrees(98.273, 0, 0) 182 180 assert num.allclose(seg_lat_long[0][0], lat_result)#lat 183 181 assert num.allclose(seg_lat_long[0][1], long_result)#long 184 182 185 186 def test_get_data_points_lat_longIII(self):187 # x,y North,east long,lat188 #for northern hemisphere189 boundary_polygon = [[419944.8, 918642.4]]190 zone = 47191 192 geo_reference = Geo_reference(zone=zone)193 geo = Geospatial_data(boundary_polygon,194 geo_reference=geo_reference)195 196 seg_lat_long = geo.get_data_points(as_lat_long=True,197 isSouthHemisphere=False)198 199 lat_result = degminsec2decimal_degrees(8.31,0,0)200 long_result = degminsec2decimal_degrees(98.273,0,0)201 #print "seg_lat_long", seg_lat_long [0]202 #print "lat_result",lat_result203 assert num.allclose(seg_lat_long[0][0], lat_result)#lat204 assert num.allclose(seg_lat_long[0][1], long_result)#long205 206 207 208 183 def test_set_geo_reference(self): 209 """test_set_georeference210 211 Test that georeference can be changed without changing the 184 '''test_set_georeference 185 186 Test that georeference can be changed without changing the 212 187 absolute values. 213 """214 215 points_ab = [[12.5, 34.7],[-4.5,-60.0]]188 ''' 189 190 points_ab = [[12.5, 34.7], [-4.5, -60.0]] 216 191 x_p = -10 217 192 y_p = -40 218 193 geo_ref = Geo_reference(56, x_p, y_p) 219 194 points_rel = geo_ref.change_points_geo_ref(points_ab) 220 195 221 196 # Create without geo_ref properly set 222 G = Geospatial_data(points_rel) 197 G = Geospatial_data(points_rel) 223 198 assert not num.allclose(points_ab, G.get_data_points(absolute=True)) 224 199 225 200 # Create the way it should be 226 201 G = Geospatial_data(points_rel, geo_reference=geo_ref) 227 202 assert num.allclose(points_ab, G.get_data_points(absolute=True)) 228 203 229 204 # Change georeference and check that absolute values are unchanged. 230 205 x_p = 10 … … 232 207 new_geo_ref = Geo_reference(56, x_p, y_p) 233 208 G.set_geo_reference(new_geo_ref) 234 assert num.allclose(points_ab, G.get_data_points(absolute=True)) 235 236 237 238 209 msg = ('points_ab=%s, but\nG.get_data_points(absolute=True)=%s' 210 % (str(points_ab), str(G.get_data_points(absolute=True)))) 211 assert num.allclose(points_ab, G.get_data_points(absolute=True)), msg 212 239 213 def test_conversions_to_points_dict(self): 240 214 #test conversions to points_dict 241 242 243 215 from anuga.coordinate_transforms.geo_reference import Geo_reference 216 244 217 points = [[1.0, 2.1], [3.0, 5.3]] 245 218 attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]} … … 248 221 default_attribute_name='a1') 249 222 250 251 223 points_dict = geospatial_data2points_dictionary(G) 252 224 253 225 assert points_dict.has_key('pointlist') 254 assert points_dict.has_key('attributelist') 226 assert points_dict.has_key('attributelist') 255 227 assert points_dict.has_key('geo_reference') 256 228 … … 260 232 assert A.has_key('a0') 261 233 assert A.has_key('a1') 262 assert A.has_key('a2') 234 assert A.has_key('a2') 263 235 264 236 assert num.allclose( A['a0'], [0, 0] ) 265 assert num.allclose( A['a1'], [2, 4] ) 237 assert num.allclose( A['a1'], [2, 4] ) 266 238 assert num.allclose( A['a2'], [79.4, -7] ) 267 268 239 269 240 geo = points_dict['geo_reference'] 270 241 assert geo is G.geo_reference 271 242 272 273 243 def test_conversions_from_points_dict(self): 274 """test conversions from points_dict 275 """ 244 '''test conversions from points_dict''' 276 245 277 246 from anuga.coordinate_transforms.geo_reference import Geo_reference 278 247 279 248 points = [[1.0, 2.1], [3.0, 5.3]] 280 249 attributes = {'a0': [0, 0], 'a1': [2, 4], 'a2': [79.4, -7]} … … 284 253 points_dict['attributelist'] = attributes 285 254 points_dict['geo_reference'] = Geo_reference(56, 100, 200) 286 287 255 288 256 G = points_dictionary2geospatial_data(points_dict) 289 290 257 P = G.get_data_points(absolute=False) 291 assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]]) 292 293 #V = G.get_attribute_values() #Get default attribute 294 #assert allclose(V, [2, 4]) 258 assert num.allclose(P, [[1.0, 2.1], [3.0, 5.3]]) 295 259 296 260 V = G.get_attributes('a0') #Get by name … … 304 268 305 269 def test_add(self): 306 """ test the addition of two geospatical objects 307 no geo_reference see next test 308 """ 270 '''test the addition of two geospatical objects 271 no geo_reference see next test 272 ''' 273 309 274 points = [[1.0, 2.1], [3.0, 5.3]] 310 attributes = {'depth':[2, 4], 'elevation':[6.1, 5]} 311 attributes1 = {'depth':[2, 4], 'elevation':[2.5, 1]} 312 G1 = Geospatial_data(points, attributes) 313 G2 = Geospatial_data(points, attributes1) 314 315 # g3 = geospatial_data2points_dictionary(G1) 316 # print 'g3=', g3 317 275 attributes = {'depth': [2, 4], 'elevation': [6.1, 5]} 276 attributes1 = {'depth': [2, 4], 'elevation': [2.5, 1]} 277 G1 = Geospatial_data(points, attributes) 278 G2 = Geospatial_data(points, attributes1) 279 318 280 G = G1 + G2 319 281 … … 326 288 327 289 def test_addII(self): 328 """ test the addition of two geospatical objects 329 no geo_reference see next test 330 """ 290 '''test the addition of two geospatical objects 291 no geo_reference see next test 292 ''' 293 331 294 points = [[1.0, 2.1], [3.0, 5.3]] 332 attributes = {'depth': [2, 4]}333 G1 = Geospatial_data(points, attributes) 334 295 attributes = {'depth': [2, 4]} 296 G1 = Geospatial_data(points, attributes) 297 335 298 points = [[5.0, 2.1], [3.0, 50.3]] 336 attributes = {'depth': [200, 400]}299 attributes = {'depth': [200, 400]} 337 300 G2 = Geospatial_data(points, attributes) 338 339 # g3 = geospatial_data2points_dictionary(G1) 340 # print 'g3=', g3 341 301 342 302 G = G1 + G2 343 303 344 assert G.attributes.has_key('depth') 304 assert G.attributes.has_key('depth') 345 305 assert G.attributes.keys(), ['depth'] 346 306 assert num.allclose(G.attributes['depth'], [2, 4, 200, 400]) 347 307 assert num.allclose(G.get_data_points(), [[1.0, 2.1], [3.0, 5.3], 348 308 [5.0, 2.1], [3.0, 50.3]]) 309 349 310 def test_add_with_geo (self): 350 """ 351 Difference in Geo_reference resolved 352 """ 311 '''Difference in Geo_reference resolved''' 312 353 313 points1 = [[1.0, 2.1], [3.0, 5.3]] 354 314 points2 = [[5.0, 6.1], [6.0, 3.3]] … … 362 322 projection='UTM', 363 323 units='m') 364 324 365 325 G1 = Geospatial_data(points1, attributes1, geo_ref1) 366 326 G2 = Geospatial_data(points2, attributes2, geo_ref2) … … 371 331 372 332 P2 = G2.get_data_points(absolute=True) 373 assert num.allclose(P2, [[5.1, 9.1], [6.1, 6.3]]) 374 333 assert num.allclose(P2, [[5.1, 9.1], [6.1, 6.3]]) 334 375 335 G = G1 + G2 376 336 … … 381 341 P = G.get_data_points(absolute=True) 382 342 383 #P_relative = G.get_data_points(absolute=False)384 #385 #assert allclose(P_relative, P - [0.1, 2.0])386 387 343 assert num.allclose(P, num.concatenate( (P1,P2) )) 344 msg = 'P=%s' % str(P) 388 345 assert num.allclose(P, [[2.0, 4.1], [4.0, 7.3], 389 [5.1, 9.1], [6.1, 6.3]]) 390 391 392 393 346 [5.1, 9.1], [6.1, 6.3]]), msg 394 347 395 348 def test_add_with_geo_absolute (self): 396 """ 397 Difference in Geo_reference resolved 398 """ 349 '''Difference in Geo_reference resolved''' 350 399 351 points1 = num.array([[2.0, 4.1], [4.0, 7.3]]) 400 points2 = num.array([[5.1, 9.1], [6.1, 6.3]]) 352 points2 = num.array([[5.1, 9.1], [6.1, 6.3]]) 401 353 attributes1 = [2, 4] 402 354 attributes2 = [5, 76] … … 404 356 geo_ref2 = Geo_reference(55, 2.0, 3.0) 405 357 406 407 408 G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()], 358 G1 = Geospatial_data(points1 - [geo_ref1.get_xllcorner(), 359 geo_ref1.get_yllcorner()], 409 360 attributes1, geo_ref1) 410 411 G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()], 361 362 G2 = Geospatial_data(points2 - [geo_ref2.get_xllcorner(), 363 geo_ref2.get_yllcorner()], 412 364 attributes2, geo_ref2) 413 365 … … 417 369 418 370 P1 = G1.get_data_points(absolute=False) 419 assert num.allclose(P1, points1 - [geo_ref1.get_xllcorner(), geo_ref1.get_yllcorner()]) 371 msg = ('P1=%s, but\npoints1 - <...>=%s' 372 % (str(P1), 373 str(points1 - [geo_ref1.get_xllcorner(), 374 geo_ref1.get_yllcorner()]))) 375 assert num.allclose(P1, points1 - [geo_ref1.get_xllcorner(), 376 geo_ref1.get_yllcorner()]), msg 420 377 421 378 P2 = G2.get_data_points(absolute=True) … … 423 380 424 381 P2 = G2.get_data_points(absolute=False) 425 assert num.allclose(P2, points2 - [geo_ref2.get_xllcorner(), geo_ref2.get_yllcorner()]) 426 382 assert num.allclose(P2, points2 - [geo_ref2.get_xllcorner(), 383 geo_ref2.get_yllcorner()]) 384 427 385 G = G1 + G2 428 429 #assert allclose(G.get_geo_reference().get_xllcorner(), 1.0)430 #assert allclose(G.get_geo_reference().get_yllcorner(), 2.0)431 432 386 P = G.get_data_points(absolute=True) 433 387 434 #P_relative = G.get_data_points(absolute=False) 435 # 436 #assert allclose(P_relative, [[1.0, 2.1], [3.0, 5.3], [4.1, 7.1], [5.1, 4.3]]) 437 438 assert num.allclose(P, num.concatenate( (points1,points2) )) 439 388 assert num.allclose(P, num.concatenate((points1, points2))) 440 389 441 390 def test_add_with_None(self): 442 """ test that None can be added to a geospatical objects 443 """ 444 391 '''test that None can be added to a geospatical objects''' 392 445 393 points1 = num.array([[2.0, 4.1], [4.0, 7.3]]) 446 points2 = num.array([[5.1, 9.1], [6.1, 6.3]]) 394 points2 = num.array([[5.1, 9.1], [6.1, 6.3]]) 447 395 448 396 geo_ref1= Geo_reference(55, 1.0, 2.0) … … 453 401 projection='UTM', 454 402 units='m') 455 456 457 attributes1 = {'depth':[2, 4.7], 'elevation':[6.1, 5]} 458 attributes2 = {'depth':[-2.3, 4], 'elevation':[2.5, 1]} 459 403 404 attributes1 = {'depth': [2, 4.7], 'elevation': [6.1, 5]} 405 attributes2 = {'depth': [-2.3, 4], 'elevation': [2.5, 1]} 460 406 461 407 G1 = Geospatial_data(points1, attributes1, geo_ref1) … … 465 411 assert G1.attributes.has_key('elevation') 466 412 assert num.allclose(G1.attributes['depth'], [2, 4.7]) 467 assert num.allclose(G1.attributes['elevation'], [6.1, 5]) 468 413 assert num.allclose(G1.attributes['elevation'], [6.1, 5]) 414 469 415 G2 = Geospatial_data(points2, attributes2, geo_ref2) 470 416 assert num.allclose(G2.get_geo_reference().get_xllcorner(), 0.1) … … 473 419 assert G2.attributes.has_key('elevation') 474 420 assert num.allclose(G2.attributes['depth'], [-2.3, 4]) 475 assert num.allclose(G2.attributes['elevation'], [2.5, 1]) 421 assert num.allclose(G2.attributes['elevation'], [2.5, 1]) 476 422 477 423 #Check that absolute values are as expected … … 480 426 481 427 P2 = G2.get_data_points(absolute=True) 482 assert num.allclose(P2, [[5.2, 12.1], [6.2, 9.3]]) 428 assert num.allclose(P2, [[5.2, 12.1], [6.2, 9.3]]) 483 429 484 430 # Normal add 485 431 G = G1 + None 486 487 432 assert G.attributes.has_key('depth') 488 433 assert G.attributes.has_key('elevation') 489 434 assert num.allclose(G.attributes['depth'], [2, 4.7]) 490 assert num.allclose(G.attributes['elevation'], [6.1, 5]) 435 assert num.allclose(G.attributes['elevation'], [6.1, 5]) 491 436 492 437 # Points are now absolute. 493 438 assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0) 494 439 assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0) 495 P = G.get_data_points(absolute=True) 496 assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]])497 440 P = G.get_data_points(absolute=True) 441 msg = 'P=%s' % str(P) 442 assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]]), msg 498 443 499 444 G = G2 + None … … 501 446 assert G.attributes.has_key('elevation') 502 447 assert num.allclose(G.attributes['depth'], [-2.3, 4]) 503 assert num.allclose(G.attributes['elevation'], [2.5, 1]) 504 448 assert num.allclose(G.attributes['elevation'], [2.5, 1]) 505 449 assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0) 506 450 assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0) 507 P = G.get_data_points(absolute=True) 451 452 P = G.get_data_points(absolute=True) 508 453 assert num.allclose(P, [[5.2, 12.1], [6.2, 9.3]]) 509 510 511 454 512 455 # Reverse add 513 456 G = None + G1 514 515 457 assert G.attributes.has_key('depth') 516 458 assert G.attributes.has_key('elevation') 517 459 assert num.allclose(G.attributes['depth'], [2, 4.7]) 518 assert num.allclose(G.attributes['elevation'], [6.1, 5]) 460 assert num.allclose(G.attributes['elevation'], [6.1, 5]) 519 461 520 462 # Points are now absolute. 521 463 assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0) 522 464 assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0) 523 P = G.get_data_points(absolute=True) 524 assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]])525 465 466 P = G.get_data_points(absolute=True) 467 assert num.allclose(P, [[3.0, 6.1], [5.0, 9.3]]) 526 468 527 469 G = None + G2 … … 529 471 assert G.attributes.has_key('elevation') 530 472 assert num.allclose(G.attributes['depth'], [-2.3, 4]) 531 assert num.allclose(G.attributes['elevation'], [2.5, 1]) 473 assert num.allclose(G.attributes['elevation'], [2.5, 1]) 532 474 533 475 assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0) 534 476 assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0) 535 P = G.get_data_points(absolute=True) 477 478 P = G.get_data_points(absolute=True) 536 479 assert num.allclose(P, [[5.2, 12.1], [6.2, 9.3]]) 537 480 538 539 540 541 542 543 481 def test_clip0(self): 544 """test_clip0(self):545 482 '''test_clip0(self): 483 546 484 Test that point sets can be clipped by a polygon 547 """548 485 ''' 486 549 487 from anuga.coordinate_transforms.geo_reference import Geo_reference 550 488 551 489 points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3], 552 490 [0, 0], [2.4, 3.3]] 553 491 G = Geospatial_data(points) 554 492 555 # First try the unit square 556 U = [[0,0], [1,0], [1,1], [0,1]] 557 assert num.allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]]) 493 # First try the unit square 494 U = [[0,0], [1,0], [1,1], [0,1]] 495 assert num.allclose(G.clip(U).get_data_points(), 496 [[0.2, 0.5], [0.4, 0.3], [0, 0]]) 558 497 559 498 # Then a more complex polygon 560 499 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] 561 points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] 500 points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], 501 [0.5, 1.5], [0.5, -0.5]] 562 502 G = Geospatial_data(points) 563 503 564 504 assert num.allclose(G.clip(polygon).get_data_points(), 565 [[0.5, 0.5], [1, -0.5], [1.5, 0]])505 [[0.5, 0.5], [1, -0.5], [1.5, 0]]) 566 506 567 507 def test_clip0_with_attributes(self): 568 """test_clip0_with_attributes(self):569 508 '''test_clip0_with_attributes(self): 509 570 510 Test that point sets with attributes can be clipped by a polygon 571 """572 511 ''' 512 573 513 from anuga.coordinate_transforms.geo_reference import Geo_reference 574 514 575 515 points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3], 576 516 [0, 0], [2.4, 3.3]] … … 579 519 att_dict = {'att1': attributes, 580 520 'att2': num.array(attributes)+1} 581 521 582 522 G = Geospatial_data(points, att_dict) 583 523 584 # First try the unit square 585 U = [[0,0], [1,0], [1,1], [0,1]] 586 assert num.allclose(G.clip(U).get_data_points(), [[0.2, 0.5], [0.4, 0.3], [0, 0]]) 524 # First try the unit square 525 U = [[0,0], [1,0], [1,1], [0,1]] 526 assert num.allclose(G.clip(U).get_data_points(), 527 [[0.2, 0.5], [0.4, 0.3], [0, 0]]) 587 528 assert num.allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1]) 588 assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1]) 529 assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1]) 589 530 590 531 # Then a more complex polygon 591 532 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] 592 points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] 533 points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], 534 [0.5, 1.5], [0.5, -0.5]] 593 535 594 536 # This time just one attribute 595 537 attributes = [2, -4, 5, 76, -2, 0.1] 596 538 G = Geospatial_data(points, attributes) 597 598 539 assert num.allclose(G.clip(polygon).get_data_points(), 599 [[0.5, 0.5], [1, -0.5], [1.5, 0]])540 [[0.5, 0.5], [1, -0.5], [1.5, 0]]) 600 541 assert num.allclose(G.clip(polygon).get_attributes(), [-4, 5, 76]) 601 602 542 603 543 def test_clip1(self): 604 """test_clip1(self):605 544 '''test_clip1(self): 545 606 546 Test that point sets can be clipped by a polygon given as 607 547 another Geospatial dataset 608 """609 548 ''' 549 610 550 from anuga.coordinate_transforms.geo_reference import Geo_reference 611 551 612 552 points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3], 613 553 [0, 0], [2.4, 3.3]] … … 616 556 'att2': num.array(attributes)+1} 617 557 G = Geospatial_data(points, att_dict) 618 619 # First try the unit square 620 U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 558 559 # First try the unit square 560 U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 621 561 assert num.allclose(G.clip(U).get_data_points(), 622 [[0.2, 0.5], [0.4, 0.3], [0, 0]]) 623 562 [[0.2, 0.5], [0.4, 0.3], [0, 0]]) 624 563 assert num.allclose(G.clip(U).get_attributes('att1'), [-4, 76, 0.1]) 625 assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1]) 626 564 assert num.allclose(G.clip(U).get_attributes('att2'), [-3, 77, 1.1]) 565 627 566 # Then a more complex polygon 628 points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] 629 attributes = [2, -4, 5, 76, -2, 0.1] 567 points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], 568 [0.5, 1.5], [0.5, -0.5]] 569 attributes = [2, -4, 5, 76, -2, 0.1] 630 570 G = Geospatial_data(points, attributes) 631 polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]])632 571 polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], 572 [2,1], [0,1]]) 633 573 634 574 assert num.allclose(G.clip(polygon).get_data_points(), 635 575 [[0.5, 0.5], [1, -0.5], [1.5, 0]]) 636 576 assert num.allclose(G.clip(polygon).get_attributes(), [-4, 5, 76]) 637 638 577 639 578 def test_clip0_outside(self): 640 """test_clip0_outside(self):641 579 '''test_clip0_outside(self): 580 642 581 Test that point sets can be clipped outside of a polygon 643 """644 582 ''' 583 645 584 from anuga.coordinate_transforms.geo_reference import Geo_reference 646 585 647 586 points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3], 648 587 [0, 0], [2.4, 3.3]] 649 attributes = [2, -4, 5, 76, -2, 0.1, 3] 588 attributes = [2, -4, 5, 76, -2, 0.1, 3] 650 589 G = Geospatial_data(points, attributes) 651 590 652 # First try the unit square 591 # First try the unit square 653 592 U = [[0,0], [1,0], [1,1], [0,1]] 654 593 assert num.allclose(G.clip_outside(U).get_data_points(), 655 594 [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]) 656 #print G.clip_outside(U).get_attributes() 657 assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3]) 658 595 assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3]) 659 596 660 597 # Then a more complex polygon 661 598 polygon = [[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]] 662 points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] 663 attributes = [2, -4, 5, 76, -2, 0.1] 599 points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], 600 [0.5, 1.5], [0.5, -0.5]] 601 attributes = [2, -4, 5, 76, -2, 0.1] 664 602 G = Geospatial_data(points, attributes) 665 666 603 assert num.allclose(G.clip_outside(polygon).get_data_points(), 667 604 [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]]) 668 assert num.allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1])669 605 assert num.allclose(G.clip_outside(polygon).get_attributes(), 606 [2, -2, 0.1]) 670 607 671 608 def test_clip1_outside(self): 672 """test_clip1_outside(self):673 609 '''test_clip1_outside(self): 610 674 611 Test that point sets can be clipped outside of a polygon given as 675 612 another Geospatial dataset 676 """677 613 ''' 614 678 615 from anuga.coordinate_transforms.geo_reference import Geo_reference 679 616 680 617 points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3], 681 618 [0, 0], [2.4, 3.3]] 682 attributes = [2, -4, 5, 76, -2, 0.1, 3] 683 G = Geospatial_data(points, attributes) 684 685 # First try the unit square 686 U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 619 attributes = [2, -4, 5, 76, -2, 0.1, 3] 620 G = Geospatial_data(points, attributes) 621 622 # First try the unit square 623 U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 687 624 assert num.allclose(G.clip_outside(U).get_data_points(), 688 625 [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]]) 689 assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1]) 626 assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1]) 690 627 691 628 # Then a more complex polygon 692 points = [ [0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], [0.5, 1.5], [0.5, -0.5]] 693 attributes = [2, -4, 5, 76, -2, 0.1] 629 points = [[0.5, 1.4], [0.5, 0.5], [1, -0.5], [1.5, 0], 630 [0.5, 1.5], [0.5, -0.5]] 631 attributes = [2, -4, 5, 76, -2, 0.1] 694 632 G = Geospatial_data(points, attributes) 695 696 polygon = Geospatial_data([[0,0], [1,0], [0.5,-1], [2, -1], [2,1], [0,1]]) 697 698 633 polygon = Geospatial_data([[0, 0], [1, 0], [0.5, -1], [2, -1], 634 [2, 1], [0, 1]]) 699 635 assert num.allclose(G.clip_outside(polygon).get_data_points(), 700 636 [[0.5, 1.4], [0.5, 1.5], [0.5, -0.5]]) 701 assert num.allclose(G.clip_outside(polygon).get_attributes(), [2, -2, 0.1]) 702 703 637 assert num.allclose(G.clip_outside(polygon).get_attributes(), 638 [2, -2, 0.1]) 704 639 705 640 def test_clip1_inside_outside(self): 706 """test_clip1_inside_outside(self):707 641 '''test_clip1_inside_outside(self): 642 708 643 Test that point sets can be clipped outside of a polygon given as 709 644 another Geospatial dataset 710 """711 645 ''' 646 712 647 from anuga.coordinate_transforms.geo_reference import Geo_reference 713 648 714 649 points = [[-1, 4], [0.2, 0.5], [1.0, 2.1], [0.4, 0.3], [3.0, 5.3], 715 650 [0, 0], [2.4, 3.3]] 716 attributes = [2, -4, 5, 76, -2, 0.1, 3] 651 attributes = [2, -4, 5, 76, -2, 0.1, 3] 717 652 G = Geospatial_data(points, attributes) 718 653 719 # First try the unit square 720 U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 654 # First try the unit square 655 U = Geospatial_data([[0,0], [1,0], [1,1], [0,1]]) 721 656 G1 = G.clip(U) 722 assert num.allclose(G1.get_data_points(),[[0.2, 0.5], [0.4, 0.3], [0, 0]]) 657 assert num.allclose(G1.get_data_points(), 658 [[0.2, 0.5], [0.4, 0.3], [0, 0]]) 723 659 assert num.allclose(G.clip(U).get_attributes(), [-4, 76, 0.1]) 724 725 660 G2 = G.clip_outside(U) 726 661 assert num.allclose(G2.get_data_points(),[[-1, 4], [1.0, 2.1], 727 662 [3.0, 5.3], [2.4, 3.3]]) 728 assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3]) 729 730 663 assert num.allclose(G.clip_outside(U).get_attributes(), [2, 5, -2, 3]) 664 731 665 # New ordering 732 new_points = [[0.2, 0.5], [0.4, 0.3], [0, 0]] +\ 733 [[-1, 4], [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]] 734 735 new_attributes = [-4, 76, 0.1, 2, 5, -2, 3] 736 666 new_points = [[0.2, 0.5], [0.4, 0.3], [0, 0], [-1, 4], 667 [1.0, 2.1], [3.0, 5.3], [2.4, 3.3]] 668 new_attributes = [-4, 76, 0.1, 2, 5, -2, 3] 669 737 670 assert num.allclose((G1+G2).get_data_points(), new_points) 738 671 assert num.allclose((G1+G2).get_attributes(), new_attributes) … … 742 675 G.export_points_file(FN) 743 676 744 745 677 # Read it back in 746 678 G3 = Geospatial_data(FN) 747 679 748 749 680 # Check result 750 assert num.allclose(G3.get_data_points(), new_points) 751 assert num.allclose(G3.get_attributes(), new_attributes) 752 681 assert num.allclose(G3.get_data_points(), new_points) 682 assert num.allclose(G3.get_attributes(), new_attributes) 683 753 684 os.remove(FN) 754 685 755 756 686 def test_load_csv(self): 757 758 import os 759 import tempfile 760 761 fileName = tempfile.mktemp(".csv") 762 file = open(fileName,"w") 763 file.write("x,y,elevation speed \n\ 687 fileName = tempfile.mktemp('.csv') 688 file = open(fileName,'w') 689 file.write('x,y,elevation speed \n\ 764 690 1.0 0.0 10.0 0.0\n\ 765 691 0.0 1.0 0.0 10.0\n\ 766 1.0 0.0 10.4 40.0\n ")692 1.0 0.0 10.4 40.0\n') 767 693 file.close() 768 #print fileName 694 769 695 results = Geospatial_data(fileName) 770 696 os.remove(fileName) 771 # print 'data', results.get_data_points() 772 assert num.allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0], 773 [1.0, 0.0]]) 697 assert num.allclose(results.get_data_points(), 698 [[1.0, 0.0],[0.0, 1.0], [1.0, 0.0]]) 774 699 assert num.allclose(results.get_attributes(attribute_name='elevation'), 775 700 [10.0, 0.0, 10.4]) … … 777 702 [0.0, 10.0, 40.0]) 778 703 779 780 ###################### .CSV ############################## 704 ################################################################################ 705 # Test CSV files 706 ################################################################################ 781 707 782 708 def test_load_csv_lat_long_bad_blocking(self): 783 """ 784 test_load_csv_lat_long_bad_blocking(self): 709 '''test_load_csv_lat_long_bad_blocking(self): 785 710 Different zones in Geo references 786 """ 787 fileName = tempfile.mktemp(".csv") 788 file = open(fileName,"w") 789 file.write("Lati,LONG,z \n\ 711 ''' 712 713 fileName = tempfile.mktemp('.csv') 714 file = open(fileName, 'w') 715 file.write('Lati,LONG,z \n\ 790 716 -25.0,180.0,452.688000\n\ 791 -34,150.0,459.126000\n ")717 -34,150.0,459.126000\n') 792 718 file.close() 793 719 794 720 results = Geospatial_data(fileName, max_read_lines=1, 795 721 load_file_now=False) 796 797 #for i in results: 798 # pass 722 799 723 try: 800 724 for i in results: … … 804 728 else: 805 729 msg = 'Different zones in Geo references not caught.' 806 raise msg807 808 os.remove(fileName) 809 730 raise Exception, msg 731 732 os.remove(fileName) 733 810 734 def test_load_csv(self): 811 812 fileName = tempfile.mktemp(".txt") 813 file = open(fileName,"w") 814 file.write(" x,y, elevation , speed \n\ 815 1.0, 0.0, 10.0, 0.0\n\ 816 0.0, 1.0, 0.0, 10.0\n\ 817 1.0, 0.0 ,10.4, 40.0\n") 735 fileName = tempfile.mktemp('.txt') 736 file = open(fileName, 'w') 737 file.write(' x,y, elevation , speed \n\ 738 1.0, 0.0, 10.0, 0.0\n\ 739 0.0, 1.0, 0.0, 10.0\n\ 740 1.0, 0.0 ,10.4, 40.0\n') 818 741 file.close() 819 742 820 743 results = Geospatial_data(fileName, max_read_lines=2) 821 744 822 823 assert num.allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]) 824 assert num.allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4]) 825 assert num.allclose(results.get_attributes(attribute_name='speed'), [0.0, 10.0, 40.0]) 745 assert num.allclose(results.get_data_points(), 746 [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]) 747 assert num.allclose(results.get_attributes(attribute_name='elevation'), 748 [10.0, 0.0, 10.4]) 749 assert num.allclose(results.get_attributes(attribute_name='speed'), 750 [0.0, 10.0, 40.0]) 826 751 827 752 # Blocking … … 829 754 for i in results: 830 755 geo_list.append(i) 831 756 832 757 assert num.allclose(geo_list[0].get_data_points(), 833 758 [[1.0, 0.0],[0.0, 1.0]]) 834 835 759 assert num.allclose(geo_list[0].get_attributes(attribute_name='elevation'), 836 760 [10.0, 0.0]) 837 761 assert num.allclose(geo_list[1].get_data_points(), 838 [[1.0, 0.0]]) 762 [[1.0, 0.0]]) 839 763 assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'), 840 764 [10.4]) 841 842 os.remove(fileName) 843 765 766 os.remove(fileName) 767 844 768 def test_load_csv_bad(self): 845 """test_load_csv_bad(self):769 '''test_load_csv_bad(self): 846 770 header column, body column missmatch 847 771 (Format error) 848 """ 849 import os 850 851 fileName = tempfile.mktemp(".txt") 852 file = open(fileName,"w") 853 file.write(" elevation , speed \n\ 854 1.0, 0.0, 10.0, 0.0\n\ 855 0.0, 1.0, 0.0, 10.0\n\ 856 1.0, 0.0 ,10.4, 40.0\n") 772 ''' 773 774 fileName = tempfile.mktemp('.txt') 775 file = open(fileName, 'w') 776 file.write(' elevation , speed \n\ 777 1.0, 0.0, 10.0, 0.0\n\ 778 0.0, 1.0, 0.0, 10.0\n\ 779 1.0, 0.0 ,10.4, 40.0\n') 857 780 file.close() 858 781 … … 862 785 # Blocking 863 786 geo_list = [] 864 #for i in results:865 # geo_list.append(i)866 787 try: 867 788 for i in results: … … 870 791 pass 871 792 else: 872 msg = ' bad file did not raise error!'873 raise msg793 msg = 'Bad file did not raise error!' 794 raise Exception, msg 874 795 os.remove(fileName) 875 796 876 797 def test_load_csv_badII(self): 877 """test_load_csv_bad(self):798 '''test_load_csv_bad(self): 878 799 header column, body column missmatch 879 800 (Format error) 880 """ 881 import os 882 883 fileName = tempfile.mktemp(".txt") 884 file = open(fileName,"w") 885 file.write(" x,y,elevation , speed \n\ 801 ''' 802 803 fileName = tempfile.mktemp('.txt') 804 file = open(fileName, 'w') 805 file.write(' x,y,elevation , speed \n\ 886 806 1.0, 0.0, 10.0, 0.0\n\ 887 807 0.0, 1.0, 0.0, 10\n\ 888 1.0, 0.0 ,10.4 yeah\n ")808 1.0, 0.0 ,10.4 yeah\n') 889 809 file.close() 890 810 … … 894 814 # Blocking 895 815 geo_list = [] 896 #for i in results:897 # geo_list.append(i)898 816 try: 899 817 for i in results: … … 902 820 pass 903 821 else: 904 msg = ' bad file did not raise error!'905 raise msg822 msg = 'Bad file did not raise error!' 823 raise Exception, msg 906 824 os.remove(fileName) 907 825 908 826 def test_load_csv_badIII(self): 909 """test_load_csv_bad(self):827 '''test_load_csv_bad(self): 910 828 space delimited 911 """ 912 import os 913 914 fileName = tempfile.mktemp(".txt") 915 file = open(fileName,"w") 916 file.write(" x y elevation speed \n\ 829 ''' 830 831 fileName = tempfile.mktemp('.txt') 832 file = open(fileName, 'w') 833 file.write(' x y elevation speed \n\ 917 834 1. 0.0 10.0 0.0\n\ 918 835 0.0 1.0 0.0 10.0\n\ 919 1.0 0.0 10.4 40.0\n ")836 1.0 0.0 10.4 40.0\n') 920 837 file.close() 921 838 … … 926 843 pass 927 844 else: 928 msg = ' bad file did not raise error!'929 raise msg930 os.remove(fileName) 931 845 msg = 'Bad file did not raise error!' 846 raise Exception, msg 847 os.remove(fileName) 848 932 849 def test_load_csv_badIV(self): 933 """test_load_csv_bad(self):850 ''' test_load_csv_bad(self): 934 851 header column, body column missmatch 935 852 (Format error) 936 """ 937 import os 938 939 fileName = tempfile.mktemp(".txt") 940 file = open(fileName,"w") 941 file.write(" x,y,elevation , speed \n\ 853 ''' 854 855 fileName = tempfile.mktemp('.txt') 856 file = open(fileName, 'w') 857 file.write(' x,y,elevation , speed \n\ 942 858 1.0, 0.0, 10.0, wow\n\ 943 859 0.0, 1.0, 0.0, ha\n\ 944 1.0, 0.0 ,10.4, yeah\n ")860 1.0, 0.0 ,10.4, yeah\n') 945 861 file.close() 946 862 … … 950 866 # Blocking 951 867 geo_list = [] 952 #for i in results:953 # geo_list.append(i)954 868 try: 955 869 for i in results: … … 958 872 pass 959 873 else: 960 msg = ' bad file did not raise error!'961 raise msg874 msg = 'Bad file did not raise error!' 875 raise Exception, msg 962 876 os.remove(fileName) 963 877 964 878 def test_load_pts_blocking(self): 965 879 #This is pts! 966 967 import os 968 969 fileName = tempfile.mktemp(".txt") 970 file = open(fileName,"w") 971 file.write(" x,y, elevation , speed \n\ 972 1.0, 0.0, 10.0, 0.0\n\ 973 0.0, 1.0, 0.0, 10.0\n\ 974 1.0, 0.0 ,10.4, 40.0\n") 880 fileName = tempfile.mktemp('.txt') 881 file = open(fileName, 'w') 882 file.write(' x,y, elevation , speed \n\ 883 1.0, 0.0, 10.0, 0.0\n\ 884 0.0, 1.0, 0.0, 10.0\n\ 885 1.0, 0.0 ,10.4, 40.0\n') 975 886 file.close() 976 887 977 pts_file = tempfile.mktemp( ".pts")978 888 pts_file = tempfile.mktemp('.pts') 889 979 890 convert = Geospatial_data(fileName) 980 891 convert.export_points_file(pts_file) 981 892 results = Geospatial_data(pts_file, max_read_lines=2) 982 893 983 assert num.allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],984 894 assert num.allclose(results.get_data_points(), 895 [[1.0, 0.0],[0.0, 1.0], [1.0, 0.0]]) 985 896 assert num.allclose(results.get_attributes(attribute_name='elevation'), 986 897 [10.0, 0.0, 10.4]) … … 991 902 geo_list = [] 992 903 for i in results: 993 geo_list.append(i) 904 geo_list.append(i) 994 905 assert num.allclose(geo_list[0].get_data_points(), 995 906 [[1.0, 0.0],[0.0, 1.0]]) … … 997 908 [10.0, 0.0]) 998 909 assert num.allclose(geo_list[1].get_data_points(), 999 [[1.0, 0.0]]) 910 [[1.0, 0.0]]) 1000 911 assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'), 1001 912 [10.4]) 1002 1003 os.remove(fileName) 1004 os.remove(pts_file) 913 914 os.remove(fileName) 915 os.remove(pts_file) 1005 916 1006 917 def verbose_test_load_pts_blocking(self): 1007 1008 import os 1009 1010 fileName = tempfile.mktemp(".txt") 1011 file = open(fileName,"w") 1012 file.write(" x,y, elevation , speed \n\ 1013 1.0, 0.0, 10.0, 0.0\n\ 1014 0.0, 1.0, 0.0, 10.0\n\ 1015 1.0, 0.0, 10.0, 0.0\n\ 1016 0.0, 1.0, 0.0, 10.0\n\ 1017 1.0, 0.0, 10.0, 0.0\n\ 1018 0.0, 1.0, 0.0, 10.0\n\ 1019 1.0, 0.0, 10.0, 0.0\n\ 1020 0.0, 1.0, 0.0, 10.0\n\ 1021 1.0, 0.0, 10.0, 0.0\n\ 1022 0.0, 1.0, 0.0, 10.0\n\ 1023 1.0, 0.0, 10.0, 0.0\n\ 1024 0.0, 1.0, 0.0, 10.0\n\ 1025 1.0, 0.0, 10.0, 0.0\n\ 1026 0.0, 1.0, 0.0, 10.0\n\ 1027 1.0, 0.0, 10.0, 0.0\n\ 1028 0.0, 1.0, 0.0, 10.0\n\ 1029 1.0, 0.0, 10.0, 0.0\n\ 1030 0.0, 1.0, 0.0, 10.0\n\ 1031 1.0, 0.0, 10.0, 0.0\n\ 1032 0.0, 1.0, 0.0, 10.0\n\ 1033 1.0, 0.0, 10.0, 0.0\n\ 1034 0.0, 1.0, 0.0, 10.0\n\ 1035 1.0, 0.0, 10.0, 0.0\n\ 1036 0.0, 1.0, 0.0, 10.0\n\ 1037 1.0, 0.0, 10.0, 0.0\n\ 1038 0.0, 1.0, 0.0, 10.0\n\ 1039 1.0, 0.0, 10.0, 0.0\n\ 1040 0.0, 1.0, 0.0, 10.0\n\ 1041 1.0, 0.0, 10.0, 0.0\n\ 1042 0.0, 1.0, 0.0, 10.0\n\ 1043 1.0, 0.0, 10.0, 0.0\n\ 1044 0.0, 1.0, 0.0, 10.0\n\ 1045 1.0, 0.0, 10.0, 0.0\n\ 1046 0.0, 1.0, 0.0, 10.0\n\ 1047 1.0, 0.0, 10.0, 0.0\n\ 1048 0.0, 1.0, 0.0, 10.0\n\ 1049 1.0, 0.0, 10.0, 0.0\n\ 1050 0.0, 1.0, 0.0, 10.0\n\ 1051 1.0, 0.0, 10.0, 0.0\n\ 1052 0.0, 1.0, 0.0, 10.0\n\ 1053 1.0, 0.0, 10.0, 0.0\n\ 1054 0.0, 1.0, 0.0, 10.0\n\ 1055 1.0, 0.0, 10.0, 0.0\n\ 1056 0.0, 1.0, 0.0, 10.0\n\ 1057 1.0, 0.0, 10.0, 0.0\n\ 1058 0.0, 1.0, 0.0, 10.0\n\ 1059 1.0, 0.0 ,10.4, 40.0\n") 918 fileName = tempfile.mktemp('.txt') 919 file = open(fileName, 'w') 920 file.write(' x,y, elevation , speed \n\ 921 1.0, 0.0, 10.0, 0.0\n\ 922 0.0, 1.0, 0.0, 10.0\n\ 923 1.0, 0.0, 10.0, 0.0\n\ 924 0.0, 1.0, 0.0, 10.0\n\ 925 1.0, 0.0, 10.0, 0.0\n\ 926 0.0, 1.0, 0.0, 10.0\n\ 927 1.0, 0.0, 10.0, 0.0\n\ 928 0.0, 1.0, 0.0, 10.0\n\ 929 1.0, 0.0, 10.0, 0.0\n\ 930 0.0, 1.0, 0.0, 10.0\n\ 931 1.0, 0.0, 10.0, 0.0\n\ 932 0.0, 1.0, 0.0, 10.0\n\ 933 1.0, 0.0, 10.0, 0.0\n\ 934 0.0, 1.0, 0.0, 10.0\n\ 935 1.0, 0.0, 10.0, 0.0\n\ 936 0.0, 1.0, 0.0, 10.0\n\ 937 1.0, 0.0, 10.0, 0.0\n\ 938 0.0, 1.0, 0.0, 10.0\n\ 939 1.0, 0.0, 10.0, 0.0\n\ 940 0.0, 1.0, 0.0, 10.0\n\ 941 1.0, 0.0, 10.0, 0.0\n\ 942 0.0, 1.0, 0.0, 10.0\n\ 943 1.0, 0.0, 10.0, 0.0\n\ 944 0.0, 1.0, 0.0, 10.0\n\ 945 1.0, 0.0, 10.0, 0.0\n\ 946 0.0, 1.0, 0.0, 10.0\n\ 947 1.0, 0.0, 10.0, 0.0\n\ 948 0.0, 1.0, 0.0, 10.0\n\ 949 1.0, 0.0, 10.0, 0.0\n\ 950 0.0, 1.0, 0.0, 10.0\n\ 951 1.0, 0.0, 10.0, 0.0\n\ 952 0.0, 1.0, 0.0, 10.0\n\ 953 1.0, 0.0, 10.0, 0.0\n\ 954 0.0, 1.0, 0.0, 10.0\n\ 955 1.0, 0.0, 10.0, 0.0\n\ 956 0.0, 1.0, 0.0, 10.0\n\ 957 1.0, 0.0, 10.0, 0.0\n\ 958 0.0, 1.0, 0.0, 10.0\n\ 959 1.0, 0.0, 10.0, 0.0\n\ 960 0.0, 1.0, 0.0, 10.0\n\ 961 1.0, 0.0, 10.0, 0.0\n\ 962 0.0, 1.0, 0.0, 10.0\n\ 963 1.0, 0.0, 10.0, 0.0\n\ 964 0.0, 1.0, 0.0, 10.0\n\ 965 1.0, 0.0, 10.0, 0.0\n\ 966 0.0, 1.0, 0.0, 10.0\n\ 967 1.0, 0.0 ,10.4, 40.0\n') 1060 968 file.close() 1061 969 1062 pts_file = tempfile.mktemp(".pts") 1063 970 pts_file = tempfile.mktemp('.pts') 1064 971 convert = Geospatial_data(fileName) 1065 972 convert.export_points_file(pts_file) … … 1069 976 geo_list = [] 1070 977 for i in results: 1071 geo_list.append(i) 978 geo_list.append(i) 1072 979 assert num.allclose(geo_list[0].get_data_points(), 1073 [[1.0, 0.0], [0.0, 1.0]])980 [[1.0, 0.0], [0.0, 1.0]]) 1074 981 assert num.allclose(geo_list[0].get_attributes(attribute_name='elevation'), 1075 982 [10.0, 0.0]) 1076 983 assert num.allclose(geo_list[1].get_data_points(), 1077 [[1.0, 0.0],[0.0, 1.0] ]) 984 [[1.0, 0.0],[0.0, 1.0] ]) 1078 985 assert num.allclose(geo_list[1].get_attributes(attribute_name='elevation'), 1079 986 [10.0, 0.0]) 1080 1081 os.remove(fileName) 987 988 os.remove(fileName) 1082 989 os.remove(pts_file) 1083 1084 1085 990 1086 991 def test_new_export_pts_file(self): … … 1089 994 att_dict['elevation'] = num.array([10.1, 0.0, 10.4]) 1090 995 att_dict['brightness'] = num.array([10.0, 1.0, 10.4]) 1091 1092 fileName = tempfile.mktemp(".pts") 1093 996 997 fileName = tempfile.mktemp('.pts') 1094 998 G = Geospatial_data(pointlist, att_dict) 1095 1096 999 G.export_points_file(fileName) 1097 1098 1000 results = Geospatial_data(file_name = fileName) 1099 1100 os.remove(fileName) 1101 1102 assert num.allclose(results.get_data_points(),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]) 1103 assert num.allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4]) 1001 os.remove(fileName) 1002 1003 assert num.allclose(results.get_data_points(), 1004 [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1005 assert num.allclose(results.get_attributes(attribute_name='elevation'), 1006 [10.1, 0.0, 10.4]) 1104 1007 answer = [10.0, 1.0, 10.4] 1105 assert num.allclose(results.get_attributes(attribute_name='brightness'), answer) 1008 assert num.allclose(results.get_attributes(attribute_name='brightness'), 1009 answer) 1106 1010 1107 1011 def test_new_export_absolute_pts_file(self): 1108 1012 att_dict = {} 1109 pointlist = num.array([[1.0, 0.0], [0.0, 1.0],[1.0, 0.0]])1013 pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1110 1014 att_dict['elevation'] = num.array([10.1, 0.0, 10.4]) 1111 1015 att_dict['brightness'] = num.array([10.0, 1.0, 10.4]) 1112 1016 geo_ref = Geo_reference(50, 25, 55) 1113 1114 fileName = tempfile.mktemp(".pts") 1115 1017 1018 fileName = tempfile.mktemp('.pts') 1116 1019 G = Geospatial_data(pointlist, att_dict, geo_ref) 1117 1118 1020 G.export_points_file(fileName, absolute=True) 1119 1120 1021 results = Geospatial_data(file_name = fileName) 1121 1122 os.remove(fileName) 1123 1124 assert num.allclose(results.get_data_points(), G.get_data_points(True)) 1125 assert num.allclose(results.get_attributes(attribute_name='elevation'), [10.1, 0.0, 10.4]) 1022 os.remove(fileName) 1023 1024 msg = ('results.get_data_points()=%s, but\nG.get_data_points(True)=%s' 1025 % (str(results.get_data_points()), str(G.get_data_points(True)))) 1026 assert num.allclose(results.get_data_points(), 1027 G.get_data_points(True)), msg 1028 msg = ("results.get_attributes(attribute_name='elevation')=%s" 1029 % str(results.get_attributes(attribute_name='elevation'))) 1030 assert num.allclose(results.get_attributes(attribute_name='elevation'), 1031 [10.1, 0.0, 10.4]), msg 1126 1032 answer = [10.0, 1.0, 10.4] 1127 assert num.allclose(results.get_attributes(attribute_name='brightness'), answer) 1033 msg = ("results.get_attributes(attribute_name='brightness')=%s, " 1034 'answer=%s' % 1035 (str(results.get_attributes(attribute_name='brightness')), 1036 str(answer))) 1037 assert num.allclose(results.get_attributes(attribute_name='brightness'), 1038 answer), msg 1128 1039 1129 1040 def test_loadpts(self): 1130 1131 1041 from Scientific.IO.NetCDF import NetCDFFile 1132 1042 1133 fileName = tempfile.mktemp( ".pts")1043 fileName = tempfile.mktemp('.pts') 1134 1044 # NetCDF file definition 1135 1045 outfile = NetCDFFile(fileName, netcdf_mode_w) 1136 1046 1137 1047 # dimension definitions 1138 outfile.createDimension('number_of_points', 3) 1139 outfile.createDimension('number_of_dimensions', 2) # This is 2d data1140 1048 outfile.createDimension('number_of_points', 3) 1049 outfile.createDimension('number_of_dimensions', 2) # This is 2d data 1050 1141 1051 # variable definitions 1142 outfile.createVariable('points', n um.Float, ('number_of_points',1143 'number_of_dimensions'))1144 outfile.createVariable('elevation', n um.Float, ('number_of_points',))1145 1052 outfile.createVariable('points', netcdf_float, ('number_of_points', 1053 'number_of_dimensions')) 1054 outfile.createVariable('elevation', netcdf_float, ('number_of_points',)) 1055 1146 1056 # Get handles to the variables 1147 1057 points = outfile.variables['points'] 1148 1058 elevation = outfile.variables['elevation'] 1149 1059 1150 1060 points[0, :] = [1.0,0.0] 1151 elevation[0] = 10.0 1061 elevation[0] = 10.0 1152 1062 points[1, :] = [0.0,1.0] 1153 elevation[1] = 0.0 1063 elevation[1] = 0.0 1154 1064 points[2, :] = [1.0,0.0] 1155 elevation[2] = 10.4 1065 elevation[2] = 10.4 1156 1066 1157 1067 outfile.close() 1158 1068 1159 1069 results = Geospatial_data(file_name = fileName) 1160 1070 os.remove(fileName) 1161 answer = [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]] 1162 assert num.allclose(results.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]) 1163 assert num.allclose(results.get_attributes(attribute_name='elevation'), [10.0, 0.0, 10.4]) 1164 1071 1072 answer = [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]] 1073 assert num.allclose(results.get_data_points(), 1074 [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1075 assert num.allclose(results.get_attributes(attribute_name='elevation'), 1076 [10.0, 0.0, 10.4]) 1077 1165 1078 def test_writepts(self): 1166 #test_writepts: Test that storage of x,y,attributes works1167 1079 '''Test that storage of x,y,attributes works''' 1080 1168 1081 att_dict = {} 1169 pointlist = num.array([[1.0, 0.0], [0.0, 1.0],[1.0, 0.0]])1082 pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1170 1083 att_dict['elevation'] = num.array([10.0, 0.0, 10.4]) 1171 1084 att_dict['brightness'] = num.array([10.0, 0.0, 10.4]) 1172 geo_reference=Geo_reference(56, 1.9,1.9)1085 geo_reference=Geo_reference(56, 1.9, 1.9) 1173 1086 1174 1087 # Test pts format 1175 fileName = tempfile.mktemp( ".pts")1088 fileName = tempfile.mktemp('.pts') 1176 1089 G = Geospatial_data(pointlist, att_dict, geo_reference) 1177 1090 G.export_points_file(fileName, False) … … 1179 1092 os.remove(fileName) 1180 1093 1181 assert num.allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]) 1182 assert num.allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4]) 1094 assert num.allclose(results.get_data_points(False), 1095 [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1096 assert num.allclose(results.get_attributes('elevation'), 1097 [10.0, 0.0, 10.4]) 1183 1098 answer = [10.0, 0.0, 10.4] 1184 1099 assert num.allclose(results.get_attributes('brightness'), answer) 1185 1100 self.failUnless(geo_reference == geo_reference, 1186 1101 'test_writepts failed. Test geo_reference') 1187 1102 1188 1103 def test_write_csv_attributes(self): 1189 #test_write : Test that storage of x,y,attributes works1190 1104 '''Test that storage of x,y,attributes works''' 1105 1191 1106 att_dict = {} 1192 pointlist = num.array([[1.0, 0.0], [0.0, 1.0],[1.0, 0.0]])1107 pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1193 1108 att_dict['elevation'] = num.array([10.0, 0.0, 10.4]) 1194 1109 att_dict['brightness'] = num.array([10.0, 0.0, 10.4]) 1195 geo_reference=Geo_reference(56,0,0) 1110 geo_reference=Geo_reference(56, 0, 0) 1111 1196 1112 # Test txt format 1197 fileName = tempfile.mktemp( ".txt")1113 fileName = tempfile.mktemp('.txt') 1198 1114 G = Geospatial_data(pointlist, att_dict, geo_reference) 1199 1115 G.export_points_file(fileName) 1200 #print "fileName", fileName1201 1116 results = Geospatial_data(file_name=fileName) 1202 1117 os.remove(fileName) 1203 assert num.allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]) 1204 assert num.allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4]) 1118 1119 assert num.allclose(results.get_data_points(False), 1120 [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1121 assert num.allclose(results.get_attributes('elevation'), 1122 [10.0, 0.0, 10.4]) 1205 1123 answer = [10.0, 0.0, 10.4] 1206 1124 assert num.allclose(results.get_attributes('brightness'), answer) 1207 1208 1125 1209 1126 def test_write_csv_attributes_lat_long(self): 1210 #test_write : Test that storage of x,y,attributes works1211 1127 '''Test that storage of x,y,attributes works''' 1128 1212 1129 att_dict = {} 1213 pointlist = num.array([[-21.5, 114.5],[-21.6,114.5],[-21.7,114.5]])1130 pointlist = num.array([[-21.5, 114.5], [-21.6, 114.5], [-21.7, 114.5]]) 1214 1131 att_dict['elevation'] = num.array([10.0, 0.0, 10.4]) 1215 1132 att_dict['brightness'] = num.array([10.0, 0.0, 10.4]) 1133 1216 1134 # Test txt format 1217 fileName = tempfile.mktemp( ".txt")1135 fileName = tempfile.mktemp('.txt') 1218 1136 G = Geospatial_data(pointlist, att_dict, points_are_lats_longs=True) 1219 1137 G.export_points_file(fileName, as_lat_long=True) 1220 #print "fileName", fileName1221 1138 results = Geospatial_data(file_name=fileName) 1222 1139 os.remove(fileName) 1140 1223 1141 assert num.allclose(results.get_data_points(False, as_lat_long=True), 1224 pointlist) 1225 assert num.allclose(results.get_attributes('elevation'), [10.0, 0.0, 10.4]) 1142 pointlist) 1143 assert num.allclose(results.get_attributes('elevation'), 1144 [10.0, 0.0, 10.4]) 1226 1145 answer = [10.0, 0.0, 10.4] 1227 1146 assert num.allclose(results.get_attributes('brightness'), answer) 1228 1147 1229 1148 def test_writepts_no_attributes(self): 1230 1231 #test_writepts_no_attributes: Test that storage of x,y alone works 1232 1149 '''Test that storage of x,y alone works''' 1150 1233 1151 att_dict = {} 1234 pointlist = num.array([[1.0, 0.0], [0.0, 1.0],[1.0, 0.0]])1235 geo_reference=Geo_reference(56, 1.9,1.9)1152 pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1153 geo_reference=Geo_reference(56, 1.9, 1.9) 1236 1154 1237 1155 # Test pts format 1238 fileName = tempfile.mktemp( ".pts")1156 fileName = tempfile.mktemp('.pts') 1239 1157 G = Geospatial_data(pointlist, None, geo_reference) 1240 1158 G.export_points_file(fileName, False) … … 1242 1160 os.remove(fileName) 1243 1161 1244 assert num.allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]) 1162 assert num.allclose(results.get_data_points(False), 1163 [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1245 1164 self.failUnless(geo_reference == geo_reference, 1246 'test_writepts failed. Test geo_reference') 1247 1248 1165 'test_writepts failed. Test geo_reference') 1166 1249 1167 def test_write_csv_no_attributes(self): 1250 #test_write txt _no_attributes: Test that storage of x,y alone works1251 1168 '''Test that storage of x,y alone works''' 1169 1252 1170 att_dict = {} 1253 pointlist = num.array([[1.0, 0.0], [0.0, 1.0],[1.0, 0.0]])1171 pointlist = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1254 1172 geo_reference=Geo_reference(56,0,0) 1173 1255 1174 # Test format 1256 fileName = tempfile.mktemp( ".txt")1175 fileName = tempfile.mktemp('.txt') 1257 1176 G = Geospatial_data(pointlist, None, geo_reference) 1258 1177 G.export_points_file(fileName) 1259 1178 results = Geospatial_data(file_name=fileName) 1260 1179 os.remove(fileName) 1261 assert num.allclose(results.get_data_points(False),[[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]]) 1262 1263 1264 1265 ########################## BAD .PTS ########################## 1180 1181 assert num.allclose(results.get_data_points(False), 1182 [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1183 1184 ################################################################################ 1185 # Check bad PTS files. 1186 ################################################################################ 1266 1187 1267 1188 def test_load_bad_no_file_pts(self): 1268 import os 1269 import tempfile 1270 1271 fileName = tempfile.mktemp(".pts") 1272 #print fileName 1189 fileName = tempfile.mktemp('.pts') 1273 1190 try: 1274 results = Geospatial_data(file_name = fileName) 1275 # dict = import_points_file(fileName) 1191 results = Geospatial_data(file_name=fileName) 1276 1192 except IOError: 1277 1193 pass 1278 1194 else: 1279 1195 msg = 'imaginary file did not raise error!' 1280 raise msg 1281 # self.failUnless(0 == 1, 1282 # 'imaginary file did not raise error!') 1283 1196 raise Exception, msg 1284 1197 1285 1198 def test_create_from_pts_file(self): 1286 1287 1199 from Scientific.IO.NetCDF import NetCDFFile 1288 1200 1289 # fileName = tempfile.mktemp(".pts") 1201 # NetCDF file definition 1290 1202 FN = 'test_points.pts' 1291 # NetCDF file definition1292 1203 outfile = NetCDFFile(FN, netcdf_mode_w) 1293 1204 1294 1205 # dimension definitions 1295 outfile.createDimension('number_of_points', 3) 1296 outfile.createDimension('number_of_dimensions', 2) # This is 2d data1297 1206 outfile.createDimension('number_of_points', 3) 1207 outfile.createDimension('number_of_dimensions', 2) # This is 2d data 1208 1298 1209 # variable definitions 1299 outfile.createVariable('points', n um.Float, ('number_of_points',1300 'number_of_dimensions'))1301 outfile.createVariable('elevation', n um.Float, ('number_of_points',))1302 1210 outfile.createVariable('points', netcdf_float, ('number_of_points', 1211 'number_of_dimensions')) 1212 outfile.createVariable('elevation', netcdf_float, ('number_of_points',)) 1213 1303 1214 # Get handles to the variables 1304 1215 points = outfile.variables['points'] 1305 1216 elevation = outfile.variables['elevation'] 1306 1217 1307 1218 points[0, :] = [1.0,0.0] 1308 elevation[0] = 10.0 1219 elevation[0] = 10.0 1309 1220 points[1, :] = [0.0,1.0] 1310 elevation[1] = 0.0 1221 elevation[1] = 0.0 1311 1222 points[2, :] = [1.0,0.0] 1312 elevation[2] = 10.4 1223 elevation[2] = 10.4 1313 1224 1314 1225 outfile.close() … … 1318 1229 assert num.allclose(G.get_geo_reference().get_xllcorner(), 0.0) 1319 1230 assert num.allclose(G.get_geo_reference().get_yllcorner(), 0.0) 1320 1321 assert num.allclose(G.get_data_points(), [[1.0, 0.0],[0.0, 1.0],[1.0, 0.0]])1231 assert num.allclose(G.get_data_points(), 1232 [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1322 1233 assert num.allclose(G.get_attributes(), [10.0, 0.0, 10.4]) 1323 1234 os.remove(FN) 1324 1235 1325 1236 def test_create_from_pts_file_with_geo(self): 1326 """This test reveals if Geospatial data is correctly instantiated from a pts file. 1327 """ 1328 1237 '''Test if Geospatial data is correctly instantiated from a pts file.''' 1238 1329 1239 from Scientific.IO.NetCDF import NetCDFFile 1330 1240 1241 # NetCDF file definition 1331 1242 FN = 'test_points.pts' 1332 # NetCDF file definition1333 1243 outfile = NetCDFFile(FN, netcdf_mode_w) 1334 1244 … … 1340 1250 1341 1251 # dimension definitions 1342 outfile.createDimension('number_of_points', 3) 1343 outfile.createDimension('number_of_dimensions', 2) # This is 2d data1344 1252 outfile.createDimension('number_of_points', 3) 1253 outfile.createDimension('number_of_dimensions', 2) # This is 2d data 1254 1345 1255 # variable definitions 1346 outfile.createVariable('points', n um.Float, ('number_of_points',1347 'number_of_dimensions'))1348 outfile.createVariable('elevation', n um.Float, ('number_of_points',))1349 1256 outfile.createVariable('points', netcdf_float, ('number_of_points', 1257 'number_of_dimensions')) 1258 outfile.createVariable('elevation', netcdf_float, ('number_of_points',)) 1259 1350 1260 # Get handles to the variables 1351 1261 points = outfile.variables['points'] … … 1353 1263 1354 1264 points[0, :] = [1.0,0.0] 1355 elevation[0] = 10.0 1265 elevation[0] = 10.0 1356 1266 points[1, :] = [0.0,1.0] 1357 elevation[1] = 0.0 1267 elevation[1] = 0.0 1358 1268 points[2, :] = [1.0,0.0] 1359 elevation[2] = 10.4 1269 elevation[2] = 10.4 1360 1270 1361 1271 outfile.close() … … 1365 1275 assert num.allclose(G.get_geo_reference().get_xllcorner(), xll) 1366 1276 assert num.allclose(G.get_geo_reference().get_yllcorner(), yll) 1367 1368 1277 assert num.allclose(G.get_data_points(), [[1.0+xll, 0.0+yll], 1369 1278 [0.0+xll, 1.0+yll], 1370 1279 [1.0+xll, 0.0+yll]]) 1371 1372 1280 assert num.allclose(G.get_attributes(), [10.0, 0.0, 10.4]) 1281 1373 1282 os.remove(FN) 1374 1283 1375 1376 1284 def test_add_(self): 1377 1285 '''test_add_(self): 1378 1286 adds an txt and pts files, reads the files and adds them 1379 1287 checking results are correct 1380 1288 ''' 1289 1381 1290 # create files 1382 1291 att_dict1 = {} 1383 pointlist1 = num.array([[1.0, 0.0], [0.0, 1.0],[1.0, 0.0]])1292 pointlist1 = num.array([[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1384 1293 att_dict1['elevation'] = num.array([-10.0, 0.0, 10.4]) 1385 1294 att_dict1['brightness'] = num.array([10.0, 0.0, 10.4]) 1386 1295 geo_reference1 = Geo_reference(56, 2.0, 1.0) 1387 1296 1388 1297 att_dict2 = {} 1389 pointlist2 = num.array([[2.0, 1.0], [1.0, 2.0],[2.0, 1.0]])1298 pointlist2 = num.array([[2.0, 1.0], [1.0, 2.0], [2.0, 1.0]]) 1390 1299 att_dict2['elevation'] = num.array([1.0, 15.0, 1.4]) 1391 1300 att_dict2['brightness'] = num.array([14.0, 1.0, -12.4]) 1392 geo_reference2 = Geo_reference(56, 1.0, 2.0) 1301 geo_reference2 = Geo_reference(56, 1.0, 2.0) 1393 1302 1394 1303 G1 = Geospatial_data(pointlist1, att_dict1, geo_reference1) 1395 1304 G2 = Geospatial_data(pointlist2, att_dict2, geo_reference2) 1396 1397 fileName1 = tempfile.mktemp( ".txt")1398 fileName2 = tempfile.mktemp( ".pts")1399 1400 # makes files1305 1306 fileName1 = tempfile.mktemp('.txt') 1307 fileName2 = tempfile.mktemp('.pts') 1308 1309 # makes files 1401 1310 G1.export_points_file(fileName1) 1402 1311 G2.export_points_file(fileName2) 1403 1312 1404 1313 # add files 1405 1406 G3 = Geospatial_data(file_name = fileName1) 1407 G4 = Geospatial_data(file_name = fileName2) 1408 1314 G3 = Geospatial_data(file_name=fileName1) 1315 G4 = Geospatial_data(file_name=fileName2) 1409 1316 G = G3 + G4 1410 1317 1411 1412 1318 #read results 1413 # print'res', G.get_data_points()1414 # print'res1', G.get_data_points(False)1415 1319 assert num.allclose(G.get_data_points(), 1416 1320 [[ 3.0, 1.0], [ 2.0, 2.0], 1417 1321 [ 3.0, 1.0], [ 3.0, 3.0], 1418 1322 [ 2.0, 4.0], [ 3.0, 3.0]]) 1419 1420 1323 assert num.allclose(G.get_attributes(attribute_name='elevation'), 1421 1324 [-10.0, 0.0, 10.4, 1.0, 15.0, 1.4]) 1422 1423 1325 answer = [10.0, 0.0, 10.4, 14.0, 1.0, -12.4] 1424 assert num.allclose(G.get_attributes(attribute_name='brightness'), answer)1425 1326 assert num.allclose(G.get_attributes(attribute_name='brightness'), 1327 answer) 1426 1328 self.failUnless(G.get_geo_reference() == geo_reference1, 1427 1428 1329 'test_writepts failed. Test geo_reference') 1330 1429 1331 os.remove(fileName1) 1430 1332 os.remove(fileName2) 1431 1333 1432 1334 def test_ensure_absolute(self): 1433 points = [[2.0, 0.0], [1.0, 1.0],1434 [2.0, 0.0],[2.0, 2.0],1435 [1.0, 3.0],[2.0, 2.0]]1335 points = [[2.0, 0.0], [1.0, 1.0], 1336 [2.0, 0.0], [2.0, 2.0], 1337 [1.0, 3.0], [2.0, 2.0]] 1436 1338 new_points = ensure_absolute(points) 1437 1339 1438 1340 assert num.allclose(new_points, points) 1439 1341 1440 points = num.array([[2.0, 0.0], [1.0, 1.0],1441 [2.0, 0.0], [2.0, 2.0],1442 [1.0, 3.0], [2.0, 2.0]])1342 points = num.array([[2.0, 0.0], [1.0, 1.0], 1343 [2.0, 0.0], [2.0, 2.0], 1344 [1.0, 3.0], [2.0, 2.0]]) 1443 1345 new_points = ensure_absolute(points) 1444 1346 1445 1347 assert num.allclose(new_points, points) 1446 1447 ab_points = num.array([[2.0, 0.0], [1.0, 1.0],1448 [2.0, 0.0], [2.0, 2.0],1449 [1.0, 3.0], [2.0, 2.0]])1450 1451 mesh_origin = (56, 290000, 618000) # zone, easting, northing1452 1453 data_points = num.zeros((ab_points.shape), num.Float) 1348 1349 ab_points = num.array([[2.0, 0.0], [1.0, 1.0], 1350 [2.0, 0.0], [2.0, 2.0], 1351 [1.0, 3.0], [2.0, 2.0]]) 1352 1353 mesh_origin = (56, 290000, 618000) # zone, easting, northing 1354 data_points = num.zeros((ab_points.shape), num.float) 1355 1454 1356 #Shift datapoints according to new origins 1455 1357 for k in range(len(ab_points)): 1456 1358 data_points[k][0] = ab_points[k][0] - mesh_origin[1] 1457 1359 data_points[k][1] = ab_points[k][1] - mesh_origin[2] 1458 #print "data_points",data_points 1459 new_points = ensure_absolute(data_points, 1460 geo_reference=mesh_origin) 1461 #print "new_points",new_points 1462 #print "ab_points",ab_points 1463 1360 new_points = ensure_absolute(data_points, geo_reference=mesh_origin) 1361 1464 1362 assert num.allclose(new_points, ab_points) 1465 1363 1466 1364 geo = Geo_reference(56,67,-56) 1467 1468 data_points = geo.change_points_geo_ref(ab_points) 1469 new_points = ensure_absolute(data_points, 1470 geo_reference=geo) 1471 #print "new_points",new_points 1472 #print "ab_points",ab_points 1473 1365 data_points = geo.change_points_geo_ref(ab_points) 1366 new_points = ensure_absolute(data_points, geo_reference=geo) 1367 1474 1368 assert num.allclose(new_points, ab_points) 1475 1476 1369 1477 1370 geo_reference = Geo_reference(56, 100, 200) … … 1479 1372 points = geo_reference.change_points_geo_ref(ab_points) 1480 1373 attributes = [2, 4] 1481 #print "geo in points", points 1482 G = Geospatial_data(points, attributes, 1483 geo_reference=geo_reference) 1484 1374 G = Geospatial_data(points, attributes, geo_reference=geo_reference) 1485 1375 new_points = ensure_absolute(G) 1486 #print "new_points",new_points 1487 #print "ab_points",ab_points 1488 1376 1489 1377 assert num.allclose(new_points, ab_points) 1490 1378 1491 1492 1493 1379 def test_ensure_geospatial(self): 1494 points = [[2.0, 0.0], [1.0, 1.0],1495 [2.0, 0.0],[2.0, 2.0],1496 [1.0, 3.0],[2.0, 2.0]]1380 points = [[2.0, 0.0], [1.0, 1.0], 1381 [2.0, 0.0], [2.0, 2.0], 1382 [1.0, 3.0], [2.0, 2.0]] 1497 1383 new_points = ensure_geospatial(points) 1498 1499 assert num.allclose(new_points.get_data_points(absolute = True), points) 1500 1501 points = num.array([[2.0, 0.0],[1.0, 1.0], 1502 [2.0, 0.0],[2.0, 2.0], 1503 [1.0, 3.0],[2.0, 2.0]]) 1384 assert num.allclose(new_points.get_data_points(absolute=True), points) 1385 1386 points = num.array([[2.0, 0.0], [1.0, 1.0], 1387 [2.0, 0.0], [2.0, 2.0], 1388 [1.0, 3.0], [2.0, 2.0]]) 1504 1389 new_points = ensure_geospatial(points) 1505 1506 assert num.allclose(new_points.get_data_points(absolute = True), points) 1507 1390 assert num.allclose(new_points.get_data_points(absolute=True), points) 1391 1508 1392 ab_points = num.array([[2.0, 0.0],[1.0, 1.0], 1509 1393 [2.0, 0.0],[2.0, 2.0], 1510 1394 [1.0, 3.0],[2.0, 2.0]]) 1511 1512 mesh_origin = (56, 290000, 618000) #zone, easting, northing 1513 1514 data_points = num.zeros((ab_points.shape), num.Float) 1395 mesh_origin = (56, 290000, 618000) # zone, easting, northing 1396 data_points = num.zeros((ab_points.shape), num.float) 1397 1515 1398 #Shift datapoints according to new origins 1516 1399 for k in range(len(ab_points)): 1517 1400 data_points[k][0] = ab_points[k][0] - mesh_origin[1] 1518 1401 data_points[k][1] = ab_points[k][1] - mesh_origin[2] 1519 #print "data_points",data_points1520 1402 new_geospatial = ensure_geospatial(data_points, 1521 1403 geo_reference=mesh_origin) 1522 1404 new_points = new_geospatial.get_data_points(absolute=True) 1523 #print "new_points",new_points1524 #print "ab_points",ab_points1525 1526 1405 assert num.allclose(new_points, ab_points) 1527 1406 1528 geo = Geo_reference(56,67,-56) 1529 1530 data_points = geo.change_points_geo_ref(ab_points) 1531 new_geospatial = ensure_geospatial(data_points, 1532 geo_reference=geo) 1407 geo = Geo_reference(56, 67, -56) 1408 data_points = geo.change_points_geo_ref(ab_points) 1409 new_geospatial = ensure_geospatial(data_points, geo_reference=geo) 1533 1410 new_points = new_geospatial.get_data_points(absolute=True) 1534 #print "new_points",new_points1535 #print "ab_points",ab_points1536 1537 1411 assert num.allclose(new_points, ab_points) 1538 1539 1412 1540 1413 geo_reference = Geo_reference(56, 100, 200) … … 1542 1415 points = geo_reference.change_points_geo_ref(ab_points) 1543 1416 attributes = [2, 4] 1544 #print "geo in points", points 1545 G = Geospatial_data(points, attributes, 1546 geo_reference=geo_reference) 1547 1548 new_geospatial = ensure_geospatial(G) 1417 G = Geospatial_data(points, attributes, geo_reference=geo_reference) 1418 new_geospatial = ensure_geospatial(G) 1549 1419 new_points = new_geospatial.get_data_points(absolute=True) 1550 #print "new_points",new_points1551 #print "ab_points",ab_points1552 1553 1420 assert num.allclose(new_points, ab_points) 1554 1421 1555 1422 def test_isinstance(self): 1556 1557 import os 1558 1559 fileName = tempfile.mktemp(".csv") 1560 file = open(fileName,"w") 1561 file.write("x,y, elevation , speed \n\ 1562 1.0, 0.0, 10.0, 0.0\n\ 1563 0.0, 1.0, 0.0, 10.0\n\ 1564 1.0, 0.0, 10.4, 40.0\n") 1423 fileName = tempfile.mktemp('.csv') 1424 file = open(fileName, 'w') 1425 file.write('x,y, elevation , speed \n\ 1426 1.0, 0.0, 10.0, 0.0\n\ 1427 0.0, 1.0, 0.0, 10.0\n\ 1428 1.0, 0.0, 10.4, 40.0\n') 1565 1429 file.close() 1566 1430 1567 1431 results = Geospatial_data(fileName) 1568 1432 assert num.allclose(results.get_data_points(absolute=True), 1569 [[1.0, 0.0], [0.0, 1.0],[1.0, 0.0]])1433 [[1.0, 0.0], [0.0, 1.0], [1.0, 0.0]]) 1570 1434 assert num.allclose(results.get_attributes(attribute_name='elevation'), 1571 1435 [10.0, 0.0, 10.4]) … … 1574 1438 1575 1439 os.remove(fileName) 1576 1577 1440 1578 1441 def test_no_constructors(self): 1579 1580 1442 try: 1581 1443 G = Geospatial_data() 1582 # results = Geospatial_data(file_name = fileName)1583 # dict = import_points_file(fileName)1584 1444 except ValueError: 1585 1445 pass 1586 1446 else: 1587 1447 msg = 'Instance must have a filename or data points' 1588 raise msg1448 raise Exception, msg 1589 1449 1590 1450 def test_load_csv_lat_long(self): 1591 """ 1592 comma delimited 1593 1594 """ 1595 fileName = tempfile.mktemp(".csv") 1596 file = open(fileName,"w") 1597 file.write("long,lat, elevation, yeah \n\ 1451 '''comma delimited''' 1452 1453 fileName = tempfile.mktemp('.csv') 1454 file = open(fileName, 'w') 1455 file.write('long,lat, elevation, yeah \n\ 1598 1456 150.916666667,-34.50,452.688000, 10\n\ 1599 150.0,-34,459.126000, 10\n ")1457 150.0,-34,459.126000, 10\n') 1600 1458 file.close() 1459 1601 1460 results = Geospatial_data(fileName) 1602 1461 os.remove(fileName) 1603 1462 points = results.get_data_points() 1604 1463 1605 1464 assert num.allclose(points[0][0], 308728.009) 1606 1465 assert num.allclose(points[0][1], 6180432.601) 1607 1466 assert num.allclose(points[1][0], 222908.705) 1608 1467 assert num.allclose(points[1][1], 6233785.284) 1609 1610 1468 1611 1469 def test_load_csv_lat_longII(self): 1612 """ 1613 comma delimited 1614 1615 """ 1616 fileName = tempfile.mktemp(".csv") 1617 file = open(fileName,"w") 1618 file.write("Lati,LONG,z \n\ 1470 '''comma delimited''' 1471 1472 fileName = tempfile.mktemp('.csv') 1473 file = open(fileName, 'w') 1474 file.write('Lati,LONG,z \n\ 1619 1475 -34.50,150.916666667,452.688000\n\ 1620 -34,150.0,459.126000\n ")1476 -34,150.0,459.126000\n') 1621 1477 file.close() 1478 1622 1479 results = Geospatial_data(fileName) 1623 1480 os.remove(fileName) 1624 1481 points = results.get_data_points() 1625 1482 1626 1483 assert num.allclose(points[0][0], 308728.009) 1627 1484 assert num.allclose(points[0][1], 6180432.601) 1628 assert num.allclose(points[1][0], 1485 assert num.allclose(points[1][0], 222908.705) 1629 1486 assert num.allclose(points[1][1], 6233785.284) 1630 1487 1631 1632 1488 def test_load_csv_lat_long_bad(self): 1633 """ 1634 comma delimited 1635 1636 """ 1637 fileName = tempfile.mktemp(".csv") 1638 file = open(fileName,"w") 1639 file.write("Lati,LONG,z \n\ 1489 '''comma delimited''' 1490 1491 fileName = tempfile.mktemp('.csv') 1492 file = open(fileName, 'w') 1493 file.write('Lati,LONG,z \n\ 1640 1494 -25.0,180.0,452.688000\n\ 1641 -34,150.0,459.126000\n ")1495 -34,150.0,459.126000\n') 1642 1496 file.close() 1497 1643 1498 try: 1644 1499 results = Geospatial_data(fileName) … … 1647 1502 else: 1648 1503 msg = 'Different zones in Geo references not caught.' 1649 raise msg1650 1651 os.remove(fileName) 1652 1504 raise Exception, msg 1505 1506 os.remove(fileName) 1507 1653 1508 def test_lat_long(self): 1654 lat_gong = degminsec2decimal_degrees(-34, 30,0.)1655 lon_gong = degminsec2decimal_degrees(150, 55,0.)1656 1657 lat_2 = degminsec2decimal_degrees(-34, 00,0.)1658 lon_2 = degminsec2decimal_degrees(150, 00,0.)1659 1509 lat_gong = degminsec2decimal_degrees(-34, 30, 0.) 1510 lon_gong = degminsec2decimal_degrees(150, 55, 0.) 1511 1512 lat_2 = degminsec2decimal_degrees(-34, 00, 0.) 1513 lon_2 = degminsec2decimal_degrees(150, 00, 0.) 1514 1660 1515 lats = [lat_gong, lat_2] 1661 1516 longs = [lon_gong, lon_2] … … 1663 1518 1664 1519 points = gsd.get_data_points(absolute=True) 1665 1520 1666 1521 assert num.allclose(points[0][0], 308728.009) 1667 1522 assert num.allclose(points[0][1], 6180432.601) 1668 assert num.allclose(points[1][0], 1523 assert num.allclose(points[1][0], 222908.705) 1669 1524 assert num.allclose(points[1][1], 6233785.284) 1670 1525 self.failUnless(gsd.get_geo_reference().get_zone() == 56, 1671 1526 'Bad zone error!') 1672 1527 1673 1528 try: 1674 1529 results = Geospatial_data(latitudes=lats) … … 1676 1531 pass 1677 1532 else: 1678 self.fail Unless(0 ==1,'Error not thrown error!')1533 self.fail('Error not thrown error!') 1679 1534 try: 1680 1535 results = Geospatial_data(latitudes=lats) … … 1682 1537 pass 1683 1538 else: 1684 self.fail Unless(0 ==1,'Error not thrown error!')1539 self.fail('Error not thrown error!') 1685 1540 try: 1686 1541 results = Geospatial_data(longitudes=lats) … … 1688 1543 pass 1689 1544 else: 1690 self.fail Unless(0 ==1,'Error not thrown error!')1545 self.fail('Error not thrown error!') 1691 1546 try: 1692 1547 results = Geospatial_data(latitudes=lats, longitudes=longs, … … 1695 1550 pass 1696 1551 else: 1697 self.fail Unless(0 ==1,'Error not thrown error!')1698 1552 self.fail('Error not thrown error!') 1553 1699 1554 try: 1700 1555 results = Geospatial_data(latitudes=lats, longitudes=longs, … … 1703 1558 pass 1704 1559 else: 1705 self.fail Unless(0 ==1,'Error not thrown error!')1560 self.fail('Error not thrown error!') 1706 1561 1707 1562 def test_lat_long2(self): 1708 lat_gong = degminsec2decimal_degrees(-34, 30,0.)1709 lon_gong = degminsec2decimal_degrees(150, 55,0.)1710 1711 lat_2 = degminsec2decimal_degrees(-34, 00,0.)1712 lon_2 = degminsec2decimal_degrees(150, 00,0.)1713 1563 lat_gong = degminsec2decimal_degrees(-34, 30, 0.) 1564 lon_gong = degminsec2decimal_degrees(150, 55, 0.) 1565 1566 lat_2 = degminsec2decimal_degrees(-34, 00, 0.) 1567 lon_2 = degminsec2decimal_degrees(150, 00, 0.) 1568 1714 1569 points = [[lat_gong, lon_gong], [lat_2, lon_2]] 1715 1570 gsd = Geospatial_data(data_points=points, points_are_lats_longs=True) 1716 1571 1717 1572 points = gsd.get_data_points(absolute=True) 1718 1573 1719 1574 assert num.allclose(points[0][0], 308728.009) 1720 1575 assert num.allclose(points[0][1], 6180432.601) 1721 assert num.allclose(points[1][0], 1576 assert num.allclose(points[1][0], 222908.705) 1722 1577 assert num.allclose(points[1][1], 6233785.284) 1723 1578 self.failUnless(gsd.get_geo_reference().get_zone() == 56, … … 1729 1584 pass 1730 1585 else: 1731 self.failUnless(0 ==1, 'Error not thrown error!') 1732 1586 self.fail('Error not thrown error!') 1733 1587 1734 1588 def test_write_urs_file(self): 1735 lat_gong = degminsec2decimal_degrees(-34, 00,0)1736 lon_gong = degminsec2decimal_degrees(150, 30,0.)1737 1738 lat_2 = degminsec2decimal_degrees(-34, 00,1)1739 lon_2 = degminsec2decimal_degrees(150, 00,0.)1589 lat_gong = degminsec2decimal_degrees(-34, 00, 0) 1590 lon_gong = degminsec2decimal_degrees(150, 30, 0.) 1591 1592 lat_2 = degminsec2decimal_degrees(-34, 00, 1) 1593 lon_2 = degminsec2decimal_degrees(150, 00, 0.) 1740 1594 p1 = (lat_gong, lon_gong) 1741 1595 p2 = (lat_2, lon_2) … … 1743 1597 gsd = Geospatial_data(data_points=list(points), 1744 1598 points_are_lats_longs=True) 1745 1746 fn = tempfile.mktemp( ".urs")1599 1600 fn = tempfile.mktemp('.urs') 1747 1601 gsd.export_points_file(fn) 1748 #print "fn", fn1749 1602 handle = open(fn) 1750 1603 lines = handle.readlines() 1751 assert lines[0], '2'1752 assert lines[1], '-34.0002778 150.0 0'1753 assert lines[2], '-34.0 150.5 1'1604 assert lines[0], '2' 1605 assert lines[1], '-34.0002778 150.0 0' 1606 assert lines[2], '-34.0 150.5 1' 1754 1607 handle.close() 1755 1608 os.remove(fn) 1756 1609 1757 1610 def test_lat_long_set(self): 1758 lat_gong = degminsec2decimal_degrees(-34, 30,0.)1759 lon_gong = degminsec2decimal_degrees(150, 55,0.)1760 1761 lat_2 = degminsec2decimal_degrees(-34, 00,0.)1762 lon_2 = degminsec2decimal_degrees(150, 00,0.)1611 lat_gong = degminsec2decimal_degrees(-34, 30, 0.) 1612 lon_gong = degminsec2decimal_degrees(150, 55, 0.) 1613 1614 lat_2 = degminsec2decimal_degrees(-34, 00, 0.) 1615 lon_2 = degminsec2decimal_degrees(150, 00, 0.) 1763 1616 p1 = (lat_gong, lon_gong) 1764 1617 p2 = (lat_2, lon_2) … … 1768 1621 1769 1622 points = gsd.get_data_points(absolute=True) 1770 #print "points[0][0]", points[0][0]1771 1623 #Note the order is unknown, due to using sets 1772 1624 # and it changes from windows to linux … … 1774 1626 assert num.allclose(points[1][0], 308728.009) 1775 1627 assert num.allclose(points[1][1], 6180432.601) 1776 assert num.allclose(points[0][0], 1628 assert num.allclose(points[0][0], 222908.705) 1777 1629 assert num.allclose(points[0][1], 6233785.284) 1778 1630 except AssertionError: 1779 1631 assert num.allclose(points[0][0], 308728.009) 1780 1632 assert num.allclose(points[0][1], 6180432.601) 1781 assert num.allclose(points[1][0], 1633 assert num.allclose(points[1][0], 222908.705) 1782 1634 assert num.allclose(points[1][1], 6233785.284) 1783 1635 1784 1636 self.failUnless(gsd.get_geo_reference().get_zone() == 56, 1785 1637 'Bad zone error!') 1786 1638 points = gsd.get_data_points(as_lat_long=True) 1787 #print "test_lat_long_set points", points1788 1639 try: 1789 1640 assert num.allclose(points[0][0], -34) … … 1794 1645 1795 1646 def test_len(self): 1796 1797 1647 points = [[1.0, 2.1], [3.0, 5.3]] 1798 1648 G = Geospatial_data(points) 1799 self.failUnless(2 == len(G),'Len error!')1800 1649 self.failUnless(2 == len(G), 'Len error!') 1650 1801 1651 points = [[1.0, 2.1]] 1802 1652 G = Geospatial_data(points) 1803 self.failUnless(1 == len(G),'Len error!')1653 self.failUnless(1 == len(G), 'Len error!') 1804 1654 1805 1655 points = [[1.0, 2.1], [3.0, 5.3], [3.0, 5.3], [3.0, 5.3]] 1806 1656 G = Geospatial_data(points) 1807 self.failUnless(4 == len(G),'Len error!')1808 1657 self.failUnless(4 == len(G), 'Len error!') 1658 1809 1659 def test_split(self): 1810 1660 """test if the results from spilt are disjoin sets""" 1811 1812 #below is a work around until the randint works on cyclones compute nodes 1813 if get_host_name()[8:9]!='0': 1814 1815 1661 1662 # below is a workaround until randint works on cyclones compute nodes 1663 if get_host_name()[8:9] != '0': 1816 1664 points = [[1.0, 1.0], [1.0, 2.0],[1.0, 3.0], [1.0, 4.0], [1.0, 5.0], 1817 1665 [2.0, 1.0], [2.0, 2.0],[2.0, 3.0], [2.0, 4.0], [2.0, 5.0], … … 1819 1667 [4.0, 1.0], [4.0, 2.0],[4.0, 3.0], [4.0, 4.0], [4.0, 5.0], 1820 1668 [5.0, 1.0], [5.0, 2.0],[5.0, 3.0], [5.0, 4.0], [5.0, 5.0]] 1821 attributes = {'depth':[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1822 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25], 1823 'speed':[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1824 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]} 1669 attributes = {'depth': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1670 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1671 21, 22, 23, 24, 25], 1672 'speed': [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 1673 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 1674 21, 22, 23, 24, 25]} 1825 1675 G = Geospatial_data(points, attributes) 1826 1676 1827 1677 factor = 0.21 1828 1829 #will return G1 with 10% of points and G2 with 90% 1830 G1, G2 = G.split(factor,100) 1831 1678 1679 # will return G1 with 10% of points and G2 with 90% 1680 G1, G2 = G.split(factor, 100) 1832 1681 assert num.allclose(len(G), len(G1)+len(G2)) 1833 1682 assert num.allclose(round(len(G)*factor), len(G1)) 1834 1683 1835 1684 P = G1.get_data_points(absolute=False) 1836 assert num.allclose(P, [[5.0,4.0],[4.0,3.0],[4.0,2.0],[3.0,1.0],[2.0,3.0]]) 1837 1685 expected = [[5.0, 4.0], [4.0, 3.0], [4.0, 2.0], 1686 [3.0, 1.0], [2.0, 3.0]] 1687 msg = 'Expected %s, but\nP=%s' % (str(expected), str(P)) 1688 assert num.allclose(P, expected), msg 1689 1838 1690 A = G1.get_attributes() 1839 assert num.allclose(A,[24, 18, 17, 11, 8]) 1840 1691 expected = [24, 18, 17, 11, 8] 1692 msg = 'expected=%s, but A=%s' % (str(expected), str(A)) 1693 assert num.allclose(A, expected), msg 1694 1841 1695 def test_split1(self): 1842 """test if the results from sp ilt are disjoinsets"""1843 #below is a work around until the randint works on cyclones compute nodes 1844 if get_host_name()[8:9]!='0':1845 1696 """test if the results from split are disjoint sets""" 1697 1698 # below is a workaround until randint works on cyclones compute nodes 1699 if get_host_name()[8:9] != '0': 1846 1700 from RandomArray import randint,seed 1701 1847 1702 seed(100,100) 1848 a_points = randint(0, 999999,(10,2))1703 a_points = randint(0, 999999, (10,2)) 1849 1704 points = a_points.tolist() 1850 # print points 1851 1705 1852 1706 G = Geospatial_data(points) 1853 1707 1854 1708 factor = 0.1 1855 1856 #will return G1 with 10% of points and G2 with 90% 1857 G1, G2 = G.split(factor,100) 1858 1859 # print 'G1',G1 1709 1710 # will return G1 with 10% of points and G2 with 90% 1711 G1, G2 = G.split(factor, 100) 1860 1712 assert num.allclose(len(G), len(G1)+len(G2)) 1861 1713 assert num.allclose(round(len(G)*factor), len(G1)) 1862 1714 1863 1715 P = G1.get_data_points(absolute=False) 1864 assert num.allclose(P, [[982420.,28233.]]) 1865 1866 1716 expected = [[982420., 28233.]] 1717 msg = 'expected=%s, but\nP=%s' % (str(expected), str(P)) 1718 assert num.allclose(P, expected), msg 1719 1867 1720 def test_find_optimal_smoothing_parameter(self): 1868 1721 """ 1869 Creates a elevation file repres ting hill (sort of) and runs1722 Creates a elevation file representing a hill (sort of) and runs 1870 1723 find_optimal_smoothing_parameter for 3 different alphas, 1871 1724 1872 1725 NOTE the random number seed is provided to control the results 1873 1726 """ 1727 1874 1728 from cmath import cos 1875 1729 1876 # below is a work around until therandint works on cyclones compute nodes1730 # below is a workaround until randint works on cyclones compute nodes 1877 1731 if get_host_name()[8:9]!='0': 1878 1879 filename = tempfile.mktemp(".csv") 1880 file = open(filename,"w") 1881 file.write("x,y,elevation \n") 1882 1883 for i in range(-5,6): 1884 for j in range(-5,6): 1885 #this equation made surface like a circle ripple 1732 filename = tempfile.mktemp('.csv') 1733 file = open(filename, 'w') 1734 file.write('x,y,elevation \n') 1735 1736 for i in range(-5, 6): 1737 for j in range(-5, 6): 1738 # this equation makes a surface like a circle ripple 1886 1739 z = abs(cos(((i*i) + (j*j))*.1)*2) 1887 # print 'x,y,f',i,j,z 1888 file.write("%s, %s, %s\n" %(i, j, z)) 1889 1740 file.write("%s, %s, %s\n" % (i, j, z)) 1741 1890 1742 file.close() 1891 1892 value, alpha = find_optimal_smoothing_parameter(data_file=filename, 1743 1744 value, alpha = find_optimal_smoothing_parameter(data_file=filename, 1893 1745 alpha_list=[0.0001, 0.01, 1], 1894 1746 mesh_file=None, … … 1901 1753 seed_num=100000, 1902 1754 verbose=False) 1903 1755 1904 1756 os.remove(filename) 1905 1757 1906 1758 # print value, alpha 1907 assert (alpha==0.01)1759 assert(alpha == 0.01) 1908 1760 1909 1761 def test_find_optimal_smoothing_parameter1(self): 1910 1762 """ 1911 Creates a elevation file repres tinghill (sort of) and1763 Creates a elevation file representing a hill (sort of) and 1912 1764 Then creates a mesh file and passes the mesh file and the elevation 1913 1765 file to find_optimal_smoothing_parameter for 3 different alphas, 1914 1766 1915 1767 NOTE the random number seed is provided to control the results 1916 1768 """ 1917 #below is a work around until the randint works on cyclones compute nodes 1769 1770 # below is a workaround until randint works on cyclones compute nodes 1918 1771 if get_host_name()[8:9]!='0': 1919 1920 1772 from cmath import cos 1921 1773 from anuga.pmesh.mesh_interface import create_mesh_from_regions 1922 1923 filename = tempfile.mktemp( ".csv")1924 file = open(filename, "w")1925 file.write( "x,y,elevation \n")1926 1927 for i in range(-5 ,6):1928 for j in range(-5, 6):1929 # this equation madesurface like a circle ripple1774 1775 filename = tempfile.mktemp('.csv') 1776 file = open(filename, 'w') 1777 file.write('x,y,elevation \n') 1778 1779 for i in range(-5 ,6): 1780 for j in range(-5, 6): 1781 # this equation makes a surface like a circle ripple 1930 1782 z = abs(cos(((i*i) + (j*j))*.1)*2) 1931 # print 'x,y,f',i,j,z 1932 file.write("%s, %s, %s\n" %(i, j, z)) 1933 1783 file.write('%s, %s, %s\n' % (i, j, z)) 1784 1934 1785 file.close() 1935 poly=[[5,5],[5,-5],[-5,-5],[-5,5]] 1936 internal_poly=[[[[1,1],[1,-1],[-1,-1],[-1,1]],.5]] 1937 mesh_filename= tempfile.mktemp(".msh") 1938 1786 1787 poly=[[5,5], [5,-5], [-5,-5], [-5,5]] 1788 internal_poly=[[[[1,1], [1,-1], [-1,-1], [-1,1]], .5]] 1789 mesh_filename= tempfile.mktemp('.msh') 1790 1939 1791 create_mesh_from_regions(poly, 1940 boundary_tags={'back': [2],1941 'side': [1,3],1942 'ocean': [0]},1943 maximum_triangle_area=3,1944 interior_regions=internal_poly,1945 filename=mesh_filename,1946 use_cache=False,1947 verbose=False)1948 1949 value, alpha = find_optimal_smoothing_parameter(data_file=filename, 1792 boundary_tags={'back': [2], 1793 'side': [1,3], 1794 'ocean': [0]}, 1795 maximum_triangle_area=3, 1796 interior_regions=internal_poly, 1797 filename=mesh_filename, 1798 use_cache=False, 1799 verbose=False) 1800 1801 value, alpha = find_optimal_smoothing_parameter(data_file=filename, 1950 1802 alpha_list=[0.0001, 0.01, 1], 1951 1803 mesh_file=mesh_filename, … … 1953 1805 seed_num=174, 1954 1806 verbose=False) 1955 1807 1956 1808 os.remove(filename) 1957 1809 os.remove(mesh_filename) 1958 1959 # print value, alpha 1960 assert (alpha==0.01) 1810 1811 assert(alpha == 0.01) 1961 1812 1962 1813 def test_find_optimal_smoothing_parameter2(self): 1963 """ 1964 Tests requirement that mesh file must exist or IOError is thrown 1965 1814 '''Tests requirement that mesh file must exist or IOError is thrown 1815 1966 1816 NOTE the random number seed is provided to control the results 1967 """ 1817 ''' 1818 1968 1819 from cmath import cos 1969 1820 from anuga.pmesh.mesh_interface import create_mesh_from_regions 1970 1971 filename = tempfile.mktemp( ".csv")1972 mesh_filename= tempfile.mktemp( ".msh")1973 1821 1822 filename = tempfile.mktemp('.csv') 1823 mesh_filename= tempfile.mktemp('.msh') 1824 1974 1825 try: 1975 value, alpha = find_optimal_smoothing_parameter(data_file=filename, 1826 value, alpha = find_optimal_smoothing_parameter(data_file=filename, 1976 1827 alpha_list=[0.0001, 0.01, 1], 1977 1828 mesh_file=mesh_filename, … … 1982 1833 pass 1983 1834 else: 1984 self.failUnless(0 ==1, 'Error not thrown error!') 1985 1986 1835 self.fail('Error not thrown error!') 1836 1837 ################################################################################ 1838 # Test the clean_line() utility function. 1839 ################################################################################ 1840 1841 # helper routine to test clean_line() 1842 def clean_line_test(self, instr, delim, expected): 1843 result = clean_line(instr, delim) 1844 self.failUnless(result == expected, 1845 "clean_line('%s', '%s'), expected %s, got %s" 1846 % (str(instr), str(delim), str(expected), str(result))) 1847 1848 def test_clean_line_01(self): 1849 self.clean_line_test('abc, ,,xyz,123', ',', ['abc', '', 'xyz', '123']) 1850 1851 def test_clean_line_02(self): 1852 self.clean_line_test(' abc , ,, xyz , 123 ', ',', 1853 ['abc', '', 'xyz', '123']) 1854 1855 def test_clean_line_03(self): 1856 self.clean_line_test('1||||2', '|', ['1', '2']) 1857 1858 def test_clean_line_04(self): 1859 self.clean_line_test('abc, ,,xyz,123, ', ',', 1860 ['abc', '', 'xyz', '123', '']) 1861 1862 def test_clean_line_05(self): 1863 self.clean_line_test('abc, ,,xyz,123, , ', ',', 1864 ['abc', '', 'xyz', '123', '', '']) 1865 1866 def test_clean_line_06(self): 1867 self.clean_line_test(',,abc, ,,xyz,123, , ', ',', 1868 ['abc', '', 'xyz', '123', '', '']) 1869 1870 def test_clean_line_07(self): 1871 self.clean_line_test('|1||||2', '|', ['1', '2']) 1872 1873 def test_clean_line_08(self): 1874 self.clean_line_test(' ,a,, , ,b,c , ,, , ', ',', 1875 ['', 'a', '', '', 'b', 'c', '', '', '']) 1876 1877 def test_clean_line_09(self): 1878 self.clean_line_test('a:b:c', ':', ['a', 'b', 'c']) 1879 1880 def test_clean_line_10(self): 1881 self.clean_line_test('a:b:c:', ':', ['a', 'b', 'c']) 1882 1883 # new version of function should leave last field if contains '\n'. 1884 # cf. test_clean_line_10() above. 1885 def test_clean_line_11(self): 1886 self.clean_line_test('a:b:c:\n', ':', ['a', 'b', 'c', '']) 1887 1987 1888 if __name__ == "__main__": 1988 1989 #suite = unittest.makeSuite(Test_Geospatial_data, 'test_write_csv_attributes_lat_long')1990 #suite = unittest.makeSuite(Test_Geospatial_data, 'test_find_optimal_smoothing_parameter')1991 #suite = unittest.makeSuite(Test_Geospatial_data, 'test_split1')1992 1889 suite = unittest.makeSuite(Test_Geospatial_data, 'test') 1993 1890 runner = unittest.TextTestRunner() #verbosity=2) 1994 1891 runner.run(suite) 1995 1892 1996 1893
Note: See TracChangeset
for help on using the changeset viewer.