source: inundation/load_mesh/loadASCII.py @ 2897

Last change on this file since 2897 was 2888, checked in by duncan, 19 years ago

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