Changeset 976 for inundation/ga/storm_surge/pyvolution/test_least_squares.py

Timestamp:

Mar 1, 2005, 12:16:51 PM (20 years ago)

Author:

ole

Message:

Made least squares work with different coordinate systems for data p;oints and mesh

File:

• inundation/ga/storm_surge/pyvolution/test_least_squares.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/test_least_squares.py

@@ -878 +878 @@
 
 
+    def test_fit_and_interpolation_with_different_origins(self):
+        """Fit a surface to one set of points. Then interpolate that surface
+        using another set of points.
+        This test tests situtaion where points and mesh belong to a different
+        coordinate system as defined by origin.
+        """
+        from mesh import Mesh
+
+        #Setup mesh used to represent fitted function
+        a = [0.0, 0.0]
+        b = [0.0, 2.0]
+        c = [2.0, 0.0]
+        d = [0.0, 4.0]
+        e = [2.0, 2.0]
+        f = [4.0, 0.0]
+
+        points = [a, b, c, d, e, f]
+        #bac, bce, ecf, dbe, daf, dae 
+        triangles = [[1,0,2], [1,2,4], [4,2,5], [3,1,4]]
+
+        #Datapoints to fit from
+        data_points1 = [[ 0.66666667, 0.66666667],
+                        [ 1.33333333, 1.33333333],
+                        [ 2.66666667, 0.66666667],
+                        [ 0.66666667, 2.66666667],
+                        [ 0.0, 1.0],
+                        [ 0.0, 3.0],
+                        [ 1.0, 0.0],
+                        [ 1.0, 1.0],
+                        [ 1.0, 2.0],
+                        [ 1.0, 3.0],                                          
+                        [ 2.0, 1.0],
+                        [ 3.0, 0.0],
+                        [ 3.0, 1.0]]
+
+
+        #First check that things are OK when using same origin
+        mesh_origin = (56, 290000, 618000) #zone, easting, northing
+        data_origin = (56, 290000, 618000) #zone, easting, northing 
+
+
+        #Fit surface to mesh
+        interp = Interpolation(points, triangles, data_points1,
+                               alpha=0.0,
+                               data_origin = data_origin,
+                               mesh_origin = mesh_origin)
+        
+        z = linear_function(data_points1) #Example z-values
+        f = interp.fit(z)                 #Fitted values at vertices
+
+
+        #New datapoints where interpolated values are sought
+        data_points2 = [[ 0.0, 0.0],
+                        [ 0.5, 0.5],
+                        [ 0.7, 0.7],
+                        [ 1.0, 0.5],
+                        [ 2.0, 0.4],
+                        [ 2.8, 1.2]]                                          
+        
+
+        #Build new A matrix based on new points
+        interp.build_interpolation_matrix_A(data_points2)
+
+        #Interpolate using fitted surface
+        z1 = interp.interpolate(f)
+
+        #Desired result
+        answer = linear_function(data_points2)
+        assert allclose(z1, answer)
+
+
+        ##############################################
+
+        #Then check situtaion where points are relative to a different
+        #origin (same zone, though, until we figure that out (FIXME))
+        
+        mesh_origin = (56, 290000, 618000) #zone, easting, northing
+        data_origin = (56, 10000, 10000) #zone, easting, northing 
+
+        #Shift datapoints according to new origin
+
+        for k in range(len(data_points1)):
+            data_points1[k][0] += mesh_origin[1] - data_origin[1]
+            data_points1[k][1] += mesh_origin[2] - data_origin[2]
+
+        for k in range(len(data_points2)):
+            data_points2[k][0] += mesh_origin[1] - data_origin[1]
+            data_points2[k][1] += mesh_origin[2] - data_origin[2]            
+
+
+        
+        #Fit surface to mesh
+        interp = Interpolation(points, triangles, data_points1,
+                               alpha=0.0,
+                               data_origin = data_origin,
+                               mesh_origin = mesh_origin)
+        
+        f1 = interp.fit(z) #Fitted values at vertices (using same z as before)
+
+        assert allclose(f,f1), 'Fit should have been unaltered'
+
+
+        #Build new A matrix based on new points
+        interp.build_interpolation_matrix_A(data_points2)
+
+        #Interpolate using fitted surface
+        z1 = interp.interpolate(f)
+        assert allclose(z1, answer)
+
+
 
     def test_fit_to_mesh_file(self):