Changeset 2939

Timestamp:

May 22, 2006, 2:15:24 PM (19 years ago)

Author:

duncan

Message:

adding fit_to_mesh_file

Location:

inundation/fit_interpolate

Files:

• fit.py (modified) (5 diffs)
• test_fit.py (modified) (1 diff)

inundation/fit_interpolate/fit.py

@@ -21 +21 @@
 """
 
-from Numeric import zeros, Float
+from Numeric import zeros, Float, ArrayType
 
 from geospatial_data.geospatial_data import Geospatial_data, ensure_absolute
@@ -215 +215 @@
         data points.
 
+        If Ata is None, the matrices AtA and Atz are created.
+
+        This function can be called again and again, with sub-sets of
+        the point coordinates.  Call fit to get the results.
+        
         Preconditions
         z and points are numeric
@@ -235 +240 @@
 
             self.AtA = Sparse(m,m)
-            #self.Atz = zeros((m,att_num), Float)
         self.point_count += point_coordinates.shape[0]
         #print "_build_matrix_AtA_Atz - self.point_count", self.point_count
@@ -256 +260 @@
             element_found, sigma0, sigma1, sigma2, k = \
                            search_tree_of_vertices(self.root, self.mesh, x)
-            #Update interpolation matrix A if necessary
+            
             if element_found is True:
-                #Assign values to matrix A
-
                 j0 = self.mesh.triangles[k,0] #Global vertex id for sigma0
                 j1 = self.mesh.triangles[k,1] #Global vertex id for sigma1
@@ -431 +433 @@
 
 
+def fit_to_mesh_file(mesh_file, point_file, mesh_output_file,
+                     alpha=DEFAULT_ALPHA, verbose= False,
+                     display_errors = True):
+    """
+    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 data_origin is not None it is assumed to be
+    a 3-tuple with geo referenced
+    UTM coordinates (zone, easting, northing)
+
+    NOTE: Throws IOErrors, for a variety of file problems.
+    
+    """
+
+    # Question
+    # should data_origin and mesh_origin be passed in?
+    # No they should be in the data structure
+    #
+    #Should the origin of the mesh be changed using this function?
+    # That is overloading this function.  Have it as a seperate
+    # method, at least initially.
+    
+    from load_mesh.loadASCII import import_mesh_file, \
+                 import_points_file, export_mesh_file, \
+                 concatinate_attributelist
+
+    # FIXME: Use geospatial instead of import_points_file
+    try:
+        mesh_dict = import_mesh_file(mesh_file)
+    except IOError,e:
+        if display_errors:
+            print "Could not load bad file. ", e
+        raise IOError  #Re-raise exception
+        
+    vertex_coordinates = mesh_dict['vertices']
+    triangles = mesh_dict['triangles']
+    if type(mesh_dict['vertex_attributes']) == ArrayType:
+        old_point_attributes = mesh_dict['vertex_attributes'].tolist()
+    else:
+        old_point_attributes = mesh_dict['vertex_attributes']
+
+    if type(mesh_dict['vertex_attribute_titles']) == ArrayType:
+        old_title_list = mesh_dict['vertex_attribute_titles'].tolist()
+    else:
+        old_title_list = mesh_dict['vertex_attribute_titles']
+
+    if verbose: print 'tsh file %s loaded' %mesh_file
+
+    # load in the .pts file
+    try:
+        point_dict = import_points_file(point_file, verbose=verbose)
+    except IOError,e:
+        if display_errors:
+            print "Could not load bad file. ", e
+        raise IOError  #Re-raise exception  
+
+    point_coordinates = point_dict['pointlist']
+    title_list,point_attributes = concatinate_attributelist(point_dict['attributelist'])
+
+    if point_dict.has_key('geo_reference') and not point_dict['geo_reference'] is None:
+        data_origin = point_dict['geo_reference'].get_origin()
+    else:
+        data_origin = (56, 0, 0) #FIXME(DSG-DSG)
+
+    if mesh_dict.has_key('geo_reference') and not mesh_dict['geo_reference'] is None:
+        mesh_origin = mesh_dict['geo_reference'].get_origin()
+    else:
+        mesh_origin = (56, 0, 0) #FIXME(DSG-DSG)
+
+    if verbose: print "points file loaded"
+    if verbose: print "fitting to mesh"
+    f = fit_to_mesh(vertex_coordinates,
+                    triangles,
+                    point_coordinates,
+                    point_attributes,
+                    alpha = alpha,
+                    verbose = verbose,
+                    data_origin = data_origin,
+                    mesh_origin = mesh_origin)
+    if verbose: print "finished fitting to mesh"
+
+    # convert array to list of lists
+    new_point_attributes = f.tolist()
+    #FIXME have this overwrite attributes with the same title - DSG
+    #Put the newer attributes last
+    if old_title_list <> []:
+        old_title_list.extend(title_list)
+        #FIXME can this be done a faster way? - DSG
+        for i in range(len(old_point_attributes)):
+            old_point_attributes[i].extend(new_point_attributes[i])
+        mesh_dict['vertex_attributes'] = old_point_attributes
+        mesh_dict['vertex_attribute_titles'] = old_title_list
+    else:
+        mesh_dict['vertex_attributes'] = new_point_attributes
+        mesh_dict['vertex_attribute_titles'] = title_list
+
+    if verbose: print "exporting to file ", mesh_output_file
+
+    try:
+        export_mesh_file(mesh_output_file, mesh_dict)
+    except IOError,e:
+        if display_errors:
+            print "Could not write file. ", e
+        raise IOError
+
+
 def _fit(*args, **kwargs):
     """Private function for use with caching. Reason is that classes

