Changeset 985

Timestamp:

Mar 1, 2005, 6:51:23 PM (20 years ago)

Author:

duncan

Message:

first cut of a netcdf .msh file

Location:

inundation/ga/storm_surge/pmesh/load_mesh

Files:

• loadASCII.py (modified) (7 diffs)
• test_loadASCII.py (modified) (5 diffs)

inundation/ga/storm_surge/pmesh/load_mesh/loadASCII.py

@@ -32 +32 @@
     segment_tags : [tag,tag,...] list of strings
     triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
-    triangle tags: [s1,s2,...] A list of strings
+    triangle tags: [[s1],[s2],...] A list of list of strings (probably not neccecary.  a list of string should be ok)
     triangle neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
         
@@ -53 +53 @@
 from string import  find, rfind
 
-from Numeric import array, Float, reshape, concatenate #, transpose
+from Numeric import array, Float, Int16, Character,reshape, concatenate
 
 from os.path import splitext
@@ -85 +85 @@
     Note: This might throw a can't load file error
     """
-    fd = open(ofile,'r')
-    dict = read_triangulation(fd)
-    dict_mesh = read_outline(fd)
-    for element in dict_mesh.keys():
-        dict[element] = dict_mesh[element]
-        
-    fd.close()
+    
+    if ofile[-4:]== ".tsh":
+        fd = open(ofile,'r')
+        dict = read_triangulation(fd)
+        dict_mesh = read_outline(fd)
+        for element in dict_mesh.keys():
+            dict[element] = dict_mesh[element]
+        fd.close()
+    elif ofile[-4:]== ".msh":
+        dict = read_tsh_file(ofile)
+    else:      
+        msg = 'Extension %s is unknown' %ofile[-4:]
+        raise IOError, msg        
     return dict
 
@@ -570 +576 @@
 
 
-def write_msh_file(filename, mesh_dict):
-    """Write .pts NetCDF file   
-    Return a Points dictionary.  see header for definition.
-
-    WARNING: This function mangles the mesh_dict data structure
+def write_msh_file(filename, mesh):
+    """Write .msh NetCDF file   
+
+    WARNING: This function mangles the mesh data structure
     """
   
-    from Scientific.IO.NetCDF import NetCDFFile  
-    
-    point_atts2array(point_atts)
+    from Scientific.IO.NetCDF import NetCDFFile
+
+    IntType = Int16
+    
+    # 
+    #the triangulation
+    mesh['vertices'] = array(mesh['vertices']).astype(Float)
+    mesh['vertex_attributes'] = array(mesh['vertex_attributes']).astype(Float)
+    mesh['vertex_attribute_titles'] = array(mesh['vertex_attribute_titles']).astype(Character) 
+    mesh['segments'] = array(mesh['segments']).astype(IntType)
+    mesh['segment_tags'] = array(mesh['segment_tags']).astype(Character)
+    mesh['triangles'] = array(mesh['triangles']).astype(IntType)  
+    mesh['triangle_tags'] = array(mesh['triangle_tags']).astype(Character)
+    mesh['triangle_neighbors'] = array(mesh['triangle_neighbors']).astype(IntType) 
+    
+    #the outline 
+    mesh['points'] = array(mesh['points']).astype(Float)
+    mesh['point_attributes'] = array(mesh['point_attributes']).astype(Float)
+    mesh['outline_segments'] = array(mesh['outline_segments']).astype(IntType)
+    mesh['outline_segment_tags'] = array(mesh['outline_segment_tags']).astype(Character)
+    mesh['holes'] = array(mesh['holes']).astype(Float)
+    mesh['regions'] = array(mesh['regions']).astype(Float)
+    mesh['region_tags'] = array(mesh['region_tags']).astype(Character)
+    mesh['region_max_areas'] = array(mesh['region_max_areas']).astype(Float)
+    
+    #mesh = mesh_dict2array(mesh)
+    #print "new_mesh",new_mesh
+    #print "mesh",mesh 
+    
     # NetCDF file definition
     outfile = NetCDFFile(filename, 'w')
@@ -588 +619 @@
                       'of spatial point data'
     
