Changeset 4535

Timestamp:

Jun 6, 2007, 5:27:22 PM (18 years ago)

Author:

duncan

Message:

getting fit 2 mesh file going again

Location:

anuga_core/source/anuga

Files:

• fit_interpolate/fit.py (modified) (2 diffs)
• fit_interpolate/test_fit.py (modified) (2 diffs)
• geospatial_data/geospatial_data.py (modified) (4 diffs)
• geospatial_data/test_geospatial_data.py (modified) (2 diffs)
• shallow_water/data_manager.py (modified) (1 diff)

anuga_core/source/anuga/fit_interpolate/fit.py

@@ -522 +522 @@
                         'mesh_origin': mesh_origin,
                         'alpha':alpha},
+                       compression = False,
                        verbose = verbose)        
         
@@ -554 +555 @@
     return Fit(*args, **kwargs)
 
+
+def fit_to_mesh_file(mesh_file, point_file, mesh_output_file,
+                     alpha=DEFAULT_ALPHA, verbose= False,
+                     expand_search = False,
+                     precrop = 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.
+
+    Note: the .xya files need titles.  If you want anuga to use the tsh file,
+    make sure the title is elevation.
+
+    NOTE: Throws IOErrors, for a variety of file problems.
+    
+    """
+
+    from load_mesh.loadASCII import import_mesh_file, \
+         export_mesh_file, concatinate_attributelist
+
+
+    try:
+        mesh_dict = import_mesh_file(mesh_file)
+    except IOError,e:
+        if display_errors:
+            print "Could not load bad file. ", e
+        raise IOError  #Could not load bad mesh file.
+    
+    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)
+        geo = Geospatial_data(point_file, verbose=verbose)
+    except IOError,e:
+        if display_errors:
+            print "Could not load bad file. ", e
+        raise IOError  #Re-raise exception  
+
+    point_coordinates = geo.get_data_points(absolute=True)
+    #point_dict['pointlist']
+    #
+    # return list of attribute titles, array of attributes 
+    # title_list,point_attributes = concatinate_attributelist(point_dict['attributelist'])
+    title_list,point_attributes = concatinate_attributelist(geo.get_all_attributes())
+
+    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 = None #(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 = None,
+                    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
+
+    #FIXME (Ole): Remember to output mesh_origin as well
+    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

anuga_core/source/anuga/fit_interpolate/test_fit.py

@@ -16 +16 @@
 from anuga.utilities.numerical_tools import ensure_numeric
 from anuga.geospatial_data.geospatial_data import Geospatial_data
+from anuga.shallow_water import Domain
 
 def distance(x, y):
@@ -728 +729 @@
 
 
+    def test_fit_to_mesh_file2domain(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')
+        domain = Domain(mesh_output_file, use_cache=True, verbose=False)
+        
+        answer = [0., 5., 5.]
+        assert allclose(domain.quantities['elevation'].vertex_values,
+                        answer)
+        #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)
+  
     def Not_yet_test_smooth_att_to_mesh_with_excess_verts(self):
 

anuga_core/source/anuga/geospatial_data/geospatial_data.py

@@ -24 +24 @@
 from anuga.utilities.anuga_exceptions import ANUGAError
 from anuga.config import points_file_block_line_size as MAX_READ_LINES    
+
+DEFAULT_ATTRIBUTE = 'elevation'
 
 class Geospatial_data:
@@ -53 +55 @@
         of lists (or arrays) of length M. In the latter case the keys
         in the dictionary represent the attribute names, in the former
-        the attribute will get the default name "attribute".
+        the attribute will get the default name "elevation".
         
         geo_reference: Object representing the origin of the data
@@ -208 +210 @@
         if not isinstance(attributes, DictType):
             #Convert single attribute into dictionary
-            attributes = {'attribute': attributes}
+            attributes = {DEFAULT_ATTRIBUTE: attributes}
 
         #Check input attributes    
@@ -514 +516 @@
         attributes = {}
         if file_name[-4:]== ".xya":
-            msg = 'Text file format is moving to comma seperated .txt files.'
-            warn(msg, DeprecationWarning) 
+            # Maybe not phase-out, so we can load in geo-ref info
+            #msg = 'Text file format is moving to comma seperated .txt files.'
+            #warn(msg, DeprecationWarning) 
             try:
                 if delimiter == None:

anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

@@ -84 +84 @@
         points = [[1.0, 2.1], [3.0, 5.3]]
         attributes = [2, 4]
-        G = Geospatial_data(points, attributes)        
-
-        assert G.attributes.keys()[0] == 'attribute'
+        G = Geospatial_data(points, attributes)       
+        assert G.attributes.keys()[0] == DEFAULT_ATTRIBUTE
         assert allclose(G.attributes.values()[0], [2, 4])
         
@@ -119 +118 @@
         assert allclose(V, [2, 4])
 
-        V = G.get_attributes('attribute') #Get by name
+        V = G.get_attributes(DEFAULT_ATTRIBUTE) #Get by name
         assert allclose(V, [2, 4])
 

anuga_core/source/anuga/shallow_water/data_manager.py

@@ -1636 +1636 @@
 
     This function returns the names of the files produced.
+
+    It will also produce as many output files as there are input sww files. 
     """