Changeset 976

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

Location:

inundation/ga/storm_surge/pyvolution

Files:

• general_mesh.py (modified) (2 diffs)
• least_squares.py (modified) (12 diffs)
• test_least_squares.py (modified) (1 diff)

inundation/ga/storm_surge/pyvolution/general_mesh.py

@@ -48 +48 @@
     """
 
-    def __init__(self, coordinates, triangles):
+    def __init__(self, coordinates, triangles, origin = None):
         """
         Build triangles from x,y coordinates (sequence of 2-tuples or
         Mx2 Numeric array of floats) and triangles (sequence of 3-tuples
         or Nx3 Numeric array of non-negative integers).
+
+        origin is a 3-tuple consisting of UTM zone, easting and northing.
+        If specified coordinates are assumed to be relative to this origin. 
         """
 
@@ -59 +62 @@
         self.triangles = array(triangles).astype(Int)
         self.coordinates = array(coordinates)
+
+        self.origin = origin
 
         #Input checks

inundation/ga/storm_surge/pyvolution/least_squares.py

@@ -59 +59 @@
 def fit_to_mesh_file(mesh_file, point_file, mesh_output_file,
                      alpha=DEFAULT_ALPHA, verbose= False,
-                     expand_search = False):
+                     expand_search = False,
+                     origin = None,
+                     mesh_origin=None):
     """
     Given a mesh file (tsh) and a point attribute file (xya), fit
     point attributes to the mesh and write a mesh file with the
     results.
+    
+    If origin is not None it is assumed to be
+    a 5-tuple with geo referenced
+    UTM coordinates (zone, easting, northing, false_easting, false_northing)
+    
+    mesh_origin is the same but refers to the input tsh file.
+    FIXME: When that format contains it own origin, this parameter can go.
     """
+    
     from load_mesh.loadASCII import mesh_file_to_mesh_dictionary, \
                  load_points_file, export_triangulation_file, \
                  concatinate_attributelist
+    
     # load in the .tsh file
     mesh_dict = mesh_file_to_mesh_dictionary(mesh_file)
@@ -78 +89 @@
    
     # load in the .xya file
+    #FIXME (Ole): Maybe data origin should be extracted here
     try:
         point_dict = load_points_file(point_file,delimiter = ',')
@@ -84 +96 @@
     point_coordinates = point_dict['pointlist']
     title_list,point_attributes = concatinate_attributelist(point_dict['attributelist'])  
-    if verbose:print "points file loaded"
+    if verbose: print "points file loaded"
+
+
+    
     if verbose:print "fitting to mesh"
     f = fit_to_mesh(vertex_coordinates,
@@ -216 +231 @@
                  verbose = False,
                  expand_search = True,
-                 max_points_per_cell = 30):
+                 max_points_per_cell = 30,
+                 mesh_origin = None, 
+                 data_origin = None):
 
         
@@ -236 +253 @@
 
           alpha: Smoothing parameter
+
+          data_origin and mesh_origin are 3-tuples consisting of
+          UTM zone, easting and northing. If specified
+          point coordinates and vertex coordinates are assumed to be
+          relative to their respective origins. 
           
         """
 
         #Convert input to Numeric arrays
+        triangles = array(triangles).astype(Int)
         vertex_coordinates = array(vertex_coordinates).astype(Float)
-        triangles = array(triangles).astype(Int)                
-        
+
         #Build underlying mesh
         if verbose: print 'Building mesh' 
-        self.mesh = General_mesh(vertex_coordinates, triangles)
+        self.mesh = General_mesh(vertex_coordinates, triangles,
+                                 origin = mesh_origin)
+        self.data_origin = data_origin
 
         #Smoothing parameter
@@ -254 +278 @@
                                         verbose = verbose,
                                         expand_search = expand_search,
-                                   max_points_per_cell = max_points_per_cell)
-        
-
-    def set_point_coordinates(self, point_coordinates):
+                                        max_points_per_cell =\
+                                        max_points_per_cell,
+                                        data_origin = data_origin)
+        
+
+    def set_point_coordinates(self, point_coordinates,
+                              data_origin = None):
         """
         A public interface to setting the point co-ordinates.
         """
