Changeset 2262


Ignore:
Timestamp:
Jan 20, 2006, 11:04:18 AM (18 years ago)
Author:
ole
Message:

Fixed missing geo reference in set_quantity and added tests

Location:
inundation
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • inundation/coordinate_transforms/geo_reference.py

    r2261 r2262  
    2121                 yllcorner = 0.0,
    2222                 datum = 'wgs84',
    23                  projection = 'UTM',
    24                  units = 'm',
     23                 projection = 'UTM',                 units = 'm',
    2524                 false_easting = 500000,
    2625                 false_northing = 10000000, #Default for southern hemisphere
  • inundation/load_mesh/loadASCII.py

    r2257 r2262  
    11701170    # There can't be an attribute called points
    11711171    # 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
    11721178   
    11731179    from Scientific.IO.NetCDF import NetCDFFile
  • inundation/pyvolution/data_manager.py

    r2256 r2262  
    10231023    """Write points and associated attribute to pts (NetCDF) format
    10241024    """
     1025
     1026    print 'WARNING: write_ptsfile is obsolete. Use export_points from load_mesh.loadASCII instead.'
    10251027
    10261028    from Numeric import Float
  • inundation/pyvolution/quantity.py

    r2233 r2262  
    167167                   quantity = None,   # Another quantity
    168168                   function = None,   # Callable object: f(x,y)
    169                    points = None, values = None, #Input for least squares
     169                   points = None, values = None, data_georef = None, #Input for least squares
    170170                   filename = None, attribute_name = None, #Input from file
    171171                   alpha = None,
     
    201201          If points is specified, values is an array of length N containing
    202202          attribute values for each point.
     203
     204        data_georef:
     205          If points is specified, geo_reference applies to each point.       
    203206         
    204207        filename:
     
    305308            assert values is not None, msg
    306309            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)
    309314        elif filename is not None:
    310315            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)
    313319        else:
    314320            raise 'This can\'t happen :-)'
     
    513519        #FIXME: Should check that function returns something sensible and
    514520        #raise a meaningfull exception if it returns None for example
     521
     522        #FIXME: Should supply absolute coordinates
    515523
    516524        from Numeric import take
     
    546554
    547555    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):
    549560        """Set quantity values from arbitray data points using least squares
    550561        """
     
    553564        from util import ensure_numeric
    554565        from least_squares import fit_to_mesh
     566        from coordinate_transforms.geo_reference import Geo_reference
     567       
    555568
    556569        points = ensure_numeric(points, Float)
     
    565578        triangles = self.domain.triangles
    566579
     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       
    567602        if use_cache is True:
    568603            try:
     
    573608                raise msg
    574609
    575             args = (coordinates, triangles, points, values)
    576             kwargs = {'alpha': alpha, 'verbose': verbose}
    577610            vertex_attributes = cache(fit_to_mesh,
    578611                                      args, kwargs,
    579612                                      verbose = verbose)
    580613        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   
    589619        self.set_values_from_array(vertex_attributes,
    590620                                   location, indices, verbose)
     
    595625
    596626    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):
    598630        """Set quantity based on arbitrary points in .pts file
    599631        using least_squares attribute_name selects name of attribute
     
    602634        """
    603635
     636        from load_mesh.loadASCII import import_points_file
    604637
    605638        from types import StringType
     
    609642
    610643        #Read from (NetCDF) file
    611         from load_mesh.loadASCII import import_points_file
    612644        points_dict = import_points_file(filename)
    613645        points = points_dict['pointlist']
     
    634666
    635667
    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
    637677        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)
    639682
    640683
  • inundation/pyvolution/test_data_manager.py

    r2074 r2262  
    439439
    440440
    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)
    490492       
    491493
  • inundation/pyvolution/test_least_squares.py

    r2260 r2262  
    14081408
    14091409
     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
    14101444    def test_fit_to_mesh_file(self):
    14111445        from load_mesh.loadASCII import import_mesh_file, \
     
    15201554        mesh_dic = {}
    15211555        mesh_dic['vertices'] = [[0.76, 0.76],
    1522                                           [0.76, 5.76],
    1523                                           [5.76, 0.76]]
     1556                                [0.76, 5.76],
     1557                                [5.76, 0.76]]
    15241558        mesh_dic['triangles'] =  [[0, 2, 1]]
    15251559        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
     
    15291563        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
    15301564        mesh_dic['segment_tags'] = ['external',
    1531                                                   'external',
    1532                                                   'external']
     1565                                    'external',
     1566                                    'external']
    15331567        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
    15341568        mesh_file = tempfile.mktemp(".tsh")
     
    15531587        ans =[[0.0, 0.0],
    15541588              [5.0, 10.0],
    1555               [5.0,10.0]]
     1589              [5.0, 10.0]]
    15561590        assert allclose(mesh_dic['vertex_attributes'],ans)
    15571591
  • inundation/pyvolution/test_quantity.py

    r1916 r2262  
    300300        z = linear_function(data_points)
    301301       
    302         #Obsoleted bit
    303         #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",f
    312         #print "answer",answer
    313         #assert allclose(f, answer)
    314 
    315 
    316302        #Use built-in least squares fit
    317303        quantity.set_values(points = data_points, values = z, alpha = 0)
     
    333319        quantity.set_values(vertex_attributes)
    334320        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
    335370
    336371
     
    357392
    358393        #Create pts file
    359         from data_manager import write_ptsfile
     394        from load_mesh.loadASCII import export_points_file
    360395        ptsfile = 'testptsfile.pts'
    361396        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)
    363401
    364402        #Check that values can be set from file
     
    366404                            attribute_name = att, alpha = 0)
    367405        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       
    369411        assert allclose(quantity.vertex_values.flat, answer)
    370412
     
    378420        os.remove(ptsfile)
    379421       
     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
    380578
    381579
     
    12301428#-------------------------------------------------------------
    12311429if __name__ == "__main__":
    1232     suite = unittest.makeSuite(Test_Quantity,'test')
     1430    suite = unittest.makeSuite(Test_Quantity, 'test')#_test_set_values_using_least_squares_w_geo')
    12331431    #print "restricted test"
    12341432    #suite = unittest.makeSuite(Test_Quantity,'test_set_values_func')
Note: See TracChangeset for help on using the changeset viewer.