source: anuga_core/source/anuga/load_mesh/loadASCII.py @ 3850

Last change on this file since 3850 was 3850, checked in by ole, 17 years ago

Introduced a covariance measure to verify Okushiri timeseries plus
some incidental work

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