-        self.build_coefficient_matrix_B(point_coordinates)
+        self.build_coefficient_matrix_B(point_coordinates, data_origin)
         
     def build_coefficient_matrix_B(self, point_coordinates=None,
                                    verbose = False, expand_search = True,
-                                   max_points_per_cell=30):
+                                   max_points_per_cell=30,
+                                   data_origin = None):
         """Build final coefficient matrix"""
         
@@ -273 +301 @@
             self.build_smoothing_matrix_D()
         
-        if point_coordinates:
+        if point_coordinates is not None:
+            
             if verbose: print 'Building interpolation matrix'         
             self.build_interpolation_matrix_A(point_coordinates,
                                               verbose = verbose,
                                               expand_search = expand_search,
-                                   max_points_per_cell = max_points_per_cell)
+                                              max_points_per_cell =\
+                                              max_points_per_cell,
+                                              data_origin = data_origin)
 
             if self.alpha <> 0:
@@ -290 +321 @@
     def build_interpolation_matrix_A(self, point_coordinates,
                                      verbose = False, expand_search = True,
-                                     max_points_per_cell=30):
+                                     max_points_per_cell=30,
+                                     data_origin = None):
         """Build n x m interpolation matrix, where
         n is the number of data points and
@@ -296 +328 @@
 
         This algorithm uses a quad tree data structure for fast binning of data points
+        origin is a 3-tuple consisting of UTM zone, easting and northing.
+        If specified coordinates are assumed to be relative to this origin. 
+
+        This one will override any data_origin that may be specified in
+        interpolation instance 
+
         """
 
         from quad import build_quadtree
-        
-        #Convert input to Numeric arrays
+
+        if data_origin is None:
+            data_origin = self.data_origin #Use the one from
+                                           #interpolation instance
+        
+        #Convert input to Numeric arrays just in case. 
         point_coordinates = array(point_coordinates).astype(Float)
+
+
+        #Shift data points to same origin as mesh (if specified)
+        mesh_origin = self.mesh.origin
+        if point_coordinates is not None:
+            if data_origin is not None:
+                if mesh_origin is not None:
+
+                    #Transformation:
+                    #
+                    #Let x_0 be the reference point of the point coordinates
+                    #and xi_0 the reference point of the mesh.
+                    #
+                    #A point coordinate (x + x_0) is then made relative
+                    #to xi_0 by
+                    #
+                    # x_new = x + x_0 - xi_0
+                    #
+                    #and similarly for eta
+                    
+                    x_offset = data_origin[1] - mesh_origin[1]
+                    y_offset = data_origin[2] - mesh_origin[2]
+                else: #Shift back to a zero origin
+                    x_offset = data_origin[1]
+                    y_offset = data_origin[2]                
+                    
+                point_coordinates[:,0] += x_offset
+                point_coordinates[:,1] += y_offset            
+            else:
+                if mesh_origin is not None:
+                    #Use mesh origin for data points
+                    point_coordinates[:,0] -= mesh_origin[1] 
+                    point_coordinates[:,1] -= mesh_origin[2]
+
+
+
+
         
         #Build n x m interpolation matrix        
@@ -307 +386 @@
         n = point_coordinates.shape[0]     #Nbr of data points
 
+        if verbose: print 'Number of datapoints: %d' %n
+        if verbose: print 'Number of basis functions: %d' %m
         
         #FIXME (Ole): We should use CSR here since mat-mat mult is now OK.
@@ -322 +403 @@
             #For each data_coordinate point
 
-            #if verbose: print 'Doing %d of %d' %(i, n)
+            if verbose and i%((n+10)/10)==0: print 'Doing %d of %d' %(i, n)
 
             x = point_coordinates[i]
@@ -656 +737 @@
             #Build n x m interpolation matrix        
             m = self.mesh.coordinates.shape[0] #Number of vertices
-            n = z.shape[1]               #Number of data points         
+            n = z.shape[1]                     #Number of data points         
 
             f = zeros((m,n), Float) #Resulting columns

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):