Changeset 1178


Ignore:
Timestamp:
Mar 31, 2005, 4:51:54 PM (20 years ago)
Author:
duncan
Message:

incorporating geo_ref object into domain class

Location:
inundation/ga/storm_surge/pyvolution
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • inundation/ga/storm_surge/pyvolution/data_manager.py

    r1171 r1178  
    5050from Numeric import concatenate
    5151
     52from coordinate_transforms.geo_reference import Geo_reference, DEFAULT_ZONE
    5253
    5354def make_filename(s):
    … …  
    248249            #FIXME: Use Georef
    249250            fid.starttime = domain.starttime
    250             fid.xllcorner = domain.xllcorner
    251             fid.yllcorner = domain.yllcorner
    252             fid.zone = domain.zone
    253251
    254252            # dimension definitions
    … …  
    267265            fid.createVariable('y', self.precision, ('number_of_points',))
    268266            fid.createVariable('elevation', self.precision, ('number_of_points',))
    269 
     267            if domain.geo_reference is not None:
     268                domain.geo_reference.write_NetCDF(fid)
    270269            #FIXME: Backwards compatibility
    271270            fid.createVariable('z', self.precision, ('number_of_points',))
    … …  
    10771076    false_northing = 10000000
    10781077    NODATA_value = -9999
    1079     #
    1080 
    10811078   
    10821079    if quantity is None:
    … …  
    11101107    number_of_points = fid.dimensions['number_of_points']
    11111108    if origin is None:
    1112         xllcorner = fid.xllcorner[0]
    1113         yllcorner = fid.yllcorner[0]
    1114         zone = fid.zone[0]
     1109       
     1110        # get geo_reference
     1111        #sww files don't have to have a geo_ref
     1112        try:
     1113            geo_reference = Geo_reference(NetCDFObject=fid)
     1114        except AttributeError, e:
     1115            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
     1116       
     1117        xllcorner = geo_reference.get_xllcorner()
     1118        yllcorner = geo_reference.get_yllcorner()
     1119        zone = geo_reference.get_zone()
    11151120    else:
    11161121        zone = origin[0]
    … …  
    18721877    ymomentum = fid.variables['ymomentum']   #Momentum in the y-direction
    18731878
    1874     xllcorner = fid.xllcorner[0]
    1875     yllcorner = fid.yllcorner[0]
    18761879    starttime = fid.starttime[0]
    1877     zone = fid.zone
    18781880    volumes = fid.variables['volumes'][:] #Connectivity
    18791881    coordinates=transpose(asarray([x.tolist(),y.tolist()]))
    … …  
    18831885    other_quantities = []
    18841886
     1887    # get geo_reference
     1888    #sww files don't have to have a geo_ref
     1889    try:
     1890        geo_reference = Geo_reference(NetCDFObject=fid)
     1891    except AttributeError, e:
     1892        geo_reference = None
     1893       
    18851894    if verbose: print '    interpolating quantities'
    18861895    for quantity in fid.variables.keys():
    … …  
    19161925    if not boundary is None:
    19171926        domain.boundary = boundary
    1918     domain.zone=zone
    1919     domain.xllcorner=float(xllcorner)
    1920     domain.yllcorner=float(yllcorner)
     1927   
     1928    domain.geo_reference = geo_reference
    19211929
    19221930    domain.starttime=float(starttime)+float(t)

  • inundation/ga/storm_surge/pyvolution/domain.py

    r1140 r1178  
    1616    def __init__(self, coordinates, vertices, boundary = None,
    1717                 conserved_quantities = None, other_quantities = None,
    18                  tagged_elements = None, zone=0, xllcorner=0, yllcorner=0):
    19 
    20         Mesh.__init__(self, coordinates, vertices, boundary, tagged_elements)
     18                 tagged_elements = None, geo_reference = None):
     19
     20        Mesh.__init__(self, coordinates, vertices, boundary,
     21                      tagged_elements, geo_reference)
    2122
    2223        from Numeric import zeros, Float
    … …  
    7374        #Origin in UTM coordinates
    7475        #FIXME: This should be set if read by a msh file
    75         self.zone = zone
    76         self.xllcorner = xllcorner
    77         self.yllcorner = yllcorner
     76        #self.zone = zone
     77        #self.xllcorner = xllcorner
     78        #self.yllcorner = yllcorner
    7879
    7980

  • inundation/ga/storm_surge/pyvolution/general_mesh.py

    r1011 r1178  
    1    
     1
     2from coordinate_transforms.geo_reference import Geo_reference,DEFAULT_ZONE
     3
    24class General_mesh:
    35    """Collection of triangular elements (purely geometric)
    … …  
    4850    """
    4951
    50     def __init__(self, coordinates, triangles, origin = None):
     52    def __init__(self, coordinates, triangles,
     53                 geo_reference=None):
    5154        """
    5255        Build triangles from x,y coordinates (sequence of 2-tuples or
    … …  
    6366        self.coordinates = array(coordinates)
    6467
    65         self.origin = origin
     68        if geo_reference is None:
     69            self.geo_reference = Geo_reference(DEFAULT_ZONE,0.0,0.0)
     70        else:
     71            self.geo_reference = geo_reference
    6672
    6773        #Input checks

  • inundation/ga/storm_surge/pyvolution/least_squares.py

    r1160 r1178  
    2828from sparse import Sparse, Sparse_CSR
    2929from cg_solve import conjugate_gradient, VectorShapeError
     30
     31from coordinate_transforms.geo_reference import Geo_reference
     32
    3033import time
     34
    3135
    3236try:
    … …  
    303307        #FIXME: Trying the normal mesh while testing precrop,
    304308        #       The functionality of boundary_polygon is needed for that
     309
     310        #FIXME - geo ref does not have to go into mesh.
     311        # Change the point co-ords to conform to the
     312        # mesh co-ords early in the code
     313        if mesh_origin == None:
     314            geo = None
     315        else:
     316            geo = Geo_reference(mesh_origin[0],mesh_origin[1],mesh_origin[2])
    305317        self.mesh = Mesh(vertex_coordinates, triangles,
    306                          origin = mesh_origin)
     318                         geo_reference = geo)
    307319        #FIXME, remove if mesh checks it.
    308320        self.mesh.check_integrity()
    … …  
    392404
    393405        #Shift data points to same origin as mesh (if specified)
    394         mesh_origin = self.mesh.origin
     406
     407        #FIXME this will shift if there was no geo_ref.
     408        #But all this should be removed amyhow.
     409        #change coords before this point
     410        mesh_origin = self.mesh.geo_reference.get_origin()
    395411        if point_coordinates is not None:
    396412            if data_origin is not None:

  • inundation/ga/storm_surge/pyvolution/mesh.py

    r1160 r1178  
    5858
    5959    def __init__(self, coordinates, triangles, boundary = None,
    60                  tagged_elements = None, origin = None):
     60                 tagged_elements = None, geo_reference = None):
    6161        """
    6262        Build triangles from x,y coordinates (sequence of 2-tuples or
    … …  
    6767        from Numeric import array, zeros, Int, Float, maximum, sqrt, sum
    6868
    69         General_mesh.__init__(self, coordinates, triangles, origin)
     69        General_mesh.__init__(self, coordinates, triangles, geo_reference)
    7070
    7171        N = self.number_of_elements

  • inundation/ga/storm_surge/pyvolution/pmesh2domain.py

    r969 r1178  
    1414    from domain import Domain
    1515
    16     vertex_coordinates, volumes, tag_dict, vertex_quantity_dict ,tagged_elements_dict = \
    17                         pmesh_to_domain(fileName,
    18                                                    setting_function=setting_function)
     16    vertex_coordinates, volumes, tag_dict, vertex_quantity_dict \
     17                        ,tagged_elements_dict, geo_reference = \
     18           pmesh_to_domain(fileName,
     19                           setting_function=setting_function)
    1920       
    20     assert issubclass(DomainClass, Domain), "DomainClass is not a subclass of Domain."
     21    assert issubclass(DomainClass, Domain),"DomainClass is not a subclass of Domain."
    2122
    2223
    2324   
    2425    domain = DomainClass(vertex_coordinates, volumes, tag_dict,
    25                          tagged_elements = tagged_elements_dict )
     26                         tagged_elements = tagged_elements_dict,
     27                         geo_reference = geo_reference )
    2628
    2729
    … …  
    3537    # This doesn't work on the domain instance.
    3638    # This is still needed so -ve elevations don't cuase 'lakes'
    37     # The fixme we discussed was to only create a quantity when its values are set.
     39    # The fixme we discussed was to only create a quantity when its values
     40    #are set.
    3841    # I think that's the way to go still
    3942   
    … …  
    7073    volumes = mesh_dict['triangles']
    7174
    72     #if setting_function:
    73     #    if not type(setting_function) is ListType:
    74     #        setting_function = [setting_function]
    75     #    for funct in setting_function:
    76     #        mesh_dict = funct(mesh_dict, vertices = mesh_vertices,
    77     #                        volumes = volumes)
    78 
    79 
    8075    vertex_quantity_dict = {}
    8176    point_atts = transpose(mesh_dict['vertex_attributes'])
    8277    point_titles  = mesh_dict['vertex_attribute_titles']
     78    geo_reference  = mesh_dict['geo_reference']
    8379    for quantity, value_vector in map (None, point_titles, point_atts):
    8480        vertex_quantity_dict[quantity] = value_vector
    8581    tag_dict = pmesh_dict_to_tag_dict(mesh_dict)
    8682    tagged_elements_dict = build_tagged_elements_dictionary(mesh_dict)
    87     return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict,tagged_elements_dict
     83    return vertex_coordinates, volumes, tag_dict, vertex_quantity_dict,tagged_elements_dict, geo_reference
    8884
    8985

  • inundation/ga/storm_surge/pyvolution/shallow_water.py

    r1087 r1178  
    6363
    6464    def __init__(self, coordinates, vertices, boundary = None,
    65                  tagged_elements = None):
     65                 tagged_elements = None, geo_reference = None):
    6666
    6767        conserved_quantities = ['stage', 'xmomentum', 'ymomentum']
    6868        other_quantities = ['elevation', 'friction']
    69 
    7069        Generic_domain.__init__(self, coordinates, vertices, boundary,
    7170                                conserved_quantities, other_quantities,
    72                                 tagged_elements)
     71                                tagged_elements, geo_reference)
    7372
    7473        from config import minimum_allowed_height, g

  • inundation/ga/storm_surge/pyvolution/test_data_manager.py

    r1171 r1178  
    1010from shallow_water import *
    1111from config import epsilon
     12
     13from coordinate_transforms.geo_reference import Geo_reference
    1214
    1315class Test_Data_Manager(unittest.TestCase):
    … …  
    604606        self.domain.set_quantity('elevation', lambda x,y: -x-y)
    605607
    606         self.domain.xllcorner = 308500
    607         self.domain.yllcorner = 6189000
    608         self.domain.zone = 56
    609        
     608        self.domain.geo_reference = Geo_reference(56,308500,6189000)
    610609       
    611610        sww = get_dataobject(self.domain)
    … …  
    746745        self.domain.set_quantity('elevation', lambda x,y: -x-y)
    747746
    748         self.domain.xllcorner = 308500
    749         self.domain.yllcorner = 6189000
    750         self.domain.zone = 56
     747        self.domain.geo_reference = Geo_reference(56,308500,6189000)
    751748       
    752749       
    … …  
    897894        domain.smooth = True
    898895
    899 
    900         domain.xllcorner = 308500
    901         domain.yllcorner = 6189000
    902         domain.zone = 56
    903        
     896        domain.geo_reference = Geo_reference(56,308500,6189000)
    904897       
    905898        sww = get_dataobject(domain)
    … …  
    919912        time = fid.variables['time'][:]
    920913
     914        try:
     915            geo_reference = Geo_reference(NetCDFObject=fid)
     916        except AttributeError, e:
     917            geo_reference = Geo_reference(DEFAULT_ZONE,0,0)
     918       
    921919        #Export to ascii/prj files
    922920        sww2asc(domain.filename,
    … …  
    13451343        yiel=0.01
    13461344        points, vertices, boundary = rectangular(10,10)
     1345       
    13471346        #Create shallow water domain
    13481347        domain = Domain(points, vertices, boundary)
     1348        domain.geo_reference = Geo_reference(56,11,11)
    13491349        domain.smooth = False
    13501350        domain.visualise = False
    … …  
    13951395        ###################
    13961396
    1397         bits = ['xllcorner','yllcorner','vertex_coordinates']
    1398 
     1397        bits = [ 'geo_reference.get_xllcorner()',
     1398                'geo_reference.get_yllcorner()',
     1399                'vertex_coordinates']
    13991400        for quantity in ['elevation']+domain.quantities_to_be_stored:
    14001401            bits.append('quantities["%s"].get_integral()'%quantity)
    … …  
    14021403
    14031404        for bit in bits:
    1404         #    print 'testing that domain.'+bit+' has been restored'
    1405         #    print bit
     1405            #print 'testing that domain.'+bit+' has been restored'
     1406            #print bit
    14061407        #print 'done'
    14071408            assert allclose(eval('domain.'+bit),eval('domain2.'+bit))
    … …  
    14521453        ##################
    14531454
    1454         bits = ['xllcorner','yllcorner','vertex_coordinates']
     1455        bits = [ 'geo_reference.get_xllcorner()',
     1456                'geo_reference.get_yllcorner()',
     1457                'vertex_coordinates']
    14551458
    14561459        for quantity in ['elevation','xmomentum','ymomentum']:#+domain.quantities_to_be_stored:
    … …  
    15281531            filler = 0
    15291532            domain2 = sww2domain(filename,None,fail_if_NaN=False,NaN_filler = filler,verbose=False)
    1530 
    1531         bits = ['xllcorner','yllcorner','vertex_coordinates']
     1533        bits = [ 'geo_reference.get_xllcorner()',
     1534                'geo_reference.get_yllcorner()',
     1535                'vertex_coordinates']
    15321536
    15331537        for quantity in ['elevation']+domain.quantities_to_be_stored:
    … …  
    15691573if __name__ == "__main__":
    15701574    suite = unittest.makeSuite(Test_Data_Manager,'test')   
    1571     #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2asc_stage')
     1575    #suite = unittest.makeSuite(Test_Data_Manager,'test_sww2domain')
    15721576    runner = unittest.TextTestRunner()
    15731577    runner.run(suite)

  • inundation/ga/storm_surge/pyvolution/test_pmesh2domain.py

    r1028 r1178  
    2222from shallow_water import Domain, Reflective_boundary, Dirichlet_boundary,\
    2323     Transmissive_boundary
     24
     25from coordinate_transforms.geo_reference import Geo_reference
    2426
    2527
    … …  
    64662 2 0 0 \n \
    65670 # <x> <y> ...Mesh Holes... \n \
    66 0 # <x> <y> <attribute>...Mesh Regions... \n")
     680 # <x> <y> <attribute>...Mesh Regions... \n \
     690 # <x> <y> <attribute>...Mesh Regions, area... \n \
     70# Geo reference \n \
     7156 \n \
     72140 \n \
     73120 \n")
    6774         file.close()
    6875
    … …  
    112119         self.failUnless( len(domain.boundary)  == 4,
    113120                          "test_pmesh2Domain Too many boundaries")
     121         #FIXME change to use get_xllcorner
     122         self.failUnless(domain.geo_reference.xllcorner  == 140.0,
     123                          "bad geo_referece")
     124         
    114125    def old_test_tags_to_boundaries (self):
    115126         meshDict = {}
Note: See TracChangeset for help on using the changeset viewer.