-    # dimension definitions
-    shape = point_atts['pointlist'].shape[0]
-    outfile.createDimension('number_of_points', shape)  
-    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
-    
-    # variable definition
-    outfile.createVariable('points', Float, ('number_of_points',
-                                             'number_of_dimensions'))
-
-    #create variables  
-    outfile.variables['points'][:] = point_atts['pointlist'] #.astype(Float32)
-    for key in point_atts['attributelist'].keys():
-        outfile.createVariable(key, Float, ('number_of_points',))
-        outfile.variables[key][:] = point_atts['attributelist'][key] #.astype(Float32)
-    outfile.close() 
-
+    # dimension definitions - fixed
+    outfile.createDimension('num_of_dimensions', 2) #This is 2d data
+    outfile.createDimension('num_of_segment_ends', 2) #Segs have two points
+    outfile.createDimension('num_of_triangle_vertices', 3)
+    outfile.createDimension('num_of_triangle_faces', 3)
+    outfile.createDimension('num_of_region_max_area', 1)
+    
+    # dimension definitions - variable
+    #trianglulation
+    outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
+    outfile.createDimension('num_of_vertex_attributes',
+                            mesh['vertex_attributes'].shape[1])
+    outfile.createDimension('num_of_vertex_attribute_title_chars',
+                            mesh['vertex_attribute_titles'].shape[1])
+    
+    outfile.createDimension('num_of_segments',
+                            mesh['segments'].shape[0])
+    outfile.createDimension('num_of_segment_tag_chars',
+                            mesh['segment_tags'].shape[1])
+    outfile.createDimension('num_of_triangles',
+                            mesh['triangles'].shape[0])
+    outfile.createDimension('num_of_triangle_tag_chars',
+                            mesh['triangle_tags'].shape[2]) 
+    outfile.createDimension('num_of_triangle_tag_per_triangle',
+                            mesh['triangle_tags'].shape[1])
+
+    #outline
+    outfile.createDimension('num_of_points', mesh['points'].shape[0])
+    outfile.createDimension('num_of_point_attributes',
+                            mesh['point_attributes'].shape[1])
+    outfile.createDimension('num_of_outline_segments',
+                            mesh['outline_segments'].shape[0])
+    outfile.createDimension('num_of_outline_segment_tag_chars',
+                            mesh['outline_segment_tags'].shape[1])
+    outfile.createDimension('num_of_holes', mesh['holes'].shape[0])
+    outfile.createDimension('num_of_regions', mesh['regions'].shape[0])
+    outfile.createDimension('num_of_region_tag_chars',
+                            mesh['region_tags'].shape[1])
+    
+    
+    # variable definitions
+    # trianglulation
+    outfile.createVariable('vertices', Float, ('num_of_vertices',
+                                             'num_of_dimensions'))
+    outfile.createVariable('vertex_attributes', Float, ('num_of_vertices',
+                                             'num_of_vertex_attributes'))
+    outfile.createVariable('vertex_attribute_titles',
+                           Character,
+                           ( 'num_of_vertex_attributes','num_of_vertex_attribute_title_chars' ))
+    outfile.createVariable('segments',
+                           IntType,
+                           ('num_of_segments', 'num_of_segment_ends'))
+    outfile.createVariable('segment_tags',
+                           Character,
+                           ('num_of_segments',
+                            'num_of_segment_tag_chars'))
+    outfile.createVariable('triangles',
+                           IntType,
+                           ('num_of_triangles', 'num_of_triangle_vertices'))
+    outfile.createVariable('triangle_tags',
+                           Character,
+                           ('num_of_triangles',
+                            'num_of_triangle_tag_per_triangle',
+                            'num_of_triangle_tag_chars'))
+    outfile.createVariable('triangle_neighbors',
+                           IntType,
+                           ('num_of_triangles', 'num_of_triangle_faces'))
+    
+    
+    # outline
+    outfile.createVariable('points', Float, ('num_of_points',
+                                             'num_of_dimensions'))
+    outfile.createVariable('point_attributes', Float, ('num_of_points',
+                                             'num_of_point_attributes'))
+    
+    outfile.createVariable('outline_segments',
+                           IntType,
+                           ('num_of_outline_segments', 'num_of_segment_ends'))
+    outfile.createVariable('outline_segment_tags',
+                           Character,
+                           ('num_of_outline_segments',
+                            'num_of_outline_segment_tag_chars'))
+    outfile.createVariable('holes', Float, ('num_of_holes',
+                                             'num_of_dimensions'))
+    outfile.createVariable('regions', Float, ('num_of_regions',
+                                             'num_of_dimensions'))
+    outfile.createVariable('region_tags',
+                           Character,
+                           ('num_of_regions',
+                            'num_of_region_tag_chars'))
+    outfile.createVariable('region_max_areas',
+                           Float,
+                           ('num_of_regions',)) #'num_of_region_max_area'))
+
+    # Setting the variables
+    #triangulation
+    outfile.variables['vertices'][:] = mesh['vertices']
+    outfile.variables['vertex_attributes'][:] = mesh['vertex_attributes']
+    outfile.variables['vertex_attribute_titles'][:] = mesh['vertex_attribute_titles']
+    outfile.variables['segments'][:] = mesh['segments']
+    outfile.variables['segment_tags'][:] = mesh['segment_tags']
+    outfile.variables['triangles'][:] = mesh['triangles']
+    outfile.variables['triangle_tags'][:] = mesh['triangle_tags']
+    outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
+    
+    # Outline
+    outfile.variables['points'][:] = mesh['points'] 
+    outfile.variables['point_attributes'][:] = mesh['point_attributes']
+    outfile.variables['outline_segments'][:] = mesh['outline_segments']
+    outfile.variables['outline_segment_tags'][:] = mesh['outline_segment_tags']
+    outfile.variables['holes'][:] = mesh['holes']
+    outfile.variables['regions'][:] = mesh['regions']
+    outfile.variables['region_tags'][:] = mesh['region_tags']
+    outfile.variables['region_max_areas'][:] = mesh['region_max_areas']
+    
+    outfile.close()
+
+
+    
+def read_msh_file(file_name):
+    """
+    Read in an msh file.
+
+    """
+        
+    from Scientific.IO.NetCDF import NetCDFFile         
+            
+    #Check contents
+    #Get NetCDF
+    fid = NetCDFFile(file_name, 'r') 
+
+    mesh = {}
+    # Get the variables
+    # the triangulation
+    mesh['vertices'] = fid.variables['vertices'][:]
+    mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:]
+    titles = fid.variables['vertex_attribute_titles'][:]
+    mesh['vertex_attribute_titles'] = []
+    for i, title in enumerate(titles):
+        mesh['vertex_attribute_titles'].append(titles[i].tostring().strip())
+    mesh['segments'] = fid.variables['segments'][:]
+    tags = fid.variables['segment_tags'][:]
+    mesh['segment_tags'] =[]
+    for i, tag in enumerate(tags):
+        mesh['segment_tags'].append(tags[i].tostring().strip())
+    mesh['triangles'] = fid.variables['triangles'][:]
+    tags = fid.variables['triangle_tags'][:]
+    mesh['triangle_tags'] =[]
+    for i, tag in enumerate(tags):
+        mesh['triangle_tags'].append(tags[i].tostring().strip())
+    mesh['triangle_neighbors'] = fid.variables['triangle_neighbors'][:]
+    
+    #the outline 
+    mesh['points'] = fid.variables['points'][:]
+    mesh['point_attributes'] = fid.variables['point_attributes'][:]
+    mesh['outline_segments'] = fid.variables['outline_segments'][:]
+    tags = fid.variables['outline_segment_tags'][:]
+    mesh['outline_segment_tags'] =[]
+    for i, tag in enumerate(tags):
+        mesh['outline_segment_tags'].append(tags[i].tostring().strip())
+    
+    mesh['holes'] = fid.variables['holes'][:]
+    mesh['regions'] = fid.variables['regions'][:]
+    tags = fid.variables['region_tags'][:]
+    mesh['region_tags'] =[]
+    for i, tag in enumerate(tags):
+        mesh['region_tags'].append(tags[i].tostring().strip())
+    mesh['region_max_areas'] = fid.variables['region_max_areas'][:]
+    #mesh[''] = fid.variables[''][:]
+    
+    fid.close()
+    
+    return mesh
+           
 
 ## used by alpha shapes
