Changeset 2690

Timestamp:

Apr 11, 2006, 1:13:32 PM (19 years ago)

Author:

duncan

Message:

added support for geospatial and geo_ref

Location:

inundation/fit_interpolate

Files:

• interpolate.py (modified) (5 diffs)
• test_interpolate.py (modified) (6 diffs)

inundation/fit_interpolate/interpolate.py

@@ -39 +39 @@
                  vertex_coordinates,
                  triangles,
-                 mesh_origin = None,
+                 mesh_origin=None,
                  verbose=False,
                  max_vertices_per_cell=30):
@@ -64 +64 @@
           max_vertices_per_cell: Number of vertices in a quad tree cell
           at which the cell is split into 4.
-
-        """
+        """
+        
         # Initialise variabels
         self._A_can_be_reused = False
@@ -72 +72 @@
         #Convert input to Numeric arrays
         triangles = ensure_numeric(triangles, Int)
-        vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
+        
+        if isinstance(vertex_coordinates,Geospatial_data):
+            vertex_coordinates = vertex_coordinates.get_data_points( \
+                absolute = True)
+            msg = "Use a Geospatial_data object or a meshorigin. Not both."
+            assert mesh_origin == None, msg
+            
+        else:
+            vertex_coordinates = ensure_numeric(vertex_coordinates, Float)
         #Build underlying mesh
         if verbose: print 'Building mesh'
@@ -78 +86 @@
         #FIXME: Trying the normal mesh while testing precrop,
         #       The functionality of boundary_polygon is needed for that
-
-        #FIXME: geo_ref can also be a geo_ref object
-        #FIXME: move this to interpolate_block
-
-        #Note: not so keen to use geospacial to represent the mesh,
-        # since the data structure is
-        # points, geo-ref and triangles, rather than just points and geo-ref
-        # this info will usually be coming from a domain instance.
         
         if mesh_origin is None:
-            geo = None
+            geo = None #Geo_reference()
         else:
-            geo = Geo_reference(mesh_origin[0],mesh_origin[1],mesh_origin[2])
-        self.mesh = Mesh(vertex_coordinates, triangles,
-                         geo_reference = geo)
-        
+            if isinstance(mesh_origin, Geo_reference):
+                geo = mesh_origin
+            else:
+                geo = Geo_reference(mesh_origin[0],
+                                    mesh_origin[1],
+                                    mesh_origin[2])
+            vertex_coordinates = geo.get_absolute(vertex_coordinates)
+        #Don't pass geo_reference to mesh.  It doesn't work.
+        self.mesh = Mesh(vertex_coordinates, triangles)
         self.mesh.check_integrity()
         self.root = build_quadtree(self.mesh,
@@ -128 +133 @@
         """
         #print "point_coordinates interpolate.interpolate",point_coordinates 
-        # Can I interpolate, based on previous point_coordinates?
         if isinstance(point_coordinates,Geospatial_data):
             point_coordinates = point_coordinates.get_data_points( \
                 absolute = True)
         
+        # Can I interpolate, based on previous point_coordinates?
         if point_coordinates is None:
             if self._A_can_be_reused is True and \

inundation/fit_interpolate/test_interpolate.py

@@ -191 +191 @@
                     [0.5, 0.0, 0.5],  #Affects vertex 0 and 2
                     [0.0, 0.5, 0.5]]
-        interp = Interpolate(points, vertices, data)
+        interp = Interpolate(points, vertices)
 
         assert allclose(interp._build_interpolation_matrix_A(data).todense(),
@@ -212 +212 @@
                    [0.0, 0.25, 0.75]]
 
-        interp = Interpolate(points, vertices, data)
+        interp = Interpolate(points, vertices)
 
         assert allclose(interp._build_interpolation_matrix_A(data).todense(),
@@ -232 +232 @@
         data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44] ]
 
-        interp = Interpolate(points, vertices, data)
+        interp = Interpolate(points, vertices)
         #print "interp.get_A()", interp.get_A()
         results = interp._build_interpolation_matrix_A(data).todense()
@@ -252 +252 @@
         data = [ [0.2, 1.5], [0.123, 1.768], [1.43, 0.44], [5.0, 7.0]]
 
-        interp = Interpolate(points, vertices, data)
+        interp = Interpolate(points, vertices)
         results = interp._build_interpolation_matrix_A(data).todense()
         assert allclose(sum(results, axis=1), [1,1,1,0])
@@ -293 +293 @@
 
         assert allclose(A, answer)
-
-        #FIXME -  tests are hopefully failing due to points outsdie the
-        # mesh
+    
+    def test_geo_ref(self):
+        v0 = [0.0, 0.0]
+        v1 = [0.0, 5.0]
+        v2 = [5.0, 0.0]
+
+        vertices_absolute = [v0, v1, v2]
+        triangles = [ [1,0,2] ]   #bac
+
+        geo = Geo_reference(57,100, 500)
+
+        vertices = geo.change_points_geo_ref(vertices_absolute)
+        #print "vertices",vertices 
+        
+        d0 = [1.0, 1.0]
+        d1 = [1.0, 2.0]
+        d2 = [3.0, 1.0]
+        point_coords = [ d0, d1, d2]
+
+        interp = Interpolate(vertices, triangles, mesh_origin=geo)
+        f = linear_function(vertices_absolute)
+        z = interp.interpolate(f, point_coords)
+        answer = linear_function(point_coords)
+
+        #print "z",z 
+        #print "answer",answer 
+        assert allclose(z, answer)
+        
+    def test_Geospatial_verts(self):
+        v0 = [0.0, 0.0]
+        v1 = [0.0, 5.0]
+        v2 = [5.0, 0.0]
+
+        vertices_absolute = [v0, v1, v2]
+        triangles = [ [1,0,2] ]   #bac
+
+        geo = Geo_reference(57,100, 500)
+        vertices = geo.change_points_geo_ref(vertices_absolute)
+        geopoints = Geospatial_data(vertices,geo_reference = geo)
+        #print "vertices",vertices 
+        
+        d0 = [1.0, 1.0]
+        d1 = [1.0, 2.0]
+        d2 = [3.0, 1.0]
+        point_coords = [ d0, d1, d2]
+
+        interp = Interpolate(geopoints, triangles)
+        f = linear_function(vertices_absolute)
+        z = interp.interpolate(f, point_coords)
+        answer = linear_function(point_coords)
+
+        #print "z",z 
+        #print "answer",answer 
+        assert allclose(z, answer)
+        
     def test_interpolate_attributes_to_points(self):
         v0 = [0.0, 0.0]
@@ -309 +361 @@
         point_coords = [ d0, d1, d2]
 
-        interp = Interpolate(vertices, triangles, point_coords)
+        interp = Interpolate(vertices, triangles)
         f = linear_function(vertices)
         z = interp.interpolate(f, point_coords)
         answer = linear_function(point_coords)
 
-
+        #print "z",z 
+        #print "answer",answer 
         assert allclose(z, answer)