Changeset 4535


Ignore:
Timestamp:
Jun 6, 2007, 5:27:22 PM (17 years ago)
Author:
duncan
Message:

getting fit 2 mesh file going again

Location:
anuga_core/source/anuga
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • anuga_core/source/anuga/fit_interpolate/fit.py

    r4425 r4535  
    522522                        'mesh_origin': mesh_origin,
    523523                        'alpha':alpha},
     524                       compression = False,
    524525                       verbose = verbose)       
    525526       
     
    554555    return Fit(*args, **kwargs)
    555556
     557
     558def fit_to_mesh_file(mesh_file, point_file, mesh_output_file,
     559                     alpha=DEFAULT_ALPHA, verbose= False,
     560                     expand_search = False,
     561                     precrop = False,
     562                     display_errors = True):
     563    """
     564    Given a mesh file (tsh) and a point attribute file (xya), fit
     565    point attributes to the mesh and write a mesh file with the
     566    results.
     567
     568    Note: the .xya files need titles.  If you want anuga to use the tsh file,
     569    make sure the title is elevation.
     570
     571    NOTE: Throws IOErrors, for a variety of file problems.
     572   
     573    """
     574
     575    from load_mesh.loadASCII import import_mesh_file, \
     576         export_mesh_file, concatinate_attributelist
     577
     578
     579    try:
     580        mesh_dict = import_mesh_file(mesh_file)
     581    except IOError,e:
     582        if display_errors:
     583            print "Could not load bad file. ", e
     584        raise IOError  #Could not load bad mesh file.
     585   
     586    vertex_coordinates = mesh_dict['vertices']
     587    triangles = mesh_dict['triangles']
     588    if type(mesh_dict['vertex_attributes']) == ArrayType:
     589        old_point_attributes = mesh_dict['vertex_attributes'].tolist()
     590    else:
     591        old_point_attributes = mesh_dict['vertex_attributes']
     592
     593    if type(mesh_dict['vertex_attribute_titles']) == ArrayType:
     594        old_title_list = mesh_dict['vertex_attribute_titles'].tolist()
     595    else:
     596        old_title_list = mesh_dict['vertex_attribute_titles']
     597
     598    if verbose: print 'tsh file %s loaded' %mesh_file
     599
     600    # load in the .pts file
     601    try:
     602        ###point_dict = import_points_file(point_file, verbose=verbose)
     603        geo = Geospatial_data(point_file, verbose=verbose)
     604    except IOError,e:
     605        if display_errors:
     606            print "Could not load bad file. ", e
     607        raise IOError  #Re-raise exception 
     608
     609    point_coordinates = geo.get_data_points(absolute=True)
     610    #point_dict['pointlist']
     611    #
     612    # return list of attribute titles, array of attributes
     613    # title_list,point_attributes = concatinate_attributelist(point_dict['attributelist'])
     614    title_list,point_attributes = concatinate_attributelist(geo.get_all_attributes())
     615
     616    if mesh_dict.has_key('geo_reference') and not mesh_dict['geo_reference'] is None:
     617        mesh_origin = mesh_dict['geo_reference'].get_origin()
     618    else:
     619        mesh_origin = None #(56, 0, 0) #FIXME(DSG-DSG)
     620
     621    if verbose: print "points file loaded"
     622    if verbose: print "fitting to mesh"
     623    f = fit_to_mesh(vertex_coordinates,
     624                    triangles,
     625                    point_coordinates,
     626                    point_attributes,
     627                    alpha = alpha,
     628                    verbose = verbose,
     629                    data_origin = None,
     630                    mesh_origin = mesh_origin)
     631    if verbose: print "finished fitting to mesh"
     632
     633    # convert array to list of lists
     634    new_point_attributes = f.tolist()
     635    #FIXME have this overwrite attributes with the same title - DSG
     636    #Put the newer attributes last
     637    if old_title_list <> []:
     638        old_title_list.extend(title_list)
     639        #FIXME can this be done a faster way? - DSG
     640        for i in range(len(old_point_attributes)):
     641            old_point_attributes[i].extend(new_point_attributes[i])
     642        mesh_dict['vertex_attributes'] = old_point_attributes
     643        mesh_dict['vertex_attribute_titles'] = old_title_list
     644    else:
     645        mesh_dict['vertex_attributes'] = new_point_attributes
     646        mesh_dict['vertex_attribute_titles'] = title_list
     647
     648    #FIXME (Ole): Remember to output mesh_origin as well
     649    if verbose: print "exporting to file ", mesh_output_file
     650
     651    try:
     652        export_mesh_file(mesh_output_file, mesh_dict)
     653    except IOError,e:
     654        if display_errors:
     655            print "Could not write file. ", e
     656        raise IOError
  • anuga_core/source/anuga/fit_interpolate/test_fit.py

    r4253 r4535  
    1616from anuga.utilities.numerical_tools import ensure_numeric
    1717from anuga.geospatial_data.geospatial_data import Geospatial_data
     18from anuga.shallow_water import Domain
    1819
    1920def distance(x, y):
     
    728729
    729730
     731    def test_fit_to_mesh_file2domain(self):
     732        from load_mesh.loadASCII import import_mesh_file, \
     733             export_mesh_file
     734        import tempfile
     735        import os
     736
     737        # create a .tsh file, no user outline
     738        mesh_dic = {}
     739        mesh_dic['vertices'] = [[0.0, 0.0],
     740                                          [0.0, 5.0],
     741                                          [5.0, 0.0]]
     742        mesh_dic['triangles'] =  [[0, 2, 1]]
     743        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
     744        mesh_dic['triangle_tags'] = ['']
     745        mesh_dic['vertex_attributes'] = [[], [], []]
     746        mesh_dic['vertiex_attribute_titles'] = []
     747        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
     748        mesh_dic['segment_tags'] = ['external',
     749                                                  'external',
     750                                                  'external']
     751        mesh_file = tempfile.mktemp(".tsh")
     752        export_mesh_file(mesh_file,mesh_dic)
     753
     754        # create an .xya file
     755        point_file = tempfile.mktemp(".xya")
     756        fd = open(point_file,'w')
     757        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")
     758        fd.close()
     759
     760        mesh_output_file = tempfile.mktemp(".tsh")
     761        fit_to_mesh_file(mesh_file,
     762                         point_file,
     763                         mesh_output_file,
     764                         alpha = 0.0)
     765        # load in the .tsh file we just wrote
     766        mesh_dic = import_mesh_file(mesh_output_file)
     767        #print "mesh_dic",mesh_dic
     768        ans =[[0.0, 0.0],
     769              [5.0, 10.0],
     770              [5.0,10.0]]
     771        assert allclose(mesh_dic['vertex_attributes'],ans)
     772
     773        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
     774                        ['elevation','stage'],
     775                        'test_fit_to_mesh_file failed')
     776        domain = Domain(mesh_output_file, use_cache=True, verbose=False)
     777       
     778        answer = [0., 5., 5.]
     779        assert allclose(domain.quantities['elevation'].vertex_values,
     780                        answer)
     781        #clean up
     782        os.remove(mesh_file)
     783        os.remove(point_file)
     784        os.remove(mesh_output_file)
     785
     786    def test_fit_to_mesh_file3(self):
     787        from load_mesh.loadASCII import import_mesh_file, \
     788             export_mesh_file
     789        import tempfile
     790        import os
     791
     792        # create a .tsh file, no user outline
     793        mesh_dic = {}
     794        mesh_dic['vertices'] = [[0.76, 0.76],
     795                                          [0.76, 5.76],
     796                                          [5.76, 0.76]]
     797        mesh_dic['triangles'] =  [[0, 2, 1]]
     798        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
     799        mesh_dic['triangle_tags'] = ['']
     800        mesh_dic['vertex_attributes'] = [[], [], []]
     801        mesh_dic['vertiex_attribute_titles'] = []
     802        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
     803        mesh_dic['segment_tags'] = ['external',
     804                                                  'external',
     805                                                  'external']
     806        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
     807        mesh_file = tempfile.mktemp(".tsh")
     808        export_mesh_file(mesh_file,mesh_dic)
     809
     810        # create an .xya file
     811        point_file = tempfile.mktemp(".xya")
     812        fd = open(point_file,'w')
     813        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")
     814        fd.close()
     815
     816        mesh_output_file = tempfile.mktemp(".tsh")
     817        fit_to_mesh_file(mesh_file,
     818                         point_file,
     819                         mesh_output_file,
     820                         alpha = 0.0)
     821        # load in the .tsh file we just wrote
     822        mesh_dic = import_mesh_file(mesh_output_file)
     823        #print "mesh_dic",mesh_dic
     824        ans =[[0.0, 0.0],
     825              [5.0, 10.0],
     826              [5.0,10.0]]
     827        assert allclose(mesh_dic['vertex_attributes'],ans)
     828
     829        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
     830                        ['elevation','stage'],
     831                        'test_fit_to_mesh_file failed')
     832
     833        #clean up
     834        os.remove(mesh_file)
     835        os.remove(point_file)
     836        os.remove(mesh_output_file)
     837
     838    def test_fit_to_mesh_file4(self):
     839        from load_mesh.loadASCII import import_mesh_file, \
     840             export_mesh_file
     841        import tempfile
     842        import os
     843
     844        # create a .tsh file, no user outline
     845        mesh_dic = {}
     846        mesh_dic['vertices'] = [[0.76, 0.76],
     847                                [0.76, 5.76],
     848                                [5.76, 0.76]]
     849        mesh_dic['triangles'] =  [[0, 2, 1]]
     850        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
     851        mesh_dic['triangle_tags'] = ['']
     852        mesh_dic['vertex_attributes'] = [[], [], []]
     853        mesh_dic['vertiex_attribute_titles'] = []
     854        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
     855        mesh_dic['segment_tags'] = ['external',
     856                                    'external',
     857                                    'external']
     858        mesh_dic['geo_reference'] = Geo_reference(56,-0.76,-0.76)
     859        mesh_file = tempfile.mktemp(".tsh")
     860        export_mesh_file(mesh_file,mesh_dic)
     861
     862        geo_ref = Geo_reference(56,-200,-400)
     863        # create an .xya file
     864        point_file = tempfile.mktemp(".xya")
     865        fd = open(point_file,'w')
     866        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")
     867        geo_ref.write_ASCII(fd)
     868        fd.close()
     869
     870        mesh_output_file = tempfile.mktemp(".tsh")
     871        fit_to_mesh_file(mesh_file,
     872                         point_file,
     873                         mesh_output_file,
     874                         alpha = 0.0)
     875        # load in the .tsh file we just wrote
     876        mesh_dic = import_mesh_file(mesh_output_file)
     877        #print "mesh_dic",mesh_dic
     878        ans =[[0.0, 0.0],
     879              [5.0, 10.0],
     880              [5.0, 10.0]]
     881        assert allclose(mesh_dic['vertex_attributes'],ans)
     882
     883        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
     884                        ['elevation','stage'],
     885                        'test_fit_to_mesh_file failed')
     886
     887        #clean up
     888        os.remove(mesh_file)
     889        os.remove(point_file)
     890        os.remove(mesh_output_file)
     891
     892    def test_fit_to_mesh_fileII(self):
     893        from load_mesh.loadASCII import import_mesh_file, \
     894             export_mesh_file
     895        import tempfile
     896        import os
     897
     898        # create a .tsh file, no user outline
     899        mesh_dic = {}
     900        mesh_dic['vertices'] = [[0.0, 0.0],
     901                                [0.0, 5.0],
     902                                [5.0, 0.0]]
     903        mesh_dic['triangles'] =  [[0, 2, 1]]
     904        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
     905        mesh_dic['triangle_tags'] = ['']
     906        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
     907        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
     908        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
     909        mesh_dic['segment_tags'] = ['external',
     910                                                  'external',
     911                                                  'external']
     912        mesh_file = tempfile.mktemp(".tsh")
     913        export_mesh_file(mesh_file,mesh_dic)
     914
     915        # create an .xya file
     916        point_file = tempfile.mktemp(".xya")
     917        fd = open(point_file,'w')
     918        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")
     919        fd.close()
     920
     921        mesh_output_file = "new_triangle.tsh"
     922        fit_to_mesh_file(mesh_file,
     923                         point_file,
     924                         mesh_output_file,
     925                         alpha = 0.0)
     926        # load in the .tsh file we just wrote
     927        mesh_dic = import_mesh_file(mesh_output_file)
     928
     929        assert allclose(mesh_dic['vertex_attributes'],
     930                        [[1.0, 2.0,0.0, 0.0],
     931                         [1.0, 2.0,5.0, 10.0],
     932                         [1.0, 2.0,5.0,10.0]])
     933
     934        self.failUnless(mesh_dic['vertex_attribute_titles']  ==
     935                        ['density', 'temp','elevation','stage'],
     936                        'test_fit_to_mesh_file failed')
     937
     938        #clean up
     939        os.remove(mesh_file)
     940        os.remove(mesh_output_file)
     941        os.remove(point_file)
     942
     943    def test_fit_to_mesh_file_errors(self):
     944        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
     945        import tempfile
     946        import os
     947
     948        # create a .tsh file, no user outline
     949        mesh_dic = {}
     950        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
     951        mesh_dic['triangles'] =  [[0, 2, 1]]
     952        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
     953        mesh_dic['triangle_tags'] = ['']
     954        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
     955        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
     956        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
     957        mesh_dic['segment_tags'] = ['external', 'external','external']
     958        mesh_file = tempfile.mktemp(".tsh")
     959        export_mesh_file(mesh_file,mesh_dic)
     960
     961        # create an .xya file
     962        point_file = tempfile.mktemp(".xya")
     963        fd = open(point_file,'w')
     964        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")
     965        fd.close()
     966
     967        mesh_output_file = "new_triangle.tsh"
     968        try:
     969            fit_to_mesh_file(mesh_file, point_file,
     970                             mesh_output_file, display_errors = False)
     971        except IOError:
     972            pass
     973        else:
     974            #self.failUnless(0 ==1,  'Bad file did not raise error!')
     975            raise 'Bad file did not raise error!'
     976           
     977        #clean up
     978        os.remove(mesh_file)
     979        os.remove(point_file)
     980
     981    def test_fit_to_mesh_file_errorsII(self):
     982        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
     983        import tempfile
     984        import os
     985
     986        # create a .tsh file, no user outline
     987        mesh_file = tempfile.mktemp(".tsh")
     988        fd = open(mesh_file,'w')
     989        fd.write("unit testing a bad .tsh file \n")
     990        fd.close()
     991
     992        # create an .xya file
     993        point_file = tempfile.mktemp(".xya")
     994        fd = open(point_file,'w')
     995        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")
     996        fd.close()
     997
     998        mesh_output_file = "new_triangle.tsh"
     999        try:
     1000            fit_to_mesh_file(mesh_file, point_file,
     1001                             mesh_output_file, display_errors = False)
     1002        except IOError:
     1003            pass
     1004        else:
     1005            raise 'Bad file did not raise error!'
     1006           
     1007        #clean up
     1008        os.remove(mesh_file)
     1009        os.remove(point_file)
     1010
     1011    def test_fit_to_mesh_file_errorsIII(self):
     1012        from load_mesh.loadASCII import import_mesh_file, export_mesh_file
     1013        import tempfile
     1014        import os
     1015
     1016        # create a .tsh file, no user outline
     1017        mesh_dic = {}
     1018        mesh_dic['vertices'] = [[0.0, 0.0],[0.0, 5.0],[5.0, 0.0]]
     1019        mesh_dic['triangles'] =  [[0, 2, 1]]
     1020        mesh_dic['segments'] = [[0, 1], [2, 0], [1, 2]]
     1021        mesh_dic['triangle_tags'] = ['']
     1022        mesh_dic['vertex_attributes'] = [[1,2], [1,2], [1,2]]
     1023        mesh_dic['vertex_attribute_titles'] = ['density', 'temp']
     1024        mesh_dic['triangle_neighbors'] = [[-1, -1, -1]]
     1025        mesh_dic['segment_tags'] = ['external', 'external','external']
     1026        mesh_file = tempfile.mktemp(".tsh")
     1027        export_mesh_file(mesh_file,mesh_dic)
     1028
     1029        # create an .xya file
     1030        point_file = tempfile.mktemp(".xya")
     1031        fd = open(point_file,'w')
     1032        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")
     1033        fd.close()
     1034
     1035        #This a deliberately illegal filename to invoke the error.
     1036        mesh_output_file = ".../\z\z:ya.tsh"       
     1037
     1038        try:
     1039            fit_to_mesh_file(mesh_file, point_file,
     1040                             mesh_output_file, display_errors = False)
     1041        except IOError:
     1042            pass
     1043        else:
     1044            raise 'Bad file did not raise error!'
     1045       
     1046        #clean up
     1047        os.remove(mesh_file)
     1048        os.remove(point_file)
     1049 
    7301050    def Not_yet_test_smooth_att_to_mesh_with_excess_verts(self):
    7311051
  • anuga_core/source/anuga/geospatial_data/geospatial_data.py

    r4513 r4535  
    2424from anuga.utilities.anuga_exceptions import ANUGAError
    2525from anuga.config import points_file_block_line_size as MAX_READ_LINES   
     26
     27DEFAULT_ATTRIBUTE = 'elevation'
    2628
    2729class Geospatial_data:
     
    5355        of lists (or arrays) of length M. In the latter case the keys
    5456        in the dictionary represent the attribute names, in the former
    55         the attribute will get the default name "attribute".
     57        the attribute will get the default name "elevation".
    5658       
    5759        geo_reference: Object representing the origin of the data
     
    208210        if not isinstance(attributes, DictType):
    209211            #Convert single attribute into dictionary
    210             attributes = {'attribute': attributes}
     212            attributes = {DEFAULT_ATTRIBUTE: attributes}
    211213
    212214        #Check input attributes   
     
    514516        attributes = {}
    515517        if file_name[-4:]== ".xya":
    516             msg = 'Text file format is moving to comma seperated .txt files.'
    517             warn(msg, DeprecationWarning)
     518            # Maybe not phase-out, so we can load in geo-ref info
     519            #msg = 'Text file format is moving to comma seperated .txt files.'
     520            #warn(msg, DeprecationWarning)
    518521            try:
    519522                if delimiter == None:
  • anuga_core/source/anuga/geospatial_data/test_geospatial_data.py

    r4501 r4535  
    8484        points = [[1.0, 2.1], [3.0, 5.3]]
    8585        attributes = [2, 4]
    86         G = Geospatial_data(points, attributes)       
    87 
    88         assert G.attributes.keys()[0] == 'attribute'
     86        G = Geospatial_data(points, attributes)       
     87        assert G.attributes.keys()[0] == DEFAULT_ATTRIBUTE
    8988        assert allclose(G.attributes.values()[0], [2, 4])
    9089       
     
    119118        assert allclose(V, [2, 4])
    120119
    121         V = G.get_attributes('attribute') #Get by name
     120        V = G.get_attributes(DEFAULT_ATTRIBUTE) #Get by name
    122121        assert allclose(V, [2, 4])
    123122
  • anuga_core/source/anuga/shallow_water/data_manager.py

    r4526 r4535  
    16361636
    16371637    This function returns the names of the files produced.
     1638
     1639    It will also produce as many output files as there are input sww files.
    16381640    """
    16391641   
Note: See TracChangeset for help on using the changeset viewer.