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

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

refactoring / changing method names
File size: 47.9 KB
Line 
1"""
2The format for a .xya file is:
31st line:     [attribute names]
4other lines:  x y [attributes]
5
6for example:
7elevation
80.6 0.7 4.9
91.9 2.8 5
102.7 2.4 5.2
11
12The first two columns are always implicitly assumed to be x, y coordinates.
13
14
15The format for a Points dictionary is:
16
17  ['pointlist'] a 2 column array describing points. 1st column x, 2nd column y.
18  ['attributelist'], a dictionary of 1D arrays, representing attribute values
19  at the point.  The dictionary key is the attribute header.
20    eg
21    dic['pointlist'] = [[1.0,2.0],[3.0,5.0]]
22    dic['attributelist']['elevation'] = [[7.0,5.0]
23
24   
25    The dict format for IO with mesh files is;
26    (the triangulation)
27    vertices: [[x1,y1],[x2,y2],...] (lists of doubles)
28    vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
29    vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings)
30    segments: [[v1,v2],[v3,v4],...] (lists of integers)
31    segment_tags : [tag,tag,...] list of strings
32    triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points
33    triangle_tags: [s1,s2,...] A list of list of strings (probably not neccecary.  a list of string should be ok)
34    triangle_neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles
35       
36    (the outline)   
37    points: [[x1,y1],[x2,y2],...] (lists of doubles)
38    point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles)
39    outline_segments: [[point1,point2],[p3,p4],...] (lists of integers)
40    outline_segment_tags : [tag1,tag2,...] list of strings
41    holes : [[x1,y1],...](List of doubles, one inside each hole region)
42    regions : [ [x1,y1],...] (List of 4 doubles)
43    region_tags : [tag1,tag2,...] (list of strings)
44    region_max_areas: [ma1,ma2,...] (A list of doubles)
45    {Convension: A -ve max area means no max area}
46
47    Points files are .xya for ascii and .pts for NetCDF
48    Mesh files are .tsh for ascii and .msh for NetCDF
49
50    Note: point att's have no titles - that's currently ok, since least
51    squares adds the point info to vertices, and they have titles
52
53    The format for a .tsh file is (FIXME update this)
54   
55    First line:  <# of vertices> <# of attributes>
56    Following lines:  <vertex #> <x> <y> [attributes]
57    One line:  <# of triangles>
58    Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
59    One line:  <# of segments>
60    Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
61"""
62##FIXME (DSG-DSG) Is the dict format mentioned above a list of a numeric array?
63#  Needs to be defined
64##FIXME (DSG-DSG) if the ascii file isn't .xya give a better error message.
65
66
67from string import  find, rfind
68from Numeric import array, Float, Int16, Int32, Character,reshape, concatenate, take
69from os.path import splitext
70
71from coordinate_transforms.geo_reference import Geo_reference
72
73import exceptions
74class TitleAmountError(exceptions.Exception): pass
75
76NOMAXAREA=-999
77
78#FIXME (DSG-DSG) this is used by least_squares -
79#                   is this the best function for it to use?
80# (DSG-DSG) put an exception handler in import_triangulation/import_mesh_file
81# Think about exceptions before doing this!
82def mesh_file_to_mesh_dictionary(fileName):
83    """Load a pmesh file.  Returning the mesh dictionary.
84    """
85    try:
86        meshdic = import_mesh_file(fileName)
87    except IOError, e:       
88        msg = 'Could not open file %s ' %fileName
89        raise IOError, msg
90    return meshdic
91
92# This is used by pmesh
93def import_mesh_file(ofile):
94    """
95    read a mesh file, either .tsh or .msh
96   
97    Note: This might throw a can't load file error
98    """
99   
100    if ofile[-4:]== ".tsh":
101        dict = _read_tsh_file(ofile)
102    elif ofile[-4:]== ".msh":
103        dict = _read_msh_file(ofile)
104    else:     
105        msg = 'Extension %s is unknown' %ofile[-4:]
106        raise IOError, msg       
107    return dict
108
109
110def export_mesh_file(ofile,mesh_dict):
111    """
112    write a file, ofile, with the format
113   
114    First line:  <# of vertices> <# of attributes>
115    Following lines:  <vertex #> <x> <y> [attributes]
116    One line:  <# of triangles>
117    Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
118    One line:  <# of segments>
119    Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
120    """
121    #FIXME DSG-DSG: automate
122    if not mesh_dict.has_key('points'):
123        mesh_dict['points'] = []
124    if not mesh_dict.has_key('point_attributes'):
125        mesh_dict['point_attributes'] = []
126    if not mesh_dict.has_key('outline_segments'):
127        mesh_dict['outline_segments'] = []
128    if not mesh_dict.has_key('outline_segment_tags'):
129        mesh_dict['outline_segment_tags'] = []
130    if not mesh_dict.has_key('holes'):
131        mesh_dict['holes'] = []
132    if not mesh_dict.has_key('regions'):
133        mesh_dict['regions'] = []
134       
135    if not mesh_dict.has_key('region_tags'):
136        mesh_dict['region_tags'] = []
137    if not mesh_dict.has_key('region_max_areas'):
138        mesh_dict['region_max_areas'] = []
139    if not mesh_dict.has_key('vertices'):
140        mesh_dict['vertices'] = []
141    if not mesh_dict.has_key('vertex_attributes'):
142        mesh_dict['vertex_attributes'] = []
143    if not mesh_dict.has_key('vertex_attribute_titles'):
144        mesh_dict['vertex_attribute_titles'] = []
145    if not mesh_dict.has_key('segments'):
146        mesh_dict['segments'] = []
147    if not mesh_dict.has_key('segment_tags'):
148        mesh_dict['segment_tags'] = []
149    if not mesh_dict.has_key('triangles'):
150        mesh_dict['triangles'] = []
151    if not mesh_dict.has_key('triangle_tags'):
152        mesh_dict['triangle_tags'] = []
153    if not mesh_dict.has_key('triangle_neighbors'):
154        mesh_dict['triangle_neighbors'] = []
155    #print "DSG************"
156    #print "mesh_dict",mesh_dict
157    #print "DSG************"
158   
159    try:
160        if (ofile[-4:] == ".tsh"):
161            _write_tsh_file(ofile,mesh_dict)
162        elif (ofile[-4:] == ".msh"):
163            _write_msh_file(ofile, mesh_dict)
164    except IOError, e:       
165        msg = 'Could not write file %s ' %fileName
166        raise IOError, msg
167
168
169def _read_tsh_file(ofile):
170    """
171    Read the text file format for meshes
172    """
173    fd = open(ofile,'r')
174    dict = _read_triangulation(fd)
175    dict_mesh = _read_outline(fd)
176    for element in dict_mesh.keys():
177        dict[element] = dict_mesh[element]
178    fd.close()
179    return dict
180
181   
182def _read_triangulation(fd):
183    """
184    Read the generated triangulation, NOT the outline.
185    """
186    delimiter = " "
187    ######### loading the point info
188    line = fd.readline()
189    #print line
190    fragments = line.split()
191    #for fragment in fragments:
192    #print fragment
193    if fragments ==[]:
194        NumOfVertices = 0
195        NumOfVertAttributes = 0
196    else:
197        NumOfVertices = fragments[0]
198        NumOfVertAttributes = fragments[1]
199    points = []
200    pointattributes = []
201    for index in range(int(NumOfVertices)):       
202        #print index
203        fragments = fd.readline().split()
204        #print fragments
205        fragments.pop(0) #pop off the index
206       
207        # pop the x y off so we're left with a list of attributes       
208        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
209        points.append(vert)
210        apointattributes  = []
211        #print fragments
212        for fragment in fragments:
213            apointattributes.append(float(fragment))
214        pointattributes.append(apointattributes)
215       
216    ######### loading the point title info
217    line = fd.readline()
218    #print "point title comments",line
219    vertTitle = []
220    for index in range(int(NumOfVertAttributes)):       
221        #print index
222        fragments = fd.readline().strip() 
223        vertTitle.append(fragments)
224        #print vertTitle
225       
226    ######### loading the triangle info
227    line = fd.readline()
228    #print "triangle comments",line
229    fragments = line.split()
230    #for fragment in fragments:
231    #    print fragment
232    NumOfTriangles = fragments[0]
233    triangles = []
234    triangleattributes = []
235    triangleneighbors = []
236    for index in range(int(NumOfTriangles)):       
237        #print index
238        line = fd.readline()
239        line.strip() # so we can get the region string
240        fragments = line.split()
241        #print "triangle info", fragments
242        fragments.pop(0) #pop off the index
243       
244        tri = [int(fragments[0]),int(fragments[1]),int(fragments[2])]
245        triangles.append(tri)
246        neighbors = [int(fragments[3]),int(fragments[4]),int(fragments[5])]
247        triangleneighbors.append(neighbors)
248        for x in range(7): # remove index [<vertex #>] [<neigbouring tri #>]
249            line = line[find(line,delimiter):] # remove index
250            line = line.lstrip()
251        stringtag = line.strip()
252        triangleattributes.append(stringtag)
253       
254    ######### loading the segment info
255    line = fd.readline()
256    #print "seg comment line",line
257    fragments = line.split()
258    #for fragment in fragments:
259    #    print fragment
260    NumOfSegments = fragments[0]
261    segments = []
262    segmenttags = []
263    for index in range(int(NumOfSegments)):       
264        #print index
265        line = fd.readline()
266        line.strip() # to get the segment string
267        fragments = line.split()
268        #print fragments
269        fragments.pop(0) #pop off the index
270        seg = [int(fragments[0]),int(fragments[1])]
271        segments.append(seg)
272        line = line[find(line,delimiter):] # remove index
273        line = line.lstrip()
274        line = line[find(line,delimiter):] # remove x
275        line = line.lstrip()
276        line = line[find(line,delimiter):] # remove y
277        stringtag = line.strip()
278        segmenttags.append(stringtag)
279
280           
281    meshDict = {}
282    meshDict['vertices'] = points
283    meshDict['vertex_attributes'] = pointattributes
284    meshDict['triangles'] = triangles
285    meshDict['triangle_tags'] = triangleattributes
286    meshDict['triangle_neighbors'] = triangleneighbors
287    meshDict['segments'] = segments
288    meshDict['segment_tags'] = segmenttags
289    meshDict['vertex_attribute_titles'] = vertTitle
290   
291    return meshDict
292   
293def _read_outline(fd):
294    """
295    Note, if a file has no mesh info, it can still be read - the meshdic
296    returned will be 'empty'.
297    """
298    delimiter = " " # warning: split() calls are using default whitespace
299   
300    ######### loading the point info
301    line = fd.readline()
302    #print 'read_outline - line',line
303    fragments = line.split()
304    #for fragment in fragments:
305    if fragments ==[]:
306        NumOfVertices = 0
307        NumOfVertAttributes = 0
308    else:
309        NumOfVertices = fragments[0]
310        NumOfVertAttributes = fragments[1]
311    points = []
312    pointattributes = []
313    for index in range(int(NumOfVertices)):       
314        #print index
315        fragments = fd.readline().split() 
316        #print fragments
317        fragments.pop(0) #pop off the index
318        # pop the x y off so we're left with a list of attributes
319        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
320        points.append(vert)
321        apointattributes  = []
322        #print fragments
323        for fragment in fragments:
324            apointattributes.append(float(fragment))
325        pointattributes.append(apointattributes)
326       
327
328    ######### loading the segment info
329    line = fd.readline()
330    fragments = line.split()
331    #for fragment in fragments:
332    #    print fragment
333    if fragments ==[]:
334        NumOfSegments = 0
335    else:
336        NumOfSegments = fragments[0]
337    segments = []
338    segmenttags = []
339    for index in range(int(NumOfSegments)): 
340        #print index
341        line = fd.readline()
342        fragments = line.split()
343        #print fragments
344        fragments.pop(0) #pop off the index
345       
346        seg = [int(fragments[0]),int(fragments[1])]
347        segments.append(seg)
348        line = line[find(line,delimiter):] # remove index
349        line = line.lstrip()
350        line = line[find(line,delimiter):] # remove x
351        line = line.lstrip()
352        line = line[find(line,delimiter):] # remove y
353        stringtag = line.strip()
354        segmenttags.append(stringtag) 
355
356    ######### loading the hole info
357    line = fd.readline()
358    #print line
359    fragments = line.split()
360    #for fragment in fragments:
361    #    print fragment
362    if fragments ==[]:
363        numOfHoles = 0
364    else:
365        numOfHoles = fragments[0]
366    holes = []
367    for index in range(int(numOfHoles)):       
368        #print index
369        fragments = fd.readline().split()
370        #print fragments
371        fragments.pop(0) #pop off the index
372        hole = [float(fragments[0]),float(fragments[1])]
373        holes.append(hole)
374
375   
376    ######### loading the region info
377    line = fd.readline()
378    #print line
379    fragments = line.split()
380    #for fragment in fragments:
381    #    print fragment
382    if fragments ==[]:
383        numOfRegions = 0
384    else:
385        numOfRegions = fragments[0]
386    regions = []
387    regionattributes = []
388    for index in range(int(numOfRegions)):
389        line = fd.readline()
390        fragments = line.split()
391        #print fragments
392        fragments.pop(0) #pop off the index
393        region = [float(fragments[0]),float(fragments[1])]
394        regions.append(region)
395
396        line = line[find(line,delimiter):] # remove index
397        line = line.lstrip()
398        line = line[find(line,delimiter):] # remove x
399        line = line.lstrip()
400        line = line[find(line,delimiter):] # remove y
401        stringtag = line.strip()
402        regionattributes.append(stringtag)
403    regionmaxareas = []
404    line = fd.readline()
405    for index in range(int(numOfRegions)): # Read in the Max area info
406        line = fd.readline()
407        fragments = line.split()
408        #print fragments
409        # The try is here for format compatibility
410        try:
411            fragments.pop(0) #pop off the index
412            if len(fragments) == 0: #no max area
413                regionmaxareas.append(None)
414            else:
415                regionmaxareas.append(float(fragments[0]))
416        except IndexError, e:
417            regionmaxareas.append(None)
418
419    try:
420        geo_reference = Geo_reference(ASCIIFile=fd)
421    except:
422        #geo_ref not compulsory
423        geo_reference = None
424       
425    meshDict = {}
426    meshDict['points'] = points
427    meshDict['point_attributes'] = pointattributes
428    meshDict['outline_segments'] = segments
429    meshDict['outline_segment_tags'] = segmenttags
430    meshDict['holes'] = holes
431    meshDict['regions'] = regions
432    meshDict['region_tags'] = regionattributes
433    meshDict['region_max_areas'] = regionmaxareas
434    meshDict['geo_reference'] = geo_reference
435   
436    #print "meshDict",meshDict
437    return meshDict
438
439def clean_line(line,delimiter):     
440    """Remove whitespace
441    """
442    #print ">%s" %line
443    line = line.strip()
444    #print "stripped>%s" %line
445    numbers = line.split(delimiter)
446    i = len(numbers) - 1
447    while i >= 0:
448        if numbers[i] == '':
449            numbers.pop(i)
450        else:
451            numbers[i] = numbers[i].strip()
452       
453        i += -1
454    #for num in numbers:
455    #    print "num>%s<" %num
456    return numbers
457
458def _write_ASCII_triangulation(fd,
459                               gen_dict):
460    vertices = gen_dict['vertices']
461    vertices_attributes = gen_dict['vertex_attributes']
462    try:
463        vertices_attribute_titles = gen_dict['vertex_attribute_titles']   
464    except KeyError, e:
465        #FIXME is this the best way?
466        if vertices_attributes == [] or vertices_attributes[0] == []:
467             vertices_attribute_titles = []
468        else:
469            raise KeyError, e
470       
471    triangles = gen_dict['triangles']
472    triangles_attributes = gen_dict['triangle_tags']
473    triangle_neighbors = gen_dict['triangle_neighbors']
474       
475    segments = gen_dict['segments']
476    segment_tags = gen_dict['segment_tags']
477     
478    numVert = str(len(vertices))
479    if (numVert == "0") or len(vertices_attributes) == 0:
480        numVertAttrib = "0"
481    else:
482        numVertAttrib = str(len(vertices_attributes[0]))
483    fd.write(numVert + " " + numVertAttrib + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Triangulation Vertices..." + "\n")
484
485    #<vertex #> <x> <y> [attributes]   
486    index = 0 
487    for vert in vertices:
488        attlist = ""
489        for att in vertices_attributes[index]:
490            attlist = attlist + str(att)+" "
491        attlist.strip()
492        fd.write(str(index) + " "
493                 + str(vert[0]) + " "
494                 + str(vert[1]) + " "
495                 + attlist + "\n")
496        index += 1
497
498    # write comments for title
499    fd.write("# attribute column titles ...Triangulation Vertex Titles..." + "\n")
500    for title in vertices_attribute_titles:
501        fd.write(title + "\n")
502       
503    #<# of triangles>
504    n = len(triangles)
505    fd.write(str(n) + " # <# of triangles>, next lines <triangle #> [<vertex #>] [<neigbouring triangle #>] [attribute of region] ...Triangulation Triangles..." + "\n")
506       
507    # <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
508    for index in range(n):
509        neighbors = ""
510        tri = triangles[index]
511        for neighbor in triangle_neighbors[index]:
512            if neighbor:
513                neighbors += str(neighbor) + " "
514            else:
515                if neighbor == 0:
516                    neighbors +=  "0 "
517                else:
518                    neighbors +=  "-1 "
519        #Warning even though a list is past, only the first value
520        #is written.  There's an assumption that the list only
521        # contains one item. This assumption is made since the
522        #dict that's being passed around is also be used to communicate
523        # with triangle, and it seems to have the option of returning
524        # more than one value for triangle attributex
525        if triangles_attributes[index] == ['']:
526            att = ""
527        else:
528            att = str(triangles_attributes[index])
529        fd.write(str(index) + " "
530                 + str(tri[0]) + " " 
531                 + str(tri[1]) + " " 
532                 + str(tri[2]) + " " 
533                 + neighbors + " "
534                 + att + "\n")
535           
536    #One line:  <# of segments>
537    fd.write(str(len(segments)) + 
538             " # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Triangulation Segments..." + "\n")
539       
540    #Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
541    for i in range(len(segments)):
542        seg = segments[i]
543        fd.write(str(i) + " "
544                 + str(seg[0]) + " " 
545                 + str(seg[1]) + " " 
546                 + str(segment_tags[i]) + "\n")
547
548
549def export_mesh_file(ofile,mesh_dict):
550    """
551    write a file, ofile, with the format
552   
553    First line:  <# of vertices> <# of attributes>
554    Following lines:  <vertex #> <x> <y> [attributes]
555    One line:  <# of triangles>
556    Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
557    One line:  <# of segments>
558    Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
559    """
560    if not mesh_dict.has_key('points'):
561        mesh_dict['points'] = []
562    if not mesh_dict.has_key('point_attributes'):
563        mesh_dict['point_attributes'] = []
564    if not mesh_dict.has_key('outline_segments'):
565        mesh_dict['outline_segments'] = []
566    if not mesh_dict.has_key('outline_segment_tags'):
567        mesh_dict['outline_segment_tags'] = []
568    if not mesh_dict.has_key('holes'):
569        mesh_dict['holes'] = []
570    if not mesh_dict.has_key('regions'):
571        mesh_dict['regions'] = []
572       
573    if not mesh_dict.has_key('region_tags'):
574        mesh_dict['region_tags'] = []
575    if not mesh_dict.has_key('region_max_areas'):
576        mesh_dict['region_max_areas'] = []
577    if not mesh_dict.has_key('vertices'):
578        mesh_dict['vertices'] = []
579    if not mesh_dict.has_key('vertex_attributes'):
580        mesh_dict['vertex_attributes'] = []
581    if not mesh_dict.has_key('vertex_attribute_titles'):
582        mesh_dict['vertex_attribute_titles'] = []
583    if not mesh_dict.has_key('segments'):
584        mesh_dict['segments'] = []
585    if not mesh_dict.has_key('segment_tags'):
586        mesh_dict['segment_tags'] = []
587    if not mesh_dict.has_key('triangles'):
588        mesh_dict['triangles'] = []
589    if not mesh_dict.has_key('triangle_tags'):
590        mesh_dict['triangle_tags'] = []
591    if not mesh_dict.has_key('triangle_neighbors'):
592        mesh_dict['triangle_neighbors'] = []
593    #print "DSG************"
594    #print "mesh_dict",mesh_dict
595    #print "DSG************"
596   
597    try:
598        if (ofile[-4:] == ".tsh"):
599            _write_tsh_file(ofile,mesh_dict)
600        elif (ofile[-4:] == ".msh"):
601            _write_msh_file(ofile, mesh_dict)
602    except IOError, e:       
603        msg = 'Could not write file %s ' %fileName
604        raise IOError, msg
605
606def _write_tsh_file(ofile,mesh_dict):
607            fd = open(ofile,'w')
608            _write_ASCII_triangulation(fd,mesh_dict)
609            _write_ASCII_outline(fd,mesh_dict)   
610            fd.close()
611
612def _write_ASCII_outline(fd,
613                     dict):
614    points = dict['points']
615    point_attributes = dict['point_attributes']
616    segments = dict['outline_segments']
617    segment_tags = dict['outline_segment_tags']
618    holes = dict['holes']
619    regions = dict['regions']
620    region_tags = dict['region_tags']
621    region_max_areas = dict['region_max_areas']
622       
623    num_points = str(len(points))
624    if (num_points == "0"):
625        num_point_atts = "0"
626    else:
627        num_point_atts = str(len(point_attributes[0]))
628    fd.write(num_points + " " + num_point_atts + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Mesh Vertices..." + "\n")
629       
630    # <x> <y> [attributes]
631    for i,point in enumerate(points):
632        attlist = ""
633        for att in point_attributes[i]:
634            attlist = attlist + str(att)+" "
635        attlist.strip()
636        fd.write(str(i) + " "
637                 +str(point[0]) + " "
638                 + str(point[1]) + " "
639                 + attlist + "\n")
640           
641    #One line:  <# of segments>
642    fd.write(str(len(segments)) + 
643             " # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Mesh Segments..." + "\n")
644       
645    #Following lines:  <vertex #>  <vertex #> [boundary tag]
646    for i,seg in enumerate(segments):
647        fd.write(str(i) + " "
648                 + str(seg[0]) + " " 
649                 + str(seg[1]) + " " 
650                 + str(segment_tags[i]) + "\n")
651
652    #One line:  <# of holes>
653    fd.write(str(len(holes)) + 
654             " # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes..." + "\n")   
655    # <x> <y>
656    for i,h in enumerate(holes):
657        fd.write(str(i) + " "
658                 + str(h[0]) + " "
659                 + str(h[1]) + "\n")
660       
661    #One line:  <# of regions>
662    fd.write(str(len(regions)) + 
663             " # <# of regions>, next lines <Region #> <x> <y> <tag>...Mesh Regions..." + "\n")   
664    # <index> <x> <y> <tag>
665    #print "regions",regions
666    for i,r in enumerate(regions):
667        #print "r",r
668        fd.write(str(i) + " " 
669                 + str(r[0]) + " "
670                 + str(r[1])+ " "
671                 + str(region_tags[i]) + "\n")
672       
673    # <index> [<MaxArea>|""]
674       
675    #One line:  <# of regions>
676    fd.write(str(len(regions)) + 
677             " # <# of regions>, next lines <Region #> [Max Area] ...Mesh Regions..." + "\n") 
678    for i,r in enumerate(regions):
679        if region_max_areas[i] <= 0 :
680            area = ""
681        else:
682            area = str(region_max_areas[i])
683           
684        fd.write(str(i) + " " + area + "\n")
685
686    # geo_reference info
687    if dict.has_key('geo_reference') and not dict['geo_reference'] is None:
688        dict['geo_reference'].write_ASCII(fd)
689
690def _write_msh_file(filename, mesh):
691    """Write .msh NetCDF file   
692
693    WARNING: This function mangles the mesh data structure
694    """
695 
696    from Scientific.IO.NetCDF import NetCDFFile
697
698    IntType = Int32
699   
700    #
701    #the triangulation
702    mesh['vertices'] = array(mesh['vertices']).astype(Float)
703    mesh['vertex_attributes'] = array(mesh['vertex_attributes']).astype(Float)
704    mesh['vertex_attribute_titles'] = array(mesh['vertex_attribute_titles']).astype(Character) 
705    mesh['segments'] = array(mesh['segments']).astype(IntType)
706    mesh['segment_tags'] = array(mesh['segment_tags']).astype(Character)
707    mesh['triangles'] = array(mesh['triangles']).astype(IntType)
708    mesh['triangle_tags'] = array(mesh['triangle_tags']) #.astype(Character)
709    mesh['triangle_neighbors'] = array(mesh['triangle_neighbors']).astype(IntType) 
710   
711    #the outline
712    mesh['points'] = array(mesh['points']).astype(Float)
713    mesh['point_attributes'] = array(mesh['point_attributes']).astype(Float)
714    mesh['outline_segments'] = array(mesh['outline_segments']).astype(IntType)
715    mesh['outline_segment_tags'] = array(mesh['outline_segment_tags']).astype(Character)
716    mesh['holes'] = array(mesh['holes']).astype(Float)
717    mesh['regions'] = array(mesh['regions']).astype(Float)
718    mesh['region_tags'] = array(mesh['region_tags']).astype(Character)
719    mesh['region_max_areas'] = array(mesh['region_max_areas']).astype(Float)
720   
721    #mesh = mesh_dict2array(mesh)
722    #print "new_mesh",new_mesh
723    #print "mesh",mesh
724   
725    # NetCDF file definition
726    outfile = NetCDFFile(filename, 'w')
727   
728    #Create new file
729    outfile.institution = 'Geoscience Australia'
730    outfile.description = 'NetCDF format for compact and portable storage ' +\
731                      'of spatial point data'
732   
733    # dimension definitions - fixed
734    outfile.createDimension('num_of_dimensions', 2) #This is 2d data
735    outfile.createDimension('num_of_segment_ends', 2) #Segs have two points
736    outfile.createDimension('num_of_triangle_vertices', 3)
737    outfile.createDimension('num_of_triangle_faces', 3)
738    outfile.createDimension('num_of_region_max_area', 1)
739
740    # Create dimensions, variables and set the variables
741   
742    # trianglulation
743    # vertices
744    if (mesh['vertices'].shape[0] > 0):
745        outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
746        outfile.createVariable('vertices', Float, ('num_of_vertices',
747                                                   'num_of_dimensions'))
748        outfile.variables['vertices'][:] = mesh['vertices']
749        if (mesh['vertex_attributes'].shape[0] > 0 and mesh['vertex_attributes'].shape[1] > 0):
750            outfile.createDimension('num_of_vertex_attributes',
751                                    mesh['vertex_attributes'].shape[1])
752            outfile.createDimension('num_of_vertex_attribute_title_chars',
753                                    mesh['vertex_attribute_titles'].shape[1])
754            outfile.createVariable('vertex_attributes',
755                                   Float,
756                                   ('num_of_vertices',
757                                    'num_of_vertex_attributes'))   
758            outfile.createVariable('vertex_attribute_titles',
759                                   Character,
760                                   ( 'num_of_vertex_attributes',
761                                     'num_of_vertex_attribute_title_chars' ))
762            outfile.variables['vertex_attributes'][:] = \
763                                                      mesh['vertex_attributes']
764            outfile.variables['vertex_attribute_titles'][:] = \
765                                     mesh['vertex_attribute_titles']
766    # segments
767    if (mesh['segments'].shape[0] > 0):       
768        outfile.createDimension('num_of_segments',
769                                mesh['segments'].shape[0])
770        outfile.createVariable('segments',
771                               IntType,
772                               ('num_of_segments', 'num_of_segment_ends'))
773        outfile.variables['segments'][:] = mesh['segments']
774        if (mesh['segment_tags'].shape[1] > 0):
775            outfile.createDimension('num_of_segment_tag_chars',
776                                    mesh['segment_tags'].shape[1])
777            outfile.createVariable('segment_tags',
778                                   Character,
779                                   ('num_of_segments',
780                                    'num_of_segment_tag_chars'))
781            outfile.variables['segment_tags'][:] = mesh['segment_tags']
782    # triangles   
783    if (mesh['triangles'].shape[0] > 0):
784        outfile.createDimension('num_of_triangles',
785                                mesh['triangles'].shape[0])
786        outfile.createVariable('triangles',
787                               IntType,
788                               ('num_of_triangles',
789                                'num_of_triangle_vertices'))
790        outfile.createVariable('triangle_neighbors',
791                               IntType,
792                               ('num_of_triangles',
793                                'num_of_triangle_faces'))
794        outfile.variables['triangles'][:] = mesh['triangles']
795        outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
796        if (mesh['triangle_tags'].shape[1] > 0):
797            outfile.createDimension('num_of_triangle_tag_chars',
798                                    mesh['triangle_tags'].shape[1]) 
799            outfile.createVariable('triangle_tags',
800                                   Character,
801                                   ('num_of_triangles',
802                                    'num_of_triangle_tag_chars'))
803            outfile.variables['triangle_tags'][:] = mesh['triangle_tags']
804   
805
806    # outline
807    # points
808    if (mesh['points'].shape[0] > 0):
809        outfile.createDimension('num_of_points', mesh['points'].shape[0])
810        outfile.createVariable('points', Float, ('num_of_points',
811                                                 'num_of_dimensions'))
812        outfile.variables['points'][:] = mesh['points'] 
813        if (mesh['point_attributes'].shape[0] > 0  and mesh['point_attributes'].shape[1] > 0):
814            outfile.createDimension('num_of_point_attributes',
815                                    mesh['point_attributes'].shape[1])
816            outfile.createVariable('point_attributes',
817                                   Float,
818                                   ('num_of_points',
819                                    'num_of_point_attributes'))
820            outfile.variables['point_attributes'][:] = mesh['point_attributes']
821    # outline_segments
822    if (mesh['outline_segments'].shape[0] > 0):
823        outfile.createDimension('num_of_outline_segments',
824                                mesh['outline_segments'].shape[0])
825        outfile.createVariable('outline_segments',
826                               IntType,
827                               ('num_of_outline_segments',
828                                'num_of_segment_ends'))
829        outfile.variables['outline_segments'][:] = mesh['outline_segments']
830        if (mesh['outline_segment_tags'].shape[1] > 0):
831            outfile.createDimension('num_of_outline_segment_tag_chars',
832                                    mesh['outline_segment_tags'].shape[1])
833            outfile.createVariable('outline_segment_tags',
834                                   Character,
835                                   ('num_of_outline_segments',
836                                    'num_of_outline_segment_tag_chars'))
837            outfile.variables['outline_segment_tags'][:] = mesh['outline_segment_tags']
838    # holes
839    if (mesh['holes'].shape[0] > 0):
840        outfile.createDimension('num_of_holes', mesh['holes'].shape[0])
841        outfile.createVariable('holes', Float, ('num_of_holes',
842                                                'num_of_dimensions'))
843        outfile.variables['holes'][:] = mesh['holes']
844    # regions
845    if (mesh['regions'].shape[0] > 0):
846        outfile.createDimension('num_of_regions', mesh['regions'].shape[0])
847        outfile.createVariable('regions', Float, ('num_of_regions',
848                                                  'num_of_dimensions'))
849        outfile.createVariable('region_max_areas',
850                               Float,
851                               ('num_of_regions',))
852        outfile.variables['regions'][:] = mesh['regions']
853        outfile.variables['region_max_areas'][:] = mesh['region_max_areas']
854        if (mesh['region_tags'].shape[1] > 0):
855            outfile.createDimension('num_of_region_tag_chars',
856                                    mesh['region_tags'].shape[1])
857            outfile.createVariable('region_tags',
858                                   Character,
859                                   ('num_of_regions',
860                                    'num_of_region_tag_chars'))
861            outfile.variables['region_tags'][:] = mesh['region_tags']
862
863    # geo_reference info
864    if mesh.has_key('geo_reference') and not mesh['geo_reference'] == None:
865        mesh['geo_reference'].write_NetCDF(outfile)
866                                                 
867    outfile.close()
868
869
870   
871def _read_msh_file(file_name):
872    """
873    Read in an msh file.
874
875    """
876       
877    from Scientific.IO.NetCDF import NetCDFFile         
878           
879    #Check contents
880    #Get NetCDF
881    fid = NetCDFFile(file_name, 'r') 
882
883    mesh = {}
884    # Get the variables
885    # the triangulation
886    try:
887        mesh['vertices'] = fid.variables['vertices'][:]
888    except KeyError:
889        mesh['vertices'] = array([])
890    try:
891        mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:]
892    except KeyError:
893        mesh['vertex_attributes'] = []
894        for ob in mesh['vertices']:
895            mesh['vertex_attributes'].append([])
896    mesh['vertex_attribute_titles'] = []
897    try:
898        titles = fid.variables['vertex_attribute_titles'][:]
899        for i, title in enumerate(titles):
900            mesh['vertex_attribute_titles'].append(titles[i].tostring().strip())
901    except KeyError:
902        pass 
903    try:
904        mesh['segments'] = fid.variables['segments'][:]
905    except KeyError:
906        mesh['segments'] = array([])
907    mesh['segment_tags'] =[]
908    try:
909        tags = fid.variables['segment_tags'][:]
910        for i, tag in enumerate(tags):
911            mesh['segment_tags'].append(tags[i].tostring().strip())
912    except KeyError:
913        for ob in mesh['segments']:
914            mesh['segment_tags'].append('')
915    try:
916        mesh['triangles'] = fid.variables['triangles'][:]
917        mesh['triangle_neighbors'] = fid.variables['triangle_neighbors'][:]
918    except KeyError:
919        mesh['triangles'] = array([])
920        mesh['triangle_neighbors'] = array([])
921    mesh['triangle_tags'] =[]
922    try:
923        tags = fid.variables['triangle_tags'][:]
924        for i, tag in enumerate(tags):
925            mesh['triangle_tags'].append(tags[i].tostring().strip())
926    except KeyError:
927        for ob in mesh['triangles']:
928            mesh['triangle_tags'].append('')
929   
930    #the outline
931    try:
932        mesh['points'] = fid.variables['points'][:]
933    except KeyError:
934        mesh['points'] = []
935    try:
936        mesh['point_attributes'] = fid.variables['point_attributes'][:]
937    except KeyError:
938        mesh['point_attributes'] = []
939        for point in mesh['points']:
940            mesh['point_attributes'].append([])
941    try:
942        mesh['outline_segments'] = fid.variables['outline_segments'][:]
943    except KeyError:
944        mesh['outline_segments'] = array([])
945    mesh['outline_segment_tags'] =[]
946    try:
947        tags = fid.variables['outline_segment_tags'][:]
948        for i, tag in enumerate(tags):
949            mesh['outline_segment_tags'].append(tags[i].tostring().strip())
950    except KeyError:
951        for ob in mesh['outline_segments']:
952            mesh['outline_segment_tags'].append('')
953    try:
954        mesh['holes'] = fid.variables['holes'][:]
955    except KeyError:
956        mesh['holes'] = array([])
957    try:
958        mesh['regions'] = fid.variables['regions'][:]
959    except KeyError:
960        mesh['regions'] = array([])
961    mesh['region_tags'] =[]
962    try:
963        tags = fid.variables['region_tags'][:]
964        for i, tag in enumerate(tags):
965            mesh['region_tags'].append(tags[i].tostring().strip())
966    except KeyError:
967        for ob in mesh['regions']:
968            mesh['region_tags'].append('')
969    try:
970        mesh['region_max_areas'] = fid.variables['region_max_areas'][:]
971    except KeyError:
972        mesh['region_max_areas'] = array([])
973    #mesh[''] = fid.variables[''][:]
974     
975    try:
976        geo_reference = Geo_reference(NetCDFObject=fid)
977        mesh['geo_reference'] = geo_reference
978    except AttributeError, e:
979        #geo_ref not compulsory
980        mesh['geo_reference'] = None
981   
982    fid.close()
983     
984    return mesh
985           
986
987## used by alpha shapes
988   
989def export_boundary_file( file_name, points, title, delimiter = ','):
990    """
991    export a file, ofile, with the format
992   
993    First line: Title variable
994    Following lines:  [point index][delimiter][point index]
995
996    file_name - the name of the new file
997    points - List of point index pairs [[p1, p2],[p3, p4]..]
998    title - info to write in the first line
999    """
1000   
1001    fd = open(file_name,'w')
1002   
1003    fd.write(title+"\n")
1004    #[point index][delimiter][point index]
1005    for point in points:
1006        fd.write( str(point[0]) + delimiter
1007                  + str(point[1]) + "\n")
1008    fd.close()
1009
1010
1011###
1012#  IMPORT/EXPORT XYA FILES
1013###
1014
1015def export_points_file(ofile,point_dict):
1016    """
1017    write a points file, ofile, as a binary (.msh) or text (.tsh) file
1018    """
1019    #this was done for all keys in the mesh file.
1020    #if not mesh_dict.has_key('points'):
1021    #    mesh_dict['points'] = []
1022    try:
1023        if (ofile[-4:] == ".xya"):
1024            _write_xya_file(ofile, point_dict)
1025        elif (ofile[-4:] == ".pts"):
1026            _write_pts_file(ofile, point_dict)
1027    except IOError, e:       
1028        msg = 'Could not write file %s ' %fileName
1029        raise IOError, msg
1030               
1031def import_points_file(ofile, delimiter = None, verbose = False):
1032    """
1033    load an .xya or .pts file   
1034    """
1035   
1036    if ofile[-4:]== ".xya":
1037        try:
1038            if delimiter == None:
1039                try:
1040                    fd = open(ofile)
1041                    points_dict = _read_xya_file(fd, ',')
1042                except SyntaxError:
1043                    fd.close()
1044                    fd = open(ofile)
1045                    points_dict = _read_xya_file(fd, ' ')
1046            else:
1047                fd = open(ofile)
1048                points_dict = _read_xya_file(fd, delimiter)
1049            fd.close()
1050            points_dict['geo_reference'] = None
1051        except IOError, e:       
1052            msg = 'Could not open file %s ' %ofile
1053            raise IOError, msg
1054       
1055    elif ofile[-4:]== ".pts":
1056        try:
1057            points_dict = _read_pts_file(ofile, verbose)       
1058        except IOError, e:       
1059            msg = 'Could not open file %s ' %ofile
1060            raise IOError, msg
1061    else:     
1062        msg = 'Extension %s is unknown' %ofile[-4:]
1063        raise IOError, msg
1064    return points_dict
1065
1066def extent_point_atts(point_atts):
1067    """
1068    Returns 4 points representing the extent
1069    This losses attribute info.
1070    """
1071    point_atts['pointlist'] = extent(point_atts['pointlist'])
1072    point_atts['attributelist'] = {}
1073    return point_atts
1074   
1075def extent(points):
1076    points = array(points).astype(Float)
1077    max_x = min_x = points[0][0]
1078    max_y = min_y = points[0][1]
1079    for point in points[1:]:
1080        x = point[0] 
1081        if x > max_x: max_x = x
1082        if x < min_x: min_x = x           
1083        y = point[1] 
1084        if y > max_y: max_y = y
1085        if y < min_y: min_y = y
1086    extent = array([[min_x,min_y],[max_x,min_y],[max_x,max_y],[min_x,max_y]])
1087    #print "extent",extent
1088    return extent
1089   
1090def reduce_pts(infile, outfile, max_points, verbose = False):
1091    point_atts = _read_pts_file(infile)
1092    while point_atts['pointlist'].shape[0] > max_points:
1093        if verbose: print "point_atts['pointlist'].shape[0]"
1094        point_atts = half_pts(point_atts)
1095    export_points_file(outfile, point_atts)
1096 
1097def produce_half_point_files(infile, max_points, delimiter, verbose = False):
1098    point_atts = _read_pts_file(infile)
1099    root, ext = splitext(infile)
1100    outfiles = []
1101    if verbose: print "# of points", point_atts['pointlist'].shape[0]
1102    while point_atts['pointlist'].shape[0] > max_points:
1103        point_atts = half_pts(point_atts)
1104        if verbose: print "# of points", point_atts['pointlist'].shape[0]
1105        outfile = root + delimiter + str(point_atts['pointlist'].shape[0]) + ext
1106        outfiles.append(outfile)
1107        export_points_file(outfile,
1108                  point_atts)
1109    return outfiles
1110   
1111def point_atts2array(point_atts):
1112    point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
1113   
1114    for key in point_atts['attributelist'].keys():
1115        point_atts['attributelist'][key]= array(point_atts['attributelist'][key]).astype(Float)
1116    return point_atts
1117   
1118def half_pts(point_atts):
1119    point_atts2array(point_atts)
1120    point_atts['pointlist'] = point_atts['pointlist'][::2]
1121   
1122    for key in point_atts['attributelist'].keys():
1123        point_atts['attributelist'][key]=  point_atts['attributelist'][key] [::2]
1124    return point_atts
1125
1126def concatinate_attributelist(dic):
1127    """
1128    giving a dic[attribute title] = attribute
1129    return list of attribute titles, array of attributes
1130    """
1131    point_attributes = array([]).astype(Float)
1132    keys = dic.keys()
1133    key = keys.pop(0)
1134    point_attributes = reshape(dic[key],(dic[key].shape[0],1))
1135    for key in keys:
1136        #point_attributes = concatenate([point_attributes, dic[key]], axis=1)
1137        reshaped = reshape(dic[key],(dic[key].shape[0],1))
1138        point_attributes = concatenate([point_attributes, reshaped], axis=1)
1139    return dic.keys(), point_attributes
1140
1141def _read_pts_file(filename, verbose = False):
1142    """Read .pts NetCDF file
1143   
1144    Return a dic of array of points, and dic of array of attribute
1145
1146    eg
1147    dic['points'] = [[1.0,2.0],[3.0,5.0]]
1148    dic['attributelist']['elevation'] = [[7.0,5.0]
1149    """
1150    #FIXME: (DSG)This format has issues.
1151    # There can't be an attribute called points
1152    # consider format change
1153
1154    from Scientific.IO.NetCDF import NetCDFFile
1155
1156    if verbose: print 'Reading ', filename
1157    fid = NetCDFFile(filename, 'r') 
1158
1159    point_atts = {} 
1160    # Get the variables
1161    point_atts['pointlist'] = array(fid.variables['points'])
1162    keys = fid.variables.keys()
1163    if verbose: print 'Got %d variables: %s' %(len(keys), keys)
1164    try:
1165        keys.remove('points')
1166    except IOError, e:       
1167        fid.close()   
1168        msg = 'Expected keyword "points" but could not find it'
1169        raise IOError, msg
1170
1171    attributes = {}
1172    for key in keys:
1173        if verbose: print "reading attribute '%s'" %key
1174       
1175        attributes[key] = array(fid.variables[key])
1176       
1177    point_atts['attributelist'] = attributes
1178   
1179    try:
1180        geo_reference = Geo_reference(NetCDFObject=fid)
1181        point_atts['geo_reference'] = geo_reference
1182    except AttributeError, e:
1183        #geo_ref not compulsory
1184        point_atts['geo_reference'] = None
1185   
1186   
1187    fid.close()
1188   
1189    #print "point_atts",point_atts
1190    return point_atts
1191
1192def _write_pts_file(filename, point_atts):
1193    """Write .pts NetCDF file   
1194
1195    WARNING: This function mangles the point_atts data structure
1196    """
1197    #FIXME: (DSG)This format has issues.
1198    # There can't be an attribute called points
1199    # consider format change
1200   
1201    from Scientific.IO.NetCDF import NetCDFFile
1202    point_atts2array(point_atts)
1203    # NetCDF file definition
1204    outfile = NetCDFFile(filename, 'w')
1205   
1206    #Create new file
1207    outfile.institution = 'Geoscience Australia'
1208    outfile.description = 'NetCDF format for compact and portable storage ' +\
1209                      'of spatial point data'
1210   
1211    # dimension definitions
1212    shape = point_atts['pointlist'].shape[0]
1213    outfile.createDimension('number_of_points', shape) 
1214    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1215   
1216    # variable definition
1217    outfile.createVariable('points', Float, ('number_of_points',
1218                                             'number_of_dimensions'))
1219
1220    #create variables 
1221    outfile.variables['points'][:] = point_atts['pointlist'] #.astype(Float32)
1222    for key in point_atts['attributelist'].keys():
1223        outfile.createVariable(key, Float, ('number_of_points',))
1224        outfile.variables[key][:] = point_atts['attributelist'][key] #.astype(Float32)
1225       
1226    if point_atts.has_key('geo_reference') and not point_atts['geo_reference'] == None:
1227        point_atts['geo_reference'].write_NetCDF(outfile)
1228       
1229    outfile.close() 
1230   
1231def _read_xya_file(fd, delimiter):
1232    lines = fd.readlines()
1233    points = []
1234    pointattributes = []
1235    if len(lines) <= 1:
1236        raise SyntaxError
1237    title = lines.pop(0) # the first (title) line
1238    att_names = clean_line(title,delimiter)
1239
1240    att_dict = {}
1241    for line in lines:
1242        #print "line >%s" %line
1243        numbers = clean_line(line,delimiter)
1244        #print "numbers >%s<" %numbers
1245        if len(numbers) < 2 and numbers != []:
1246            raise SyntaxError
1247        if numbers != []:
1248            try:
1249                x = float(numbers[0])
1250                y = float(numbers[1])
1251                points.append([x,y])
1252                numbers.pop(0)
1253                numbers.pop(0)
1254                #attributes = []
1255                #print "att_names",att_names
1256                #print "numbers",numbers
1257                if len(att_names) != len(numbers):
1258                    fd.close()
1259                    raise TitleAmountError
1260                   
1261                for i,num in enumerate(numbers):
1262                    num.strip()
1263                    if num != '\n' and num != '':
1264                        #attributes.append(float(num))
1265                        att_dict.setdefault(att_names[i],[]).append(float(num))
1266                       
1267            except ValueError:
1268                raise SyntaxError
1269            #pointattributes.append(attributes)
1270    xya_dict = {}
1271    xya_dict['pointlist'] = array(points).astype(Float)
1272   
1273    for key in att_dict.keys():
1274        att_dict[key] = array(att_dict[key]).astype(Float)
1275    xya_dict['attributelist'] = att_dict
1276
1277   
1278    ##print "xya_dict",xya_dict
1279    return xya_dict
1280
1281#FIXME(dsg), turn this dict plus methods into a class?
1282def take_points(dict,indices_to_keep):
1283    dict = point_atts2array(dict)
1284    #FIXME maybe the points data structure should become a class?
1285    dict['pointlist'] = take(dict['pointlist'],indices_to_keep)
1286
1287    for key in dict['attributelist'].keys():
1288        dict['attributelist'][key]= take(dict['attributelist'][key],
1289                                         indices_to_keep)
1290    return dict
1291   
1292def add_point_dictionaries (dict1, dict2):
1293    """
1294    """
1295    dict1 = point_atts2array(dict1)
1296    dict2 = point_atts2array(dict2)
1297   
1298    combined = {} 
1299    combined['pointlist'] = concatenate((dict2['pointlist'],
1300                                         dict1['pointlist']),axis=0)   
1301    atts = {}   
1302    for key in dict2['attributelist'].keys():
1303        atts[key]= concatenate((dict2['attributelist'][key],
1304                                dict1['attributelist'][key]), axis=0)
1305    combined['attributelist']=atts
1306    combined['geo_reference'] = dict1['geo_reference']
1307    return combined
1308   
1309#FIXME (DSG) need an export_points_file method..               
1310def export_xya_file( file_name, xya_dict, title, delimiter = ','):
1311    """
1312    export a file, ofile, with the format
1313   
1314    First line: Title variable
1315    Following lines:  <vertex #> <x> <y> [attributes]
1316
1317    file_name - the name of the new file
1318    xya_dict - point and point attribute info in a dictionary
1319    title - info to write in the first line
1320    """
1321    #FIXME, move the test for this from meshharness to loadasciiharness
1322    points = xya_dict['pointlist'] 
1323    pointattributes = xya_dict['pointattributelist']
1324   
1325    fd = open(file_name,'w')
1326   
1327    fd.write(title+"\n")
1328    #<vertex #> <x> <y> [attributes]
1329    for vert, vertatts in map(None, points, pointattributes):
1330        attlist = ""
1331        for att in vertatts:
1332            attlist = attlist + str(att)+ delimiter
1333        attlist = attlist[0:-len(delimiter)] # remove the last delimiter
1334        attlist.strip()
1335        fd.write( str(vert[0]) + delimiter
1336                  + str(vert[1]) + delimiter
1337                  + attlist + "\n")
1338    fd.close()
1339     
1340if __name__ == "__main__":
1341    pass
Note: See TracBrowser for help on using the repository browser.