@@ -779 +968 @@
         outfile.variables[key][:] = point_atts['attributelist'][key] #.astype(Float32)
     outfile.close() 
-        
-# used to be load_xya_file            
+                  
 def load_points_file(ofile,delimiter = ','):
     """
-    load a file, ofile, with the format 1st line header,
-    x,y, [attributes]
+    load an .xya or .pts file   
     """
     
@@ -805 +992 @@
         msg = 'Extension %s is unknown' %ofile[-4:]
         raise IOError, msg
-        
-    
-        
     return points_dict
 

inundation/ga/storm_surge/pmesh/load_mesh/test_loadASCII.py

@@ -26 +26 @@
 class loadASCIITestCase(unittest.TestCase):
     def setUp(self):
-        pass
-    
+        self.dict ={}
+        self.dict['outline_segments'] = [(0, 1), (1, 2), (0, 2), (0, 3)]
+        self.dict['outline_segment_tags'] = ['50', '40', '30', '20']
+        self.dict['holes'] = [(0.2, 0.6)]
+        self.dict['point_attributes'] = [[5, 2], [4, 2], [3, 2], [2,2]]
+        self.dict['regions'] = [(0.3, 0.3),(0.3, 0.4)]
+        self.dict['region_tags'] = ['1.3', 'yeah']
+        self.dict['region_max_areas'] = [36.0,-7.1]
+        self.dict['points'] = [(0.0, 0.0), (0.0, 4.0), (4.0, 0.0), (1.0, 1.0)]
+
+        self.dict['vertices'] = [(0.0, 0.0), (0.0, 4.0),
+                                          (4.0, 0.0), (1.0, 1.0), (2.0, 2.0)]
+        self.dict['triangles'] = [(3, 2, 4), (1, 0, 3),
+                                             (3, 4,1), (2, 3, 0)]
+        self.dict['segments'] = [(0, 1), (1, 4), (2, 0),
+                                            (0, 3), (4, 2)]
+        self.dict['triangle_tags'] = [['1.3'], ['1.3'],
+                                                      ['1.3'], ['1.3']]
+        self.dict['vertex_attributes'] = [[1.2,2.], [1.2,2.], [1.2,2.], [1.2,2.], [1.2,3.]]
+        self.dict['triangle_neighbors'] = [[-1, 2, 3], [3, 2, -1],
+                                                     [-1, 1, 0], [1, -1, 0]]
+        self.dict['segment_tags'] = ['50', '40', '30', '20', '40']
+        self.dict['vertex_attribute_titles'] = ['bed elevation', 'height']
     def tearDown(self):
         pass
