source: inundation/ga/storm_surge/pmesh/load_mesh/loadASCII.py @ 668

Last change on this file since 668 was 666, checked in by duncan, 20 years ago

syntax fix

File size: 17.4 KB
Line 
1
2from string import  find, rfind
3
4def mesh_file_to_mesh_dictionary(fileName):
5    """Load a pmesh file.  Returning the mesh dictionary.
6    """
7    try:
8        meshdic = import_trianglulation(fileName)
9    except IOError, e:       
10        msg = 'Could not open file %s ' %fileName
11        raise IOError, msg
12    return meshdic
13
14def import_trianglulation(ofile):
15    """
16    reading triangulation and meshoutline info.
17    import a file, ofile, with the format
18   
19    First line:  <# of vertices> <# of attributes>
20    Following lines:  <vertex #> <x> <y> [attributes]
21    One line:  <# of triangles>
22    Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
23    One line:  <# of segments>
24    Following lines:  <segment #> <vertex #>  <vertex #> [boundary marker]
25    Note: This might throw a can't load file error
26    """
27    fd = open(ofile,'r')
28    dict = read_trianglulation(fd)
29    dict_mesh = read_mesh(fd)
30    for element in dict_mesh.keys():
31        dict[element] = dict_mesh[element]
32       
33    fd.close()
34    return dict
35
36
37def read_trianglulation(fd):
38    delimiter = " "
39    ######### loading the point info
40    line = fd.readline()
41    #print line
42    fragments = line.split()
43    #for fragment in fragments:
44    #    print fragment
45    if fragments ==[]:
46        NumOfVertices = 0
47        NumOfVertAttributes = 0
48    else:
49        NumOfVertices = fragments[0]
50        NumOfVertAttributes = fragments[1]
51    points = []
52    pointattributes = []
53    for index in range(int(NumOfVertices)):       
54        #print index
55        fragments = fd.readline().split()
56        #print fragments
57        fragments.pop(0) #pop off the index
58       
59        # pop the x y off so we're left with a list of attributes       
60        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
61        points.append(vert)
62        apointattributes  = []
63        #print fragments
64        for fragment in fragments:
65            apointattributes.append(float(fragment))
66        pointattributes.append(apointattributes)
67       
68    ######### loading the point title info
69    line = fd.readline()
70    #print "point title comments",line
71    vertTitle = []
72    for index in range(int(NumOfVertAttributes)):       
73        #print index
74        fragments = fd.readline().strip() 
75        vertTitle.append(fragments)
76        #print vertTitle
77       
78    ######### loading the triangle info
79    line = fd.readline()
80    #print "triangle comments",line
81    fragments = line.split()
82    #for fragment in fragments:
83    #    print fragment
84    NumOfTriangles = fragments[0]
85    triangles = []
86    triangleattributes = []
87    triangleneighbors = []
88    for index in range(int(NumOfTriangles)):       
89        #print index
90        line = fd.readline()
91        line.strip() # so we can get the region string
92        fragments = line.split()
93        #print "triangle info", fragments
94        fragments.pop(0) #pop off the index
95       
96        tri = [int(fragments[0]),int(fragments[1]),int(fragments[2])]
97        triangles.append(tri)
98        neighbors = [int(fragments[3]),int(fragments[4]),int(fragments[5])]
99        triangleneighbors.append(neighbors)
100        for x in range(7): # remove index [<vertex #>] [<neigbouring tri #>]
101            line = line[find(line,delimiter):] # remove index
102            line = line.lstrip()
103        stringmarker = line.strip()
104        triangleattributes.append([stringmarker])
105       
106    ######### loading the segment info
107    line = fd.readline()
108    #print "seg comment line",line
109    fragments = line.split()
110    #for fragment in fragments:
111    #    print fragment
112    NumOfSegments = fragments[0]
113    segments = []
114    segmentmarkers = []
115    for index in range(int(NumOfSegments)):       
116        #print index
117        line = fd.readline()
118        line.strip() # to get the segment string
119        fragments = line.split()
120        #print fragments
121        fragments.pop(0) #pop off the index
122        seg = [int(fragments[0]),int(fragments[1])]
123        segments.append(seg)
124        line = line[find(line,delimiter):] # remove index
125        line = line.lstrip()
126        line = line[find(line,delimiter):] # remove x
127        line = line.lstrip()
128        line = line[find(line,delimiter):] # remove y
129        stringmarker = line.strip()
130        segmentmarkers.append(stringmarker)
131    meshDict = {}
132    meshDict['generatedpointlist'] = points
133    meshDict['generatedpointattributelist'] = pointattributes
134    meshDict['generatedtrianglelist'] = triangles
135    meshDict['generatedtriangleattributelist'] = triangleattributes
136    meshDict['generatedtriangleneighborlist'] = triangleneighbors
137    meshDict['generatedsegmentlist'] = segments
138    meshDict['generatedsegmentmarkerlist'] = segmentmarkers
139    meshDict['generatedpointattributetitlelist'] = vertTitle
140    return meshDict
141   
142def import_mesh(ofile):
143    """
144    import a file, ofile, with the format
145   
146    First line:  <# of vertices> <# of attributes>
147    Following lines: <x> <y> [attributes]
148    One line:  <# of segments>
149    Following lines:  <vertex #>  <vertex #> [boundary marker]
150    Note: This might throw a can't load file error
151    """
152    fd = open(ofile,'r')
153    meshDict = read_mesh(fd)
154    fd.close()
155    return meshDict
156
157def read_mesh(fd):
158    """
159    Note, if a file has no mesh info, it can still be read - the meshdic
160    returned will be 'empty'.
161    """
162    delimiter = " " # warning: split() calls are using default whitespace
163   
164    ######### loading the point info
165    line = fd.readline()
166    #print line
167    fragments = line.split()
168    #for fragment in fragments:
169    if fragments ==[]:
170        NumOfVertices = 0
171        NumOfVertAttributes = 0
172    else:
173        NumOfVertices = fragments[0]
174        NumOfVertAttributes = fragments[1]
175    points = []
176    pointattributes = []
177    for index in range(int(NumOfVertices)):       
178        #print index
179        fragments = fd.readline().split() 
180        #print fragments
181        fragments.pop(0) #pop off the index
182        # pop the x y off so we're left with a list of attributes
183        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
184        points.append(vert)
185        apointattributes  = []
186        #print fragments
187        for fragment in fragments:
188            apointattributes.append(float(fragment))
189        pointattributes.append(apointattributes)
190       
191
192    ######### loading the segment info
193    line = fd.readline()
194    #print line
195    fragments = line.split()
196    #for fragment in fragments:
197    #    print fragment
198    if fragments ==[]:
199        NumOfSegments = 0
200    else:
201        NumOfSegments = fragments[0]
202    segments = []
203    segmentmarkers = []
204    for index in range(int(NumOfSegments)): 
205        #print index
206        line = fd.readline()
207        fragments = line.split()
208        #print fragments
209        fragments.pop(0) #pop off the index
210       
211        seg = [int(fragments[0]),int(fragments[1])]
212        segments.append(seg)
213        line = line[find(line,delimiter):] # remove index
214        line = line.lstrip()
215        line = line[find(line,delimiter):] # remove x
216        line = line.lstrip()
217        line = line[find(line,delimiter):] # remove y
218        stringmarker = line.strip()
219        segmentmarkers.append(stringmarker) 
220
221    ######### loading the hole info
222    line = fd.readline()
223    #print line
224    fragments = line.split()
225    #for fragment in fragments:
226    #    print fragment
227    if fragments ==[]:
228        numOfHoles = 0
229    else:
230        numOfHoles = fragments[0]
231    holes = []
232    for index in range(int(numOfHoles)):       
233        #print index
234        fragments = fd.readline().split()
235        #print fragments
236        fragments.pop(0) #pop off the index
237        hole = [float(fragments[0]),float(fragments[1])]
238        holes.append(hole)
239
240   
241    ######### loading the region info
242    line = fd.readline()
243    #print line
244    fragments = line.split()
245    #for fragment in fragments:
246    #    print fragment
247    if fragments ==[]:
248        numOfRegions = 0
249    else:
250        numOfRegions = fragments[0]
251    regions = []
252    regionattributes = []
253    for index in range(int(numOfRegions)):
254        line = fd.readline()
255        fragments = line.split()
256        #print fragments
257        fragments.pop(0) #pop off the index
258        region = [float(fragments[0]),float(fragments[1])]
259        regions.append(region)
260
261        line = line[find(line,delimiter):] # remove index
262        line = line.lstrip()
263        line = line[find(line,delimiter):] # remove x
264        line = line.lstrip()
265        line = line[find(line,delimiter):] # remove y
266        stringmarker = line.strip()
267        regionattributes.append(stringmarker)
268    regionmaxareas = []
269    line = fd.readline()
270    for index in range(int(numOfRegions)): # Read in the Max area info
271        line = fd.readline()
272        fragments = line.split()
273        #print fragments
274        # The try is here for format compatibility
275        try:
276            fragments.pop(0) #pop off the index
277            if len(fragments) == 0: #no max area
278                regionmaxareas.append(None)
279            else:
280                regionmaxareas.append(float(fragments[0]))
281        except IndexError, e:
282            regionmaxareas.append(None)
283       
284    meshDict = {}
285    meshDict['pointlist'] = points
286    meshDict['pointattributelist'] = pointattributes
287    meshDict['segmentlist'] = segments
288    meshDict['segmentmarkerlist'] = segmentmarkers
289    meshDict['holelist'] = holes
290    meshDict['regionlist'] = regions
291    meshDict['regionattributelist'] = regionattributes
292    meshDict['regionmaxarealist'] = regionmaxareas
293
294    return meshDict
295
296def clean_line(line,delimiter):     
297    """Remove whitespace
298    """
299    #print ">%s" %line
300    line = line.strip()
301    #print "stripped>%s" %line
302    numbers = line.split(delimiter)
303    i = len(numbers) - 1
304    while i >= 0:
305        if numbers[i] == '':
306            numbers.pop(i)
307        i += -1
308    #for num in numbers:
309    #    print "num>%s<" %num
310    return numbers
311
312def write_ASCII_trianglulation(fd,
313                               gen_dict):
314    vertices = gen_dict['generatedpointlist']
315    vertices_attributes = gen_dict['generatedpointattributelist']
316    try:
317        vertices_attribute_titles = gen_dict['generatedpointattributetitlelist']
318   
319    except KeyError, e:
320        #FIXME is this the best way?
321        if vertices_attributes == [] or vertices_attributes[0] == []:
322             vertices_attribute_titles = []
323        else:
324            raise KeyError, e
325       
326    triangles = gen_dict['generatedtrianglelist']
327    triangles_attributes = gen_dict['generatedtriangleattributelist']
328    triangle_neighbors = gen_dict['generatedtriangleneighborlist']
329       
330    segments = gen_dict['generatedsegmentlist']
331    segment_markers = gen_dict['generatedsegmentmarkerlist']
332     
333    numVert = str(len(vertices))
334    if (numVert == "0"):
335        numVertAttrib = "0"
336    else:
337        numVertAttrib = str(len(vertices_attributes[0]))
338    fd.write(numVert + " " + numVertAttrib + " # <vertex #> <x> <y> [attributes] ...Triangulation Vertices..." + "\n")
339
340    #<vertex #> <x> <y> [attributes]   
341    index = 0 
342    for vert in vertices:
343        attlist = ""
344        for att in vertices_attributes[index]:
345            attlist = attlist + str(att)+" "
346        attlist.strip()
347        fd.write(str(index) + " "
348                 + str(vert[0]) + " "
349                 + str(vert[1]) + " "
350                 + attlist + "\n")
351        index += 1
352
353    # write comments for title
354    fd.write("# attribute column titles ...Triangulation Vertex Titles..." + "\n")
355    for title in vertices_attribute_titles:
356        fd.write(title + "\n")
357       
358    #<# of triangles>
359    n = len(triangles)
360    fd.write(str(n) + " # <triangle #> [<vertex #>] [<neigbouring triangle #>] [attribute of region] ...Triangulation Triangles..." + "\n")
361       
362    # <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
363    for index in range(n):
364        neighbors = ""
365        tri = triangles[index]
366        for neighbor in triangle_neighbors[index]:
367            if neighbor:
368                neighbors += str(neighbor) + " "
369            else:
370                if neighbor == 0:
371                    neighbors +=  "0 "
372                else:
373                    neighbors +=  "-1 "
374        if triangles_attributes[index] == ['']:
375            att = ""
376        else:
377            att = str(triangles_attributes[index])
378        fd.write(str(index) + " "
379                 + str(tri[0]) + " " 
380                 + str(tri[1]) + " " 
381                 + str(tri[2]) + " " 
382                 + neighbors + " "
383                 + att + "\n")
384           
385    #One line:  <# of segments>
386    fd.write(str(len(segments)) + 
387             " # <segment #> <vertex #>  <vertex #> [boundary marker] ...Triangulation Segments..." + "\n")
388       
389    #Following lines:  <segment #> <vertex #>  <vertex #> [boundary marker]
390    for i in range(len(segments)):
391        seg = segments[i]
392        fd.write(str(i) + " "
393                 + str(seg[0]) + " " 
394                 + str(seg[1]) + " " 
395                 + str(segment_markers[i]) + "\n")
396
397
398def export_trianglulation_file(ofile,gen_dict):
399    """
400    write a file, ofile, with the format
401   
402    First line:  <# of vertices> <# of attributes>
403    Following lines:  <vertex #> <x> <y> [attributes]
404    One line:  <# of triangles>
405    Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
406    One line:  <# of segments>
407    Following lines:  <segment #> <vertex #>  <vertex #> [boundary marker]
408    """
409    try:
410        fd = open(ofile,'w')
411        write_ASCII_trianglulation(fd,gen_dict)
412        fd.close()
413    except IOError, e:       
414        msg = 'Could not write file %s ' %fileName
415        raise IOError, msg
416
417## used by alpha shapes
418   
419def export_boundary_file( file_name, points, title, delimiter = ','):
420    """
421    export a file, ofile, with the format
422   
423    First line: Title variable
424    Following lines:  [point index][delimiter][point index]
425
426    file_name - the name of the new file
427    points - List of point index pairs [[p1, p2],[p3, p4]..]
428    title - info to write in the first line
429    """
430   
431    fd = open(file_name,'w')
432   
433    fd.write(title+"\n")
434    #[point index][delimiter][point index]
435    for point in points:
436        fd.write( str(point[0]) + delimiter
437                  + str(point[1]) + "\n")
438    fd.close()
439
440
441###
442#  LOADING XYA FILES
443###
444 
445           
446def load_xya_file(ofile,delimiter = ','):
447    """
448    load a file, ofile, with the format
449    x,y, [attributes]
450    """
451    try:
452        fd = open(ofile)
453        #print "ofile",ofile
454        xya_dic = read_xya_file(fd, delimiter)
455        fd.close()
456    except IOError, e:       
457        msg = 'Could not open file %s ' %fileName
458        raise IOError, msg
459    return xya_dic
460
461def read_xya_file(fd, delimiter):
462    lines = fd.readlines()
463    points = []
464    pointattributes = []
465    if len(lines) <= 1:
466        raise SyntaxError
467    title = lines.pop(0) # the first (title) line
468    attLength = len(clean_line(lines[0],delimiter))-2
469   
470    #print "initlegth"
471    #print attLength
472    for line in lines:
473        #print "line >%s" %line
474        numbers = clean_line(line,delimiter)
475        #print "numbers >%s<" %numbers
476        if len(numbers) < 2 and numbers != []:
477            raise SyntaxError
478        if numbers != []:
479            try:
480                x = float(numbers[0])
481                y = float(numbers[1])
482                points.append([x,y])
483                numbers.pop(0)
484                numbers.pop(0)
485                attributes = []
486                if attLength != len(numbers):
487                    raise SyntaxError
488                   
489                attLength = len(numbers)
490
491                for num in numbers:
492                    num.strip()
493                    if num != '\n' and num != '':
494                        attributes.append(float(num))
495            except ValueError:
496                raise SyntaxError
497            pointattributes.append(attributes)
498    xya_dict = {}
499    xya_dict['pointlist'] = points
500    xya_dict['pointattributelist'] = pointattributes
501    xya_dict['title'] = title
502    xya_dict['segmentlist'] = []
503    xya_dict['segmentmarkerlist'] = []
504    xya_dict['regionlist'] = []
505    xya_dict['regionattributelist'] = []
506    xya_dict['regionmaxarealist'] = []
507    xya_dict['holelist'] = []
508   
509    return xya_dict
510
511               
512def export_xya_file( file_name, xya_dict, title, delimiter = ','):
513    """
514    export a file, ofile, with the format
515   
516    First line: Title variable
517    Following lines:  <vertex #> <x> <y> [attributes]
518
519    file_name - the name of the new file
520    xya_dict - point and point attribute info in a dictionary
521    title - info to write in the first line
522    """
523    #FIXME, move the test for this from meshharness to loadasciiharness
524    points = xya_dict['pointlist'] 
525    pointattributes = xya_dict['pointattributelist']
526   
527    fd = open(file_name,'w')
528   
529    fd.write(title+"\n")
530    #<vertex #> <x> <y> [attributes]
531    for vert, vertatts in map(None, points, pointattributes):
532        attlist = ""
533        for att in vertatts:
534            attlist = attlist + str(att)+ delimiter
535        attlist = attlist[0:-len(delimiter)] # remove the last delimiter
536        attlist.strip()
537        fd.write( str(vert[0]) + delimiter
538                  + str(vert[1]) + delimiter
539                  + attlist + "\n")
540    fd.close()
541     
542if __name__ == "__main__":
543    m = import_mesh("tee.txt")
544    print m
Note: See TracBrowser for help on using the repository browser.