Changeset 2262
- Timestamp:
- Jan 20, 2006, 11:04:18 AM (19 years ago)
- Location:
- inundation
- Files:
-
- 7 edited
-
coordinate_transforms/geo_reference.py (modified) (1 diff)
-
load_mesh/loadASCII.py (modified) (1 diff)
-
pyvolution/data_manager.py (modified) (1 diff)
-
pyvolution/quantity.py (modified) (12 diffs)
-
pyvolution/test_data_manager.py (modified) (1 diff)
-
pyvolution/test_least_squares.py (modified) (4 diffs)
-
pyvolution/test_quantity.py (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
inundation/coordinate_transforms/geo_reference.py
r2261 r2262 21 21 yllcorner = 0.0, 22 22 datum = 'wgs84', 23 projection = 'UTM', 24 units = 'm', 23 projection = 'UTM', units = 'm', 25 24 false_easting = 500000, 26 25 false_northing = 10000000, #Default for southern hemisphere -
inundation/load_mesh/loadASCII.py
r2257 r2262 1170 1170 # There can't be an attribute called points 1171 1171 # consider format change 1172 1173 1174 legal_keys = ['pointlist', 'attributelist', 'geo_reference'] 1175 for key in point_atts.keys(): 1176 msg = 'Key %s is illegal. Valid keys are %s' %(key, legal_keys) 1177 assert key in legal_keys, msg 1172 1178 1173 1179 from Scientific.IO.NetCDF import NetCDFFile -
inundation/pyvolution/data_manager.py
r2256 r2262 1023 1023 """Write points and associated attribute to pts (NetCDF) format 1024 1024 """ 1025 1026 print 'WARNING: write_ptsfile is obsolete. Use export_points from load_mesh.loadASCII instead.' 1025 1027 1026 1028 from Numeric import Float -
inundation/pyvolution/quantity.py
r2233 r2262 167 167 quantity = None, # Another quantity 168 168 function = None, # Callable object: f(x,y) 169 points = None, values = None, #Input for least squares169 points = None, values = None, data_georef = None, #Input for least squares 170 170 filename = None, attribute_name = None, #Input from file 171 171 alpha = None, … … 201 201 If points is specified, values is an array of length N containing 202 202 attribute values for each point. 203 204 data_georef: 205 If points is specified, geo_reference applies to each point. 203 206 204 207 filename: … … 305 308 assert values is not None, msg 306 309 self.set_values_from_points(points, values, alpha, 307 location, indices, verbose, 308 use_cache) 310 location, indices, 311 data_georef = data_georef, 312 verbose = verbose, 313 use_cache = use_cache) 309 314 elif filename is not None: 310 315 self.set_values_from_file(filename, attribute_name, alpha, 311 location, indices, verbose, 312 use_cache) 316 location, indices, 317 verbose = verbose, 318 use_cache = use_cache) 313 319 else: 314 320 raise 'This can\'t happen :-)' … … 513 519 #FIXME: Should check that function returns something sensible and 514 520 #raise a meaningfull exception if it returns None for example 521 522 #FIXME: Should supply absolute coordinates 515 523 516 524 from Numeric import take … … 546 554 547 555 def set_values_from_points(self, points, values, alpha, 548 location, indices, verbose, use_cache): 556 location, indices, 557 data_georef = None, 558 verbose = False, 559 use_cache = False): 549 560 """Set quantity values from arbitray data points using least squares 550 561 """ … … 553 564 from util import ensure_numeric 554 565 from least_squares import fit_to_mesh 566 from coordinate_transforms.geo_reference import Geo_reference 567 555 568 556 569 points = ensure_numeric(points, Float) … … 565 578 triangles = self.domain.triangles 566 579 580 581 #Take care of georeferencing 582 if data_georef is None: 583 data_georef = Geo_reference() 584 585 586 mesh_georef = self.domain.geo_reference 587 588 #print mesh_georef 589 #print data_georef 590 #print points 591 592 593 #Call least squares method 594 args = (coordinates, triangles, points, values) 595 kwargs = {'data_origin': data_georef.get_origin(), 596 'mesh_origin': mesh_georef.get_origin(), 597 'alpha': alpha, 598 'verbose': verbose} 599 600 #print kwargs 601 567 602 if use_cache is True: 568 603 try: … … 573 608 raise msg 574 609 575 args = (coordinates, triangles, points, values)576 kwargs = {'alpha': alpha, 'verbose': verbose}577 610 vertex_attributes = cache(fit_to_mesh, 578 611 args, kwargs, 579 612 verbose = verbose) 580 613 else: 581 vertex_attributes = fit_to_mesh(coordinates, 582 triangles, 583 points, 584 values, 585 alpha = alpha, 586 verbose = verbose) 587 588 614 615 vertex_attributes = apply(fit_to_mesh, 616 args, kwargs) 617 618 #Call underlying method using array values 589 619 self.set_values_from_array(vertex_attributes, 590 620 location, indices, verbose) … … 595 625 596 626 def set_values_from_file(self, filename, attribute_name, alpha, 597 location, indices, verbose, use_cache): 627 location, indices, 628 verbose = False, 629 use_cache = False): 598 630 """Set quantity based on arbitrary points in .pts file 599 631 using least_squares attribute_name selects name of attribute … … 602 634 """ 603 635 636 from load_mesh.loadASCII import import_points_file 604 637 605 638 from types import StringType … … 609 642 610 643 #Read from (NetCDF) file 611 from load_mesh.loadASCII import import_points_file612 644 points_dict = import_points_file(filename) 613 645 points = points_dict['pointlist'] … … 634 666 635 667 636 #Call least squares method 668 #Take care of georeferencing 669 if points_dict.has_key('geo_reference') and \ 670 points_dict['geo_reference'] is not None: 671 data_georef = points_dict['geo_reference'] 672 else: 673 data_georef = None 674 675 676 #Call underlying method for points 637 677 self.set_values_from_points(points, z, alpha, 638 location, indices, verbose, use_cache) 678 location, indices, 679 data_georef = data_georef, 680 verbose = verbose, 681 use_cache = use_cache) 639 682 640 683 -
inundation/pyvolution/test_data_manager.py
r2074 r2262 439 439 440 440 441 def test_write_pts(self): 442 #Get (enough) datapoints 443 444 from Numeric import array 445 points = array([[ 0.66666667, 0.66666667], 446 [ 1.33333333, 1.33333333], 447 [ 2.66666667, 0.66666667], 448 [ 0.66666667, 2.66666667], 449 [ 0.0, 1.0], 450 [ 0.0, 3.0], 451 [ 1.0, 0.0], 452 [ 1.0, 1.0], 453 [ 1.0, 2.0], 454 [ 1.0, 3.0], 455 [ 2.0, 1.0], 456 [ 3.0, 0.0], 457 [ 3.0, 1.0]]) 458 459 z = points[:,0] + 2*points[:,1] 460 461 ptsfile = 'testptsfile.pts' 462 write_ptsfile(ptsfile, points, z, 463 attribute_name = 'linear_combination') 464 465 #Check contents 466 #Get NetCDF 467 from Scientific.IO.NetCDF import NetCDFFile 468 fid = NetCDFFile(ptsfile, 'r') 469 470 # Get the variables 471 #print fid.variables.keys() 472 points1 = fid.variables['points'] 473 z1 = fid.variables['linear_combination'] 474 475 #Check values 476 477 #print points[:] 478 #print ref_points 479 assert allclose(points, points1) 480 481 #print attributes[:] 482 #print ref_elevation 483 assert allclose(z, z1) 484 485 #Cleanup 486 fid.close() 487 488 import os 489 os.remove(ptsfile) 441 #def test_write_pts(self): 442 # #Obsolete 443 # 444 # #Get (enough) datapoints 445 # 446 # from Numeric import array 447 # points = array([[ 0.66666667, 0.66666667], 448 # [ 1.33333333, 1.33333333], 449 # [ 2.66666667, 0.66666667], 450 # [ 0.66666667, 2.66666667], 451 # [ 0.0, 1.0], 452 # [ 0.0, 3.0], 453 # [ 1.0, 0.0], 454 # [ 1.0, 1.0], 455 # [ 1.0, 2.0], 456 # [ 1.0, 3.0], 457 # [ 2.0, 1.0], 458 # [ 3.0, 0.0], 459 # [ 3.0, 1.0]]) 460 # 461 # z = points[:,0] + 2*points[:,1] 462 # 463 # ptsfile = 'testptsfile.pts' 464 # write_ptsfile(ptsfile, points, z, 465 # attribute_name = 'linear_combination') 466 # 467 # #Check contents 468 # #Get NetCDF 469 # from Scientific.IO.NetCDF import NetCDFFile 470 # fid = NetCDFFile(ptsfile, 'r') 471 # 472 # # Get the variables 473 # #print fid.variables.keys() 474 # points1 = fid.variables['points'] 475 # z1 = fid.variables['linear_combination'] 476 # 477 # #Check values# 478 # 479 # #print points[:] 480 # #print ref_points 481 # assert allclose(points, points1) 482 # 483 # #print attributes[:] 484 # #print ref_elevation 485 # assert allclose(z, z1) 486 # 487 # #Cleanup 488 # fid.close() 489 # 490 # import os 491 # os.remove(ptsfile) 490 492 491 493 -
inundation/pyvolution/test_least_squares.py
r2260 r2262 1408 1408 1409 1409 1410 1411 def test_fit_to_mesh_w_georef(self): 1412 """Simple check that georef works at the fit_to_mesh level 1413 """ 1414 1415 from coordinate_transforms.geo_reference import Geo_reference 1416 1417 #Mesh 1418 vertex_coordinates = [[0.76, 0.76], 1419 [0.76, 5.76], 1420 [5.76, 0.76]] 1421 triangles = [[0,2,1]] 1422 1423 mesh_geo = Geo_reference(56,-0.76,-0.76) 1424 1425 1426 #Data 1427 data_points = [[ 201.0, 401.0], 1428 [ 201.0, 403.0], 1429 [ 203.0, 401.0]] 1430 1431 z = [2, 4, 4] 1432 1433 data_geo = Geo_reference(56,-200,-400) 1434 1435 #Fit 1436 zz = fit_to_mesh(vertex_coordinates, triangles, data_points, z, 1437 data_origin = data_geo.get_origin(), 1438 mesh_origin = mesh_geo.get_origin(), 1439 alpha = 0) 1440 1441 assert allclose( zz, [0,5,5] ) 1442 1443 1410 1444 def test_fit_to_mesh_file(self): 1411 1445 from load_mesh.loadASCII import import_mesh_file, \ … … 1520 1554 mesh_dic = {} 1521 1555 mesh_dic['vertices'] = [[0.76, 0.76], 1522 1523 1556 [0.76, 5.76], 1557 [5.76, 0.76]] 1524 1558 mesh_dic['triangles'] = [[0, 2, 1]] 1525 1559 mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]] … … 1529 1563 mesh_dic['triangle_neighbors'] = [[-1, -1, -1]] 1530 1564 mesh_dic['segment_tags'] = ['external', 1531 1532 1565 'external', 1566 'external'] 1533 1567 mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76) 1534 1568 mesh_file = tempfile.mktemp(".tsh") … … 1553 1587 ans =[[0.0, 0.0], 1554 1588 [5.0, 10.0], 1555 [5.0, 10.0]]1589 [5.0, 10.0]] 1556 1590 assert allclose(mesh_dic['vertex_attributes'],ans) 1557 1591 -
inundation/pyvolution/test_quantity.py
r1916 r2262 300 300 z = linear_function(data_points) 301 301 302 #Obsoleted bit303 #interp = Interpolation(quantity.domain.coordinates,304 # quantity.domain.triangles,305 # data_points, alpha=0.0)306 #307 #answer = linear_function(quantity.domain.coordinates)308 #309 #f = interp.fit(z)310 #311 #print "f",f312 #print "answer",answer313 #assert allclose(f, answer)314 315 316 302 #Use built-in least squares fit 317 303 quantity.set_values(points = data_points, values = z, alpha = 0) … … 333 319 quantity.set_values(vertex_attributes) 334 320 assert allclose(quantity.vertex_values.flat, answer) 321 322 323 324 325 326 def test_test_set_values_using_least_squares_w_geo(self): 327 328 from domain import Domain 329 from coordinate_transforms.geo_reference import Geo_reference 330 from least_squares import fit_to_mesh 331 332 #Mesh 333 vertex_coordinates = [[0.76, 0.76], 334 [0.76, 5.76], 335 [5.76, 0.76]] 336 triangles = [[0,2,1]] 337 338 mesh_georef = Geo_reference(56,-0.76,-0.76) 339 mesh1 = Domain(vertex_coordinates, triangles, 340 geo_reference = mesh_georef) 341 mesh1.check_integrity() 342 343 #Quantity 344 quantity = Quantity(mesh1) 345 346 #Data 347 data_points = [[ 201.0, 401.0], 348 [ 201.0, 403.0], 349 [ 203.0, 401.0]] 350 351 z = [2, 4, 4] 352 353 data_georef = Geo_reference(56,-200,-400) 354 355 356 #Reference 357 ref = fit_to_mesh(vertex_coordinates, triangles, data_points, z, 358 data_origin = data_georef.get_origin(), 359 mesh_origin = mesh_georef.get_origin(), 360 alpha = 0) 361 362 assert allclose( ref, [0,5,5] ) 363 364 #Test set_values 365 quantity.set_values(points = data_points, values = z, data_georef = data_georef, 366 alpha = 0) 367 368 assert allclose(quantity.vertex_values.flat, ref) 369 335 370 336 371 … … 357 392 358 393 #Create pts file 359 from data_manager import write_ptsfile394 from load_mesh.loadASCII import export_points_file 360 395 ptsfile = 'testptsfile.pts' 361 396 att = 'spam_and_eggs' 362 write_ptsfile(ptsfile, data_points, z, attribute_name = att) 397 points_dict = {'pointlist': data_points, 398 'attributelist': {att: z}} 399 400 export_points_file(ptsfile, points_dict) 363 401 364 402 #Check that values can be set from file … … 366 404 attribute_name = att, alpha = 0) 367 405 answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True)) 368 #print quantity.vertex_values, answer 406 407 #print quantity.vertex_values.flat 408 #print answer 409 410 369 411 assert allclose(quantity.vertex_values.flat, answer) 370 412 … … 378 420 os.remove(ptsfile) 379 421 422 423 def test_set_values_from_file_with_georef1(self): 424 425 #Mesh in zone 56 (absolute coords) 426 from domain import Domain 427 from coordinate_transforms.geo_reference import Geo_reference 428 429 x0 = 314036.58727982 430 y0 = 6224951.2960092 431 432 a = [x0+0.0, y0+0.0] 433 b = [x0+0.0, y0+2.0] 434 c = [x0+2.0, y0+0.0] 435 d = [x0+0.0, y0+4.0] 436 e = [x0+2.0, y0+2.0] 437 f = [x0+4.0, y0+0.0] 438 439 points = [a, b, c, d, e, f] 440 441 #bac, bce, ecf, dbe 442 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 443 444 mesh4 = Domain(points, elements, 445 geo_reference = Geo_reference(56, 0, 0)) 446 mesh4.check_integrity() 447 quantity = Quantity(mesh4) 448 449 #Get (enough) datapoints (relative to georef) 450 data_points = [[ 0.66666667, 0.66666667], 451 [ 1.33333333, 1.33333333], 452 [ 2.66666667, 0.66666667], 453 [ 0.66666667, 2.66666667], 454 [ 0.0, 1.0], 455 [ 0.0, 3.0], 456 [ 1.0, 0.0], 457 [ 1.0, 1.0], 458 [ 1.0, 2.0], 459 [ 1.0, 3.0], 460 [ 2.0, 1.0], 461 [ 3.0, 0.0], 462 [ 3.0, 1.0]] 463 464 z = linear_function(data_points) 465 466 467 #Create pts file 468 from data_manager import write_ptsfile 469 from load_mesh.loadASCII import export_points_file 470 471 ptsfile = 'testptsfile.pts' 472 att = 'spam_and_eggs' 473 474 points_dict = {'pointlist': data_points, 475 'attributelist': {att: z}, 476 'geo_reference': Geo_reference(zone = 56, 477 xllcorner = x0, 478 yllcorner = y0)} 479 480 export_points_file(ptsfile, points_dict) 481 482 483 #Check that values can be set from file 484 quantity.set_values(filename = ptsfile, 485 attribute_name = att, alpha = 0) 486 answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True) - [x0, y0]) 487 488 assert allclose(quantity.vertex_values.flat, answer) 489 490 491 #Check that values can be set from file using default attribute 492 quantity.set_values(filename = ptsfile, alpha = 0) 493 assert allclose(quantity.vertex_values.flat, answer) 494 495 #Cleanup 496 import os 497 os.remove(ptsfile) 498 499 500 def test_set_values_from_file_with_georef2(self): 501 502 #Mesh in zone 56 (relative coords) 503 from domain import Domain 504 from coordinate_transforms.geo_reference import Geo_reference 505 506 x0 = 314036.58727982 507 y0 = 6224951.2960092 508 509 a = [0.0, 0.0] 510 b = [0.0, 2.0] 511 c = [2.0, 0.0] 512 d = [0.0, 4.0] 513 e = [2.0, 2.0] 514 f = [4.0, 0.0] 515 516 points = [a, b, c, d, e, f] 517 518 #bac, bce, ecf, dbe 519 elements = [ [1,0,2], [1,2,4], [4,2,5], [3,1,4] ] 520 521 mesh4 = Domain(points, elements, 522 geo_reference = Geo_reference(56, x0, y0)) 523 mesh4.check_integrity() 524 quantity = Quantity(mesh4) 525 526 #Get (enough) datapoints (relative to georef) 527 data_points = [[ x0+0.66666667, y0+0.66666667], 528 [ x0+1.33333333, y0+1.33333333], 529 [ x0+2.66666667, y0+0.66666667], 530 [ x0+0.66666667, y0+2.66666667], 531 [ x0+0.0, y0+1.0], 532 [ x0+0.0, y0+3.0], 533 [ x0+1.0, y0+0.0], 534 [ x0+1.0, y0+1.0], 535 [ x0+1.0, y0+2.0], 536 [ x0+1.0, y0+3.0], 537 [ x0+2.0, y0+1.0], 538 [ x0+3.0, y0+0.0], 539 [ x0+3.0, y0+1.0]] 540 541 z = linear_function(data_points) 542 543 544 #Create pts file 545 from data_manager import write_ptsfile 546 from load_mesh.loadASCII import export_points_file 547 548 ptsfile = 'testptsfile.pts' 549 att = 'spam_and_eggs' 550 551 points_dict = {'pointlist': data_points, 552 'attributelist': {att: z}, 553 'geo_reference': Geo_reference(zone = 56, 554 xllcorner = 0, 555 yllcorner = 0)} 556 557 export_points_file(ptsfile, points_dict) 558 559 560 #Check that values can be set from file 561 quantity.set_values(filename = ptsfile, 562 attribute_name = att, alpha = 0) 563 answer = linear_function(quantity.domain.get_vertex_coordinates(obj = True) + [x0, y0]) 564 565 566 assert allclose(quantity.vertex_values.flat, answer) 567 568 569 #Check that values can be set from file using default attribute 570 quantity.set_values(filename = ptsfile, alpha = 0) 571 assert allclose(quantity.vertex_values.flat, answer) 572 573 #Cleanup 574 import os 575 os.remove(ptsfile) 576 577 380 578 381 579 … … 1230 1428 #------------------------------------------------------------- 1231 1429 if __name__ == "__main__": 1232 suite = unittest.makeSuite(Test_Quantity, 'test')1430 suite = unittest.makeSuite(Test_Quantity, 'test')#_test_set_values_using_least_squares_w_geo') 1233 1431 #print "restricted test" 1234 1432 #suite = unittest.makeSuite(Test_Quantity,'test_set_values_func')
Note: See TracChangeset
for help on using the changeset viewer.