@@ -33 +54 @@
     
     def test_import_mesh(self):
-        dict ={}
-        dict['outline_segments'] = [(0, 1), (1, 2), (0, 2), (0, 3)]
-        dict['outline_segment_tags'] = ['50', '40', '30', '20']
-        dict['holes'] = [(0.2, 0.6)]
-        dict['point_attributes'] = [[5, 2], [4, 2], [3, 2], [2,2]]
-        dict['regions'] = [(0.3, 0.3)]
-        dict['region_tags'] = ['1.3']
-        dict['region_max_areas'] = [36.0]
-        dict['points'] = [(0.0, 0.0), (0.0, 4.0), (4.0, 0.0), (1.0, 1.0)]
-
+        
+        dict = self.dict
         fileName = tempfile.mktemp(".txt")
         fd = open(fileName,'w')
@@ -70 +83 @@
                         array(dict['outline_segment_tags']),
                          'test_import_mesh failed. Test 4')
-        #print "loaded_dict['regions']",loaded_dict['regions']
-        #print "dict['regions']",dict['regions'] 
         self.failUnless(array(loaded_dict['regions'])  ==
                         array(dict['regions']),
@@ -89 +100 @@
         import os
         import tempfile
-
-        meshDict = {} 
-        meshDict['vertices'] = [(0.0, 0.0), (0.0, 4.0),
-                                          (4.0, 0.0), (1.0, 1.0), (2.0, 2.0)]
-        meshDict['triangles'] = [(3, 2, 4), (1, 0, 3),
-                                             (3, 4,1), (2, 3, 0)]
-        meshDict['segments'] = [(0, 1), (1, 4), (2, 0),
-                                            (0, 3), (4, 2)]
-        meshDict['triangle_tags'] = [['1.3'], ['1.3'],
-                                                      ['1.3'], ['1.3']]
-        meshDict['vertex_attributes'] = [[1.2,2.], [1.2,2.], [1.2,2.], [1.2,2.], [1.2,3.]]
-        meshDict['triangle_neighbors'] = [[-1, 2, 3], [3, 2, -1],
-                                                     [-1, 1, 0], [1, -1, 0]]
-        meshDict['segment_tags'] = ['50', '40', '30', '20', '40']
-        meshDict['vertex_attribute_titles'] = ['bed elevation', 'height']
-        
+        
+        meshDict = self.dict
         fileName = tempfile.mktemp(".tsh")
         export_triangulation_file(fileName, meshDict)
@@ -142 +139 @@
             
         os.remove(fileName)
-
+ 
+    def test_read_write_msh_file(self):
+        
+        dict = self.dict.copy()
+        fileName = tempfile.mktemp(".txt")
+        write_msh_file(fileName,dict)
+        loaded_dict = read_msh_file(fileName)
+        os.remove(fileName)
+
+        
+        dict = self.dict
+        #print "*********************"
+        #print dict
+        #print "**loaded_dict*******************"
+        #print loaded_dict
+        #print "*********************"
+
+        
+        self.failUnless(array(loaded_dict['points'])  ==
+                        array(dict['points']),
+                         'test_import_mesh failed. Test 1')
+        self.failUnless(array(loaded_dict['point_attributes'])  ==
+                        array(dict['point_attributes']),
+                         'test_import_mesh failed. Test 2')
+        self.failUnless(array(loaded_dict['outline_segments'])  ==
+                        array(dict['outline_segments']),
+                         'test_import_mesh failed. Test 3')
+        self.failUnless(array(loaded_dict['outline_segment_tags'])  ==
+                        array(dict['outline_segment_tags']),
+                         'test_import_mesh failed. Test 4')
+        self.failUnless(array(loaded_dict['regions'])  ==
+                        array(dict['regions']),
+                         'test_import_mesh failed. Test 5')
+        self.failUnless(array(loaded_dict['region_tags'])  ==
+                        array(dict['region_tags']),
+                         'test_import_mesh failed. Test 5')
+        self.failUnless(array(loaded_dict['region_max_areas'])  ==
+                        array(dict['region_max_areas']),
+                         'test_import_mesh failed. Test 5')
+
+        self.failUnless(array(loaded_dict['holes'])  ==
+                        array(dict['holes']),
+                         'test_import_mesh failed. Test 6')
+   
+        self.failUnless(array(dict['vertices'])  ==
+                        array(loaded_dict['vertices']),
+                         'test_export_triangulation_file failed. Test 1')
+        self.failUnless(array(dict['triangles'])  ==
+                        array(loaded_dict['triangles']),
+                         'test_export_triangulation_file failed. Test 2')
+        self.failUnless(array(dict['segments'])  ==
+                        array(loaded_dict['segments']),
+                         'test_export_triangulation_file failed. Test 3')
+        self.failUnless(array(dict['triangle_tags'])  ==
+                        array(loaded_dict['triangle_tags']),
+                         'test_export_triangulation_file failed. Test 4')
+        
+        self.failUnless(dict['vertex_attributes']  ==
+                        loaded_dict['vertex_attributes'],
+                         'test_export_triangulation_file failed. Test 5')
+        self.failUnless(array(dict['triangle_neighbors'])  ==
+                        array(loaded_dict['triangle_neighbors']),
+                         'test_export_triangulation_file failed. Test 6')
+        self.failUnless(array(dict['segment_tags'])  ==
+                        array(loaded_dict['segment_tags']),
+                         'test_export_triangulation_file failed. Test 7')
+        self.failUnless(array(dict['vertex_attribute_titles'])  ==
+                        array(loaded_dict['vertex_attribute_titles']),
+                         'test_export_triangulation_file failed. Test 8')
+        
     def test_loadpts(self):