inundation/fit_interpolate/test_fit.py

@@ -488 +488 @@
 
 
+    def test_fit_to_mesh_file(self):
+        from load_mesh.loadASCII import import_mesh_file, \
+             export_mesh_file
+        import tempfile
+        import os
+
+        # create a .tsh file, no user outline
+        mesh_dic = {}
+        mesh_dic['vertices'] = [[0.0, 0.0],
+                                          [0.0, 5.0],
+                                          [5.0, 0.0]]
+        mesh_dic['triangles'] =  [[0, 2, 1]]
+        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
+        mesh_dic['triangle_tags'] = ['']
+        mesh_dic['vertex_attributes'] = [[], [], []]
+        mesh_dic['vertiex_attribute_titles'] = []
+        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
+        mesh_dic['segment_tags'] = ['external',
+                                                  'external',
+                                                  'external']
+        mesh_file = tempfile.mktemp(".tsh")
+        export_mesh_file(mesh_file,mesh_dic)
+
+        # create an .xya file
+        point_file = tempfile.mktemp(".xya")
+        fd = open(point_file,'w')
+        fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
+        fd.close()
+
+        mesh_output_file = tempfile.mktemp(".tsh") 
+        fit_to_mesh_file(mesh_file,
+                         point_file,
+                         mesh_output_file,
+                         alpha = 0.0)
+        # load in the .tsh file we just wrote
+        mesh_dic = import_mesh_file(mesh_output_file)
+        #print "mesh_dic",mesh_dic
+        ans =[[0.0, 0.0],
+              [5.0, 10.0],
+              [5.0,10.0]]
+        assert allclose(mesh_dic['vertex_attributes'],ans)
+
+        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
+                        ['elevation','stage'],
+                        'test_fit_to_mesh_file failed')
+
+        #clean up
+        os.remove(mesh_file)
+        os.remove(point_file)
+        os.remove(mesh_output_file)
+
+
+    def test_fit_to_mesh_file3(self):
+        from load_mesh.loadASCII import import_mesh_file, \
+             export_mesh_file
+        import tempfile
+        import os
+
+        # create a .tsh file, no user outline
+        mesh_dic = {}
+        mesh_dic['vertices'] = [[0.76, 0.76],
+                                          [0.76, 5.76],
+                                          [5.76, 0.76]]
+        mesh_dic['triangles'] =  [[0, 2, 1]]
+        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
+        mesh_dic['triangle_tags'] = ['']
+        mesh_dic['vertex_attributes'] = [[], [], []]
+        mesh_dic['vertiex_attribute_titles'] = []
+        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
+        mesh_dic['segment_tags'] = ['external',
+                                                  'external',
+                                                  'external']
+        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
+        mesh_file = tempfile.mktemp(".tsh")
+        export_mesh_file(mesh_file,mesh_dic)
+
+        # create an .xya file
+        point_file = tempfile.mktemp(".xya")
+        fd = open(point_file,'w')
+        fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
+        fd.close()
+
+        mesh_output_file = tempfile.mktemp(".tsh")
+        fit_to_mesh_file(mesh_file,
+                         point_file,
+                         mesh_output_file,
+                         alpha = 0.0)
+        # load in the .tsh file we just wrote
+        mesh_dic = import_mesh_file(mesh_output_file)
+        #print "mesh_dic",mesh_dic
+        ans =[[0.0, 0.0],
+              [5.0, 10.0],
+              [5.0,10.0]]
+        assert allclose(mesh_dic['vertex_attributes'],ans)
+
+        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
+                        ['elevation','stage'],
+                        'test_fit_to_mesh_file failed')
+
+        #clean up
+        os.remove(mesh_file)
+        os.remove(point_file)
+        os.remove(mesh_output_file)
+
+    def test_fit_to_mesh_file4(self):
+        from load_mesh.loadASCII import import_mesh_file, \
+             export_mesh_file
+        import tempfile
+        import os
+
+        # create a .tsh file, no user outline
+        mesh_dic = {}
+        mesh_dic['vertices'] = [[0.76, 0.76],
+                                [0.76, 5.76],
+                                [5.76, 0.76]]
+        mesh_dic['triangles'] =  [[0, 2, 1]]
+        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
+        mesh_dic['triangle_tags'] = ['']
+        mesh_dic['vertex_attributes'] = [[], [], []]
+        mesh_dic['vertiex_attribute_titles'] = []
+        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
+        mesh_dic['segment_tags'] = ['external',
+                                    'external',
+                                    'external']
+        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
+        mesh_file = tempfile.mktemp(".tsh")
+        export_mesh_file(mesh_file,mesh_dic)
+
+        geo_ref = Geo_reference(56,-200,-400)
+        # create an .xya file
+        point_file = tempfile.mktemp(".xya")
+        fd = open(point_file,'w')
+        fd.write("elevation, stage \n 201.0, 401.0,2.,4 \n 201.0, 403.0,4,8 \n 203.0, 401.0,4.,8 \n")
+        geo_ref.write_ASCII(fd)
+        fd.close()
+
+        mesh_output_file = tempfile.mktemp(".tsh")
+        fit_to_mesh_file(mesh_file,
+                         point_file,
+                         mesh_output_file,
+                         alpha = 0.0)
+        # load in the .tsh file we just wrote
+        mesh_dic = import_mesh_file(mesh_output_file)
+        #print "mesh_dic",mesh_dic
+        ans =[[0.0, 0.0],
+              [5.0, 10.0],
+              [5.0, 10.0]]
+        assert allclose(mesh_dic['vertex_attributes'],ans)
+
+        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
+                        ['elevation','stage'],
+                        'test_fit_to_mesh_file failed')
+
+        #clean up
+        os.remove(mesh_file)
+        os.remove(point_file)
+        os.remove(mesh_output_file)
+
+    def test_fit_to_mesh_fileII(self):
+        from load_mesh.loadASCII import import_mesh_file, \
+             export_mesh_file
+        import tempfile
+        import os
+
+        # create a .tsh file, no user outline
+        mesh_dic = {}
+        mesh_dic['vertices'] = [[0.0, 0.0],
+                                [0.0, 5.0],
+                                [5.0, 0.0]]
+        mesh_dic['triangles'] =  [[0, 2, 1]]
+        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
+        mesh_dic['triangle_tags'] = ['']
+        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
+        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
+        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
+        mesh_dic['segment_tags'] = ['external',
+                                                  'external',
+                                                  'external']
+        mesh_file = tempfile.mktemp(".tsh")
+        export_mesh_file(mesh_file,mesh_dic)
+
+        # create an .xya file
+        point_file = tempfile.mktemp(".xya")
+        fd = open(point_file,'w')
+        fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
+        fd.close()
+
+        mesh_output_file = "new_triangle.tsh"
+        fit_to_mesh_file(mesh_file,
+                         point_file,
+                         mesh_output_file,
+                         alpha = 0.0)
+        # load in the .tsh file we just wrote
+        mesh_dic = import_mesh_file(mesh_output_file)
+
+        assert allclose(mesh_dic['vertex_attributes'],
+                        [[1.0, 2.0,0.0, 0.0],
+                         [1.0, 2.0,5.0, 10.0],
+                         [1.0, 2.0,5.0,10.0]])
+
+        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
+                        ['density', 'temp','elevation','stage'],
+                        'test_fit_to_mesh_file failed')
+
+        #clean up
+        os.remove(mesh_file)
+        os.remove(mesh_output_file)
+        os.remove(point_file)
+
+    def test_fit_to_mesh_file_errors(self):
+        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
+        import tempfile
+        import os
+
+        # create a .tsh file, no user outline
+        mesh_dic = {}
+        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
+        mesh_dic['triangles'] =  [[0, 2, 1]]
+        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
+        mesh_dic['triangle_tags'] = ['']
+        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
+        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
+        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
+        mesh_dic['segment_tags'] = ['external', 'external','external']
+        mesh_file = tempfile.mktemp(".tsh")
+        export_mesh_file(mesh_file,mesh_dic)
+
+        # create an .xya file
+        point_file = tempfile.mktemp(".xya")
+        fd = open(point_file,'w')
+        fd.write("elevation stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
+        fd.close()
+
+        mesh_output_file = "new_triangle.tsh"
+        try:
+            fit_to_mesh_file(mesh_file, point_file,
+                             mesh_output_file, display_errors = False)
+        except IOError:
+            pass
+        else:
+            #self.failUnless(0 ==1,  'Bad file did not raise error!')
+            raise 'Bad file did not raise error!'
+            
+        #clean up
+        os.remove(mesh_file)
+        os.remove(point_file)
+
+    def test_fit_to_mesh_file_errorsII(self):
+        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
+        import tempfile
+        import os
+
+        # create a .tsh file, no user outline
+        mesh_file = tempfile.mktemp(".tsh")
+        fd = open(mesh_file,'w')
+        fd.write("unit testing a bad .tsh file \n")
+        fd.close()
+
+        # create an .xya file
+        point_file = tempfile.mktemp(".xya")
+        fd = open(point_file,'w')
+        fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
+        fd.close()
+
+        mesh_output_file = "new_triangle.tsh"
+        try:
+            fit_to_mesh_file(mesh_file, point_file,
+                             mesh_output_file, display_errors = False)
+        except IOError:
+            pass
+        else:
+            raise 'Bad file did not raise error!'
+            
+        #clean up
+        os.remove(mesh_file)
+        os.remove(point_file)
+
+    def test_fit_to_mesh_file_errorsIII(self):
+        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
+        import tempfile
+        import os
+
+        # create a .tsh file, no user outline
+        mesh_dic = {}
+        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
+        mesh_dic['triangles'] =  [[0, 2, 1]]
+        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
+        mesh_dic['triangle_tags'] = ['']
+        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
+        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
+        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
+        mesh_dic['segment_tags'] = ['external', 'external','external']
+        mesh_file = tempfile.mktemp(".tsh")
+        export_mesh_file(mesh_file,mesh_dic)
+
+        # create an .xya file
+        point_file = tempfile.mktemp(".xya")
+        fd = open(point_file,'w')
+        fd.write("elevation, stage \n 1.0, 1.0,2.,4 \n 1.0, 3.0,4,8 \n 3.0,1.0,4.,8 \n")
+        fd.close()
+
+        #This a deliberately illegal filename to invoke the error.
+        mesh_output_file = ".../\z\z:ya.tsh"        
+
+        try:
+            fit_to_mesh_file(mesh_file, point_file,
+                             mesh_output_file, display_errors = False)
+        except IOError:
+            pass
+        else:
+            raise 'Bad file did not raise error!'
+        
+        #clean up
+        os.remove(mesh_file)
+        os.remove(point_file)
 
 #-------------------------------------------------------------