Changeset 2690


Ignore:
Timestamp:
Apr 11, 2006, 1:13:32 PM (18 years ago)
Author:
duncan
Message:

added support for geospatial and geo_ref

Location:
inundation/fit_interpolate
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • inundation/fit_interpolate/interpolate.py

    r2684 r2690  
    3939                 vertex_coordinates,
    4040                 triangles,
    41                  mesh_origin = None,
     41                 mesh_origin=None,
    4242                 verbose=False,
    4343                 max_vertices_per_cell=30):
     
    6464          max_vertices_per_cell: Number of vertices in a quad tree cell
    6565          at which the cell is split into 4.
    66 
    67         """
     66        """
     67       
    6868        # Initialise variabels
    6969        self._A_can_be_reused = False
     
    7272        #Convert input to Numeric arrays
    7373        triangles = ensure_numeric(triangles, Int)
    74         vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
     74       
     75        if isinstance(vertex_coordinates,Geospatial_data):
     76            vertex_coordinates = vertex_coordinates.get_data_points( \
     77                absolute = True)
     78            msg = "Use a Geospatial_data object or a meshorigin. Not both."
     79            assert mesh_origin == None, msg
     80           
     81        else:
     82            vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
    7583        #Build underlying mesh
    7684        if verbose: print 'Building mesh'
     
    7886        #FIXME: Trying the normal mesh while testing precrop,
    7987        #       The functionality of boundary_polygon is needed for that
    80 
    81         #FIXME: geo_ref can also be a geo_ref object
    82         #FIXME: move this to interpolate_block
    83 
    84         #Note: not so keen to use geospacial to represent the mesh,
    85         # since the data structure is
    86         # points, geo-ref and triangles, rather than just points and geo-ref
    87         # this info will usually be coming from a domain instance.
    8888       
    8989        if mesh_origin is None:
    90             geo = None
     90            geo = None #Geo_reference()
    9191        else:
    92             geo = Geo_reference(mesh_origin[0],mesh_origin[1],mesh_origin[2])
    93         self.mesh = Mesh(vertex_coordinates, triangles,
    94                          geo_reference = geo)
    95        
     92            if isinstance(mesh_origin, Geo_reference):
     93                geo = mesh_origin
     94            else:
     95                geo = Geo_reference(mesh_origin[0],
     96                                    mesh_origin[1],
     97                                    mesh_origin[2])
     98            vertex_coordinates = geo.get_absolute(vertex_coordinates)
     99        #Don't pass geo_reference to mesh.  It doesn't work.
     100        self.mesh = Mesh(vertex_coordinates, triangles)
    96101        self.mesh.check_integrity()
    97102        self.root = build_quadtree(self.mesh,
     
    128133        """
    129134        #print "point_coordinates interpolate.interpolate",point_coordinates
    130         # Can I interpolate, based on previous point_coordinates?
    131135        if isinstance(point_coordinates,Geospatial_data):
    132136            point_coordinates = point_coordinates.get_data_points( \
    133137                absolute = True)
    134138       
     139        # Can I interpolate, based on previous point_coordinates?
    135140        if point_coordinates is None:
    136141            if self._A_can_be_reused is True and \
  • inundation/fit_interpolate/test_interpolate.py

    r2684 r2690  
    191191                    [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
    192192                    [0.0, 0.5, 0.5]]
    193         interp = Interpolate(points, vertices, data)
     193        interp = Interpolate(points, vertices)
    194194
    195195        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
     
    212212                   [0.0, 0.25, 0.75]]
    213213
    214         interp = Interpolate(points, vertices, data)
     214        interp = Interpolate(points, vertices)
    215215
    216216        assert allclose(interp._build_interpolation_matrix_A(data).todense(),
     
    232232        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
    233233
    234         interp = Interpolate(points, vertices, data)
     234        interp = Interpolate(points, vertices)
    235235        #print "interp.get_A()", interp.get_A()
    236236        results = interp._build_interpolation_matrix_A(data).todense()
     
    252252        data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
    253253
    254         interp = Interpolate(points, vertices, data)
     254        interp = Interpolate(points, vertices)
    255255        results = interp._build_interpolation_matrix_A(data).todense()
    256256        assert allclose(sum(results, axis=1), [1,1,1,0])
     
    293293
    294294        assert allclose(A, answer)
    295 
    296         #FIXME -  tests are hopefully failing due to points outsdie the
    297         # mesh
     295   
     296    def test_geo_ref(self):
     297        v0 = [0.0, 0.0]
     298        v1 = [0.0, 5.0]
     299        v2 = [5.0, 0.0]
     300
     301        vertices_absolute = [v0, v1, v2]
     302        triangles = [ [1,0,2] ]   #bac
     303
     304        geo = Geo_reference(57,100, 500)
     305
     306        vertices = geo.change_points_geo_ref(vertices_absolute)
     307        #print "vertices",vertices
     308       
     309        d0 = [1.0, 1.0]
     310        d1 = [1.0, 2.0]
     311        d2 = [3.0, 1.0]
     312        point_coords = [ d0, d1, d2]
     313
     314        interp = Interpolate(vertices, triangles, mesh_origin=geo)
     315        f = linear_function(vertices_absolute)
     316        z = interp.interpolate(f, point_coords)
     317        answer = linear_function(point_coords)
     318
     319        #print "z",z
     320        #print "answer",answer
     321        assert allclose(z, answer)
     322       
     323    def test_Geospatial_verts(self):
     324        v0 = [0.0, 0.0]
     325        v1 = [0.0, 5.0]
     326        v2 = [5.0, 0.0]
     327
     328        vertices_absolute = [v0, v1, v2]
     329        triangles = [ [1,0,2] ]   #bac
     330
     331        geo = Geo_reference(57,100, 500)
     332        vertices = geo.change_points_geo_ref(vertices_absolute)
     333        geopoints = Geospatial_data(vertices,geo_reference = geo)
     334        #print "vertices",vertices
     335       
     336        d0 = [1.0, 1.0]
     337        d1 = [1.0, 2.0]
     338        d2 = [3.0, 1.0]
     339        point_coords = [ d0, d1, d2]
     340
     341        interp = Interpolate(geopoints, triangles)
     342        f = linear_function(vertices_absolute)
     343        z = interp.interpolate(f, point_coords)
     344        answer = linear_function(point_coords)
     345
     346        #print "z",z
     347        #print "answer",answer
     348        assert allclose(z, answer)
     349       
    298350    def test_interpolate_attributes_to_points(self):
    299351        v0 = [0.0, 0.0]
     
    309361        point_coords = [ d0, d1, d2]
    310362
    311         interp = Interpolate(vertices, triangles, point_coords)
     363        interp = Interpolate(vertices, triangles)
    312364        f = linear_function(vertices)
    313365        z = interp.interpolate(f, point_coords)
    314366        answer = linear_function(point_coords)
    315367
    316 
     368        #print "z",z
     369        #print "answer",answer
    317370        assert allclose(z, answer)
    318371
Note: See TracChangeset for help on using the changeset viewer.