source: inundation/load_mesh/loadASCII.py @ 3381

Last change on this file since 3381 was 3133, checked in by duncan, 19 years ago

small change
File size: 47.7 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, TitleError
85
86from Scientific.IO.NetCDF import NetCDFFile
87   
88import exceptions
89class TitleAmountError(exceptions.Exception): pass
90
91NOMAXAREA=-999
92
93# This is used by pmesh
94def import_mesh_file(ofile):
95    """
96    read a mesh file, either .tsh or .msh
97   
98    Note: will throw an IOError if it can't load the file.
99    Catch these!
100    """
101    try:
102        if ofile[-4:]== ".tsh":
103            dict = _read_tsh_file(ofile)
104        elif ofile[-4:]== ".msh":
105            dict = _read_msh_file(ofile)
106        else:     
107            msg = 'Extension .%s is unknown' %ofile[-4:]
108            raise IOError, msg
109    except (TitleError,SyntaxError,IndexError, ValueError): #FIXME No test for ValueError
110            msg = 'File could not be opened'
111            raise IOError, msg 
112    return dict
113
114
115def export_mesh_file(ofile,mesh_dict):
116    """
117    write a file, ofile, with the format
118   
119    First line:  <# of vertices> <# of attributes>
120    Following lines:  <vertex #> <x> <y> [attributes]
121    One line:  <# of triangles>
122    Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
123    One line:  <# of segments>
124    Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
125    """
126    #FIXME DSG-anyone: automate
127    if not mesh_dict.has_key('points'):
128        mesh_dict['points'] = []
129    if not mesh_dict.has_key('point_attributes'):
130        mesh_dict['point_attributes'] = []
131    if not mesh_dict.has_key('outline_segments'):
132        mesh_dict['outline_segments'] = []
133    if not mesh_dict.has_key('outline_segment_tags'):
134        mesh_dict['outline_segment_tags'] = []
135    if not mesh_dict.has_key('holes'):
136        mesh_dict['holes'] = []
137    if not mesh_dict.has_key('regions'):
138        mesh_dict['regions'] = []
139       
140    if not mesh_dict.has_key('region_tags'):
141        mesh_dict['region_tags'] = []
142    if not mesh_dict.has_key('region_max_areas'):
143        mesh_dict['region_max_areas'] = []
144    if not mesh_dict.has_key('vertices'):
145        mesh_dict['vertices'] = []
146    if not mesh_dict.has_key('vertex_attributes'):
147        mesh_dict['vertex_attributes'] = []
148    if not mesh_dict.has_key('vertex_attribute_titles'):
149        mesh_dict['vertex_attribute_titles'] = []
150    if not mesh_dict.has_key('segments'):
151        mesh_dict['segments'] = []
152    if not mesh_dict.has_key('segment_tags'):
153        mesh_dict['segment_tags'] = []
154    if not mesh_dict.has_key('triangles'):
155        mesh_dict['triangles'] = []
156    if not mesh_dict.has_key('triangle_tags'):
157        mesh_dict['triangle_tags'] = []
158    if not mesh_dict.has_key('triangle_neighbors'):
159        mesh_dict['triangle_neighbors'] = []
160    #print "DSG************"
161    #print "mesh_dict",mesh_dict
162    #print "DSG************"
163   
164    if (ofile[-4:] == ".tsh"):
165        _write_tsh_file(ofile,mesh_dict)
166    elif (ofile[-4:] == ".msh"):
167        _write_msh_file(ofile, mesh_dict)
168    else:
169        msg = 'Unknown file type %s ' %ofile
170        raise IOError, msg
171
172def _read_tsh_file(ofile):
173    """
174    Read the text file format for meshes
175    """
176    fd = open(ofile,'r')
177    dict = _read_triangulation(fd)
178    dict_mesh = _read_outline(fd)
179    for element in dict_mesh.keys():
180        dict[element] = dict_mesh[element]
181    fd.close()
182    return dict
183
184   
185def _read_triangulation(fd):
186    """
187    Read the generated triangulation, NOT the outline.
188    """
189    delimiter = " "
190    ######### loading the point info
191    line = fd.readline()
192    #print line
193    fragments = line.split()
194    #for fragment in fragments:
195    #print fragment
196    if fragments ==[]:
197        NumOfVertices = 0
198        NumOfVertAttributes = 0
199    else:
200        NumOfVertices = fragments[0]
201        NumOfVertAttributes = fragments[1]
202    points = []
203    pointattributes = []
204    for index in range(int(NumOfVertices)):       
205        #print index
206        fragments = fd.readline().split()
207        #print fragments
208        fragments.pop(0) #pop off the index
209       
210        # pop the x y off so we're left with a list of attributes       
211        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
212        points.append(vert)
213        apointattributes  = []
214        #print fragments
215        for fragment in fragments:
216            apointattributes.append(float(fragment))
217        pointattributes.append(apointattributes)
218       
219    ######### loading the point title info
220    line = fd.readline()
221    #print "point title comments",line
222    vertTitle = []
223    for index in range(int(NumOfVertAttributes)):       
224        #print index
225        fragments = fd.readline().strip() 
226        vertTitle.append(fragments)
227        #print vertTitle
228       
229    ######### loading the triangle info
230    line = fd.readline()
231    #print "triangle comments",line
232    fragments = line.split()
233    #for fragment in fragments:
234    #    print fragment
235    NumOfTriangles = fragments[0]
236    triangles = []
237    triangleattributes = []
238    triangleneighbors = []
239    for index in range(int(NumOfTriangles)):       
240        #print index
241        line = fd.readline()
242        line.strip() # so we can get the region string
243        fragments = line.split()
244        #print "triangle info", fragments
245        fragments.pop(0) #pop off the index
246       
247        tri = [int(fragments[0]),int(fragments[1]),int(fragments[2])]
248        triangles.append(tri)
249        neighbors = [int(fragments[3]),int(fragments[4]),int(fragments[5])]
250        triangleneighbors.append(neighbors)
251        for x in range(7): # remove index [<vertex #>] [<neigbouring tri #>]
252            line = line[find(line,delimiter):] # remove index
253            line = line.lstrip()
254        stringtag = line.strip()
255        triangleattributes.append(stringtag)
256       
257    ######### loading the segment info
258    line = fd.readline()
259    #print "seg comment line",line
260    fragments = line.split()
261    #for fragment in fragments:
262    #    print fragment
263    NumOfSegments = fragments[0]
264    segments = []
265    segmenttags = []
266    for index in range(int(NumOfSegments)):       
267        #print index
268        line = fd.readline()
269        line.strip() # to get the segment string
270        fragments = line.split()
271        #print fragments
272        fragments.pop(0) #pop off the index
273        seg = [int(fragments[0]),int(fragments[1])]
274        segments.append(seg)
275        line = line[find(line,delimiter):] # remove index
276        line = line.lstrip()
277        line = line[find(line,delimiter):] # remove x
278        line = line.lstrip()
279        line = line[find(line,delimiter):] # remove y
280        stringtag = line.strip()
281        segmenttags.append(stringtag)
282
283           
284    meshDict = {}
285    meshDict['vertices'] = points
286    meshDict['vertex_attributes'] = pointattributes
287    meshDict['triangles'] = triangles
288    meshDict['triangle_tags'] = triangleattributes
289    meshDict['triangle_neighbors'] = triangleneighbors
290    meshDict['segments'] = segments
291    meshDict['segment_tags'] = segmenttags
292    meshDict['vertex_attribute_titles'] = vertTitle
293   
294    return meshDict
295   
296def _read_outline(fd):
297    """
298    Note, if a file has no mesh info, it can still be read - the meshdic
299    returned will be 'empty'.
300    """
301    delimiter = " " # warning: split() calls are using default whitespace
302   
303    ######### loading the point info
304    line = fd.readline()
305    #print 'read_outline - line',line
306    fragments = line.split()
307    #for fragment in fragments:
308    if fragments ==[]:
309        NumOfVertices = 0
310        NumOfVertAttributes = 0
311    else:
312        NumOfVertices = fragments[0]
313        NumOfVertAttributes = fragments[1]
314    points = []
315    pointattributes = []
316    for index in range(int(NumOfVertices)):       
317        #print index
318        fragments = fd.readline().split() 
319        #print fragments
320        fragments.pop(0) #pop off the index
321        # pop the x y off so we're left with a list of attributes
322        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
323        points.append(vert)
324        apointattributes  = []
325        #print fragments
326        for fragment in fragments:
327            apointattributes.append(float(fragment))
328        pointattributes.append(apointattributes)
329       
330
331    ######### loading the segment info
332    line = fd.readline()
333    fragments = line.split()
334    #for fragment in fragments:
335    #    print fragment
336    if fragments ==[]:
337        NumOfSegments = 0
338    else:
339        NumOfSegments = fragments[0]
340    segments = []
341    segmenttags = []
342    for index in range(int(NumOfSegments)): 
343        #print index
344        line = fd.readline()
345        fragments = line.split()
346        #print fragments
347        fragments.pop(0) #pop off the index
348       
349        seg = [int(fragments[0]),int(fragments[1])]
350        segments.append(seg)
351        line = line[find(line,delimiter):] # remove index
352        line = line.lstrip()
353        line = line[find(line,delimiter):] # remove x
354        line = line.lstrip()
355        line = line[find(line,delimiter):] # remove y
356        stringtag = line.strip()
357        segmenttags.append(stringtag) 
358
359    ######### loading the hole info
360    line = fd.readline()
361    #print line
362    fragments = line.split()
363    #for fragment in fragments:
364    #    print fragment
365    if fragments ==[]:
366        numOfHoles = 0
367    else:
368        numOfHoles = fragments[0]
369    holes = []
370    for index in range(int(numOfHoles)):       
371        #print index
372        fragments = fd.readline().split()
373        #print fragments
374        fragments.pop(0) #pop off the index
375        hole = [float(fragments[0]),float(fragments[1])]
376        holes.append(hole)
377
378   
379    ######### loading the region info
380    line = fd.readline()
381    #print line
382    fragments = line.split()
383    #for fragment in fragments:
384    #    print fragment
385    if fragments ==[]:
386        numOfRegions = 0
387    else:
388        numOfRegions = fragments[0]
389    regions = []
390    regionattributes = []
391    for index in range(int(numOfRegions)):
392        line = fd.readline()
393        fragments = line.split()
394        #print fragments
395        fragments.pop(0) #pop off the index
396        region = [float(fragments[0]),float(fragments[1])]
397        regions.append(region)
398
399        line = line[find(line,delimiter):] # remove index
400        line = line.lstrip()
401        line = line[find(line,delimiter):] # remove x
402        line = line.lstrip()
403        line = line[find(line,delimiter):] # remove y
404        stringtag = line.strip()
405        regionattributes.append(stringtag)
406    regionmaxareas = []
407    line = fd.readline()
408    for index in range(int(numOfRegions)): # Read in the Max area info
409        line = fd.readline()
410        fragments = line.split()
411        #print fragments
412        # The try is here for format compatibility
413        try:
414            fragments.pop(0) #pop off the index
415            if len(fragments) == 0: #no max area
416                regionmaxareas.append(None)
417            else:
418                regionmaxareas.append(float(fragments[0]))
419        except IndexError, e:
420            regionmaxareas.append(None)
421
422    try:
423        geo_reference = Geo_reference(ASCIIFile=fd)
424    except:
425        #geo_ref not compulsory
426        geo_reference = None
427       
428    meshDict = {}
429    meshDict['points'] = points
430    meshDict['point_attributes'] = pointattributes
431    meshDict['outline_segments'] = segments
432    meshDict['outline_segment_tags'] = segmenttags
433    meshDict['holes'] = holes
434    meshDict['regions'] = regions
435    meshDict['region_tags'] = regionattributes
436    meshDict['region_max_areas'] = regionmaxareas
437    meshDict['geo_reference'] = geo_reference
438   
439    #print "meshDict",meshDict
440    return meshDict
441
442def clean_line(line,delimiter):     
443    """Remove whitespace
444    """
445    #print ">%s" %line
446    line = line.strip()
447    #print "stripped>%s" %line
448    numbers = line.split(delimiter)
449    i = len(numbers) - 1
450    while i >= 0:
451        if numbers[i] == '':
452            numbers.pop(i)
453        else:
454            numbers[i] = numbers[i].strip()
455       
456        i += -1
457    #for num in numbers:
458    #    print "num>%s<" %num
459    return numbers
460
461def _write_ASCII_triangulation(fd,
462                               gen_dict):
463    vertices = gen_dict['vertices']
464    vertices_attributes = gen_dict['vertex_attributes']
465    try:
466        vertices_attribute_titles = gen_dict['vertex_attribute_titles']   
467    except KeyError, e:
468        #FIXME is this the best way?
469        if vertices_attributes == [] or vertices_attributes[0] == []:
470             vertices_attribute_titles = []
471        else:
472            raise KeyError, e
473       
474    triangles = gen_dict['triangles']
475    triangles_attributes = gen_dict['triangle_tags']
476    triangle_neighbors = gen_dict['triangle_neighbors']
477       
478    segments = gen_dict['segments']
479    segment_tags = gen_dict['segment_tags']
480     
481    numVert = str(len(vertices))
482    if (numVert == "0") or len(vertices_attributes) == 0:
483        numVertAttrib = "0"
484    else:
485        numVertAttrib = str(len(vertices_attributes[0]))
486    fd.write(numVert + " " + numVertAttrib + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Triangulation Vertices..." + "\n")
487
488    #<vertex #> <x> <y> [attributes]   
489    index = 0 
490    for vert in vertices:
491        attlist = ""
492        for att in vertices_attributes[index]:
493            attlist = attlist + str(att)+" "
494        attlist.strip()
495        fd.write(str(index) + " "
496                 + str(vert[0]) + " "
497                 + str(vert[1]) + " "
498                 + attlist + "\n")
499        index += 1
500
501    # write comments for title
502    fd.write("# attribute column titles ...Triangulation Vertex Titles..." + "\n")
503    for title in vertices_attribute_titles:
504        fd.write(title + "\n")
505       
506    #<# of triangles>
507    n = len(triangles)
508    fd.write(str(n) + " # <# of triangles>, next lines <triangle #> [<vertex #>] [<neigbouring triangle #>] [attribute of region] ...Triangulation Triangles..." + "\n")
509       
510    # <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
511    for index in range(n):
512        neighbors = ""
513        tri = triangles[index]
514        for neighbor in triangle_neighbors[index]:
515            if neighbor:
516                neighbors += str(neighbor) + " "
517            else:
518                if neighbor == 0:
519                    neighbors +=  "0 "
520                else:
521                    neighbors +=  "-1 "
522        #Warning even though a list is past, only the first value
523        #is written.  There's an assumption that the list only
524        # contains one item. This assumption is made since the
525        #dict that's being passed around is also be used to communicate
526        # with triangle, and it seems to have the option of returning
527        # more than one value for triangle attributex
528        if triangles_attributes[index] == ['']:
529            att = ""
530        else:
531            att = str(triangles_attributes[index])
532        fd.write(str(index) + " "
533                 + str(tri[0]) + " " 
534                 + str(tri[1]) + " " 
535                 + str(tri[2]) + " " 
536                 + neighbors + " "
537                 + att + "\n")
538           
539    #One line:  <# of segments>
540    fd.write(str(len(segments)) + 
541             " # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Triangulation Segments..." + "\n")
542       
543    #Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
544    for i in range(len(segments)):
545        seg = segments[i]
546        fd.write(str(i) + " "
547                 + str(seg[0]) + " " 
548                 + str(seg[1]) + " " 
549                 + str(segment_tags[i]) + "\n")
550
551
552def _write_tsh_file(ofile,mesh_dict):
553            fd = open(ofile,'w')
554            _write_ASCII_triangulation(fd,mesh_dict)
555            _write_ASCII_outline(fd,mesh_dict)   
556            fd.close()
557
558def _write_ASCII_outline(fd,
559                     dict):
560    points = dict['points']
561    point_attributes = dict['point_attributes']
562    segments = dict['outline_segments']
563    segment_tags = dict['outline_segment_tags']
564    holes = dict['holes']
565    regions = dict['regions']
566    region_tags = dict['region_tags']
567    region_max_areas = dict['region_max_areas']
568       
569    num_points = str(len(points))
570    if (num_points == "0"):
571        num_point_atts = "0"
572    else:
573        num_point_atts = str(len(point_attributes[0]))
574    fd.write(num_points + " " + num_point_atts + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Mesh Vertices..." + "\n")
575       
576    # <x> <y> [attributes]
577    for i,point in enumerate(points):
578        attlist = ""
579        for att in point_attributes[i]:
580            attlist = attlist + str(att)+" "
581        attlist.strip()
582        fd.write(str(i) + " "
583                 +str(point[0]) + " "
584                 + str(point[1]) + " "
585                 + attlist + "\n")
586           
587    #One line:  <# of segments>
588    fd.write(str(len(segments)) + 
589             " # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Mesh Segments..." + "\n")
590       
591    #Following lines:  <vertex #>  <vertex #> [boundary tag]
592    for i,seg in enumerate(segments):
593        fd.write(str(i) + " "
594                 + str(seg[0]) + " " 
595                 + str(seg[1]) + " " 
596                 + str(segment_tags[i]) + "\n")
597
598    #One line:  <# of holes>
599    fd.write(str(len(holes)) + 
600             " # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes..." + "\n")   
601    # <x> <y>
602    for i,h in enumerate(holes):
603        fd.write(str(i) + " "
604                 + str(h[0]) + " "
605                 + str(h[1]) + "\n")
606       
607    #One line:  <# of regions>
608    fd.write(str(len(regions)) + 
609             " # <# of regions>, next lines <Region #> <x> <y> <tag>...Mesh Regions..." + "\n")   
610    # <index> <x> <y> <tag>
611    #print "regions",regions
612    for i,r in enumerate(regions):
613        #print "r",r
614        fd.write(str(i) + " " 
615                 + str(r[0]) + " "
616                 + str(r[1])+ " "
617                 + str(region_tags[i]) + "\n")
618       
619    # <index> [<MaxArea>|""]
620       
621    #One line:  <# of regions>
622    fd.write(str(len(regions)) + 
623             " # <# of regions>, next lines <Region #> [Max Area] ...Mesh Regions..." + "\n") 
624    for i,r in enumerate(regions):
625        area = str(region_max_areas[i])
626           
627        fd.write(str(i) + " " + area + "\n")
628
629    # geo_reference info
630    if dict.has_key('geo_reference') and not dict['geo_reference'] is None:
631        dict['geo_reference'].write_ASCII(fd)
632
633def _write_msh_file(file_name, mesh):
634    """Write .msh NetCDF file   
635
636    WARNING: This function mangles the mesh data structure
637    """
638 
639
640    IntType = Int32
641   
642    #
643    #the triangulation
644    mesh['vertices'] = array(mesh['vertices']).astype(Float)
645    mesh['vertex_attributes'] = array(mesh['vertex_attributes']).astype(Float)
646    mesh['vertex_attribute_titles'] = array(mesh['vertex_attribute_titles']).astype(Character) 
647    mesh['segments'] = array(mesh['segments']).astype(IntType)
648    mesh['segment_tags'] = array(mesh['segment_tags']).astype(Character)
649    mesh['triangles'] = array(mesh['triangles']).astype(IntType)
650    mesh['triangle_tags'] = array(mesh['triangle_tags']) #.astype(Character)
651    mesh['triangle_neighbors'] = array(mesh['triangle_neighbors']).astype(IntType) 
652   
653    #the outline
654    mesh['points'] = array(mesh['points']).astype(Float)
655    mesh['point_attributes'] = array(mesh['point_attributes']).astype(Float)
656    mesh['outline_segments'] = array(mesh['outline_segments']).astype(IntType)
657    mesh['outline_segment_tags'] = array(mesh['outline_segment_tags']).astype(Character)
658    mesh['holes'] = array(mesh['holes']).astype(Float)
659    mesh['regions'] = array(mesh['regions']).astype(Float)
660    mesh['region_tags'] = array(mesh['region_tags']).astype(Character)
661    mesh['region_max_areas'] = array(mesh['region_max_areas']).astype(Float)
662   
663    #mesh = mesh_dict2array(mesh)
664    #print "new_mesh",new_mesh
665    #print "mesh",mesh
666   
667    # NetCDF file definition
668    try:
669        outfile = NetCDFFile(file_name, 'w')
670    except IOError:
671        msg = 'File %s could not be created' %file_name
672        raise msg
673       
674    #Create new file
675    outfile.institution = 'Geoscience Australia'
676    outfile.description = 'NetCDF format for compact and portable storage ' +\
677                      'of spatial point data'
678   
679    # dimension definitions - fixed
680    outfile.createDimension('num_of_dimensions', 2) #This is 2d data
681    outfile.createDimension('num_of_segment_ends', 2) #Segs have two points
682    outfile.createDimension('num_of_triangle_vertices', 3)
683    outfile.createDimension('num_of_triangle_faces', 3)
684    outfile.createDimension('num_of_region_max_area', 1)
685
686    # Create dimensions, variables and set the variables
687   
688    # trianglulation
689    # vertices
690    if (mesh['vertices'].shape[0] > 0):
691        outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
692        outfile.createVariable('vertices', Float, ('num_of_vertices',
693                                                   'num_of_dimensions'))
694        outfile.variables['vertices'][:] = mesh['vertices']
695        if (mesh['vertex_attributes'].shape[0] > 0 and mesh['vertex_attributes'].shape[1] > 0):
696            outfile.createDimension('num_of_vertex_attributes',
697                                    mesh['vertex_attributes'].shape[1])
698            outfile.createDimension('num_of_vertex_attribute_title_chars',
699                                    mesh['vertex_attribute_titles'].shape[1])
700            outfile.createVariable('vertex_attributes',
701                                   Float,
702                                   ('num_of_vertices',
703                                    'num_of_vertex_attributes'))   
704            outfile.createVariable('vertex_attribute_titles',
705                                   Character,
706                                   ( 'num_of_vertex_attributes',
707                                     'num_of_vertex_attribute_title_chars' ))
708            outfile.variables['vertex_attributes'][:] = \
709                                                      mesh['vertex_attributes']
710            outfile.variables['vertex_attribute_titles'][:] = \
711                                     mesh['vertex_attribute_titles']
712    # segments
713    if (mesh['segments'].shape[0] > 0):       
714        outfile.createDimension('num_of_segments',
715                                mesh['segments'].shape[0])
716        outfile.createVariable('segments',
717                               IntType,
718                               ('num_of_segments', 'num_of_segment_ends'))
719        outfile.variables['segments'][:] = mesh['segments']
720        if (mesh['segment_tags'].shape[1] > 0):
721            outfile.createDimension('num_of_segment_tag_chars',
722                                    mesh['segment_tags'].shape[1])
723            outfile.createVariable('segment_tags',
724                                   Character,
725                                   ('num_of_segments',
726                                    'num_of_segment_tag_chars'))
727            outfile.variables['segment_tags'][:] = mesh['segment_tags']
728    # triangles   
729    if (mesh['triangles'].shape[0] > 0):
730        outfile.createDimension('num_of_triangles',
731                                mesh['triangles'].shape[0])
732        outfile.createVariable('triangles',
733                               IntType,
734                               ('num_of_triangles',
735                                'num_of_triangle_vertices'))
736        outfile.createVariable('triangle_neighbors',
737                               IntType,
738                               ('num_of_triangles',
739                                'num_of_triangle_faces'))
740        outfile.variables['triangles'][:] = mesh['triangles']
741        outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
742        if (mesh['triangle_tags'].shape[1] > 0):
743            outfile.createDimension('num_of_triangle_tag_chars',
744                                    mesh['triangle_tags'].shape[1]) 
745            outfile.createVariable('triangle_tags',
746                                   Character,
747                                   ('num_of_triangles',
748                                    'num_of_triangle_tag_chars'))
749            outfile.variables['triangle_tags'][:] = mesh['triangle_tags']
750   
751
752    # outline
753    # points
754    if (mesh['points'].shape[0] > 0):
755        outfile.createDimension('num_of_points', mesh['points'].shape[0])
756        outfile.createVariable('points', Float, ('num_of_points',
757                                                 'num_of_dimensions'))
758        outfile.variables['points'][:] = mesh['points'] 
759        if (mesh['point_attributes'].shape[0] > 0  and mesh['point_attributes'].shape[1] > 0):
760            outfile.createDimension('num_of_point_attributes',
761                                    mesh['point_attributes'].shape[1])
762            outfile.createVariable('point_attributes',
763                                   Float,
764                                   ('num_of_points',
765                                    'num_of_point_attributes'))
766            outfile.variables['point_attributes'][:] = mesh['point_attributes']
767    # outline_segments
768    if (mesh['outline_segments'].shape[0] > 0):
769        outfile.createDimension('num_of_outline_segments',
770                                mesh['outline_segments'].shape[0])
771        outfile.createVariable('outline_segments',
772                               IntType,
773                               ('num_of_outline_segments',
774                                'num_of_segment_ends'))
775        outfile.variables['outline_segments'][:] = mesh['outline_segments']
776        if (mesh['outline_segment_tags'].shape[1] > 0):
777            outfile.createDimension('num_of_outline_segment_tag_chars',
778                                    mesh['outline_segment_tags'].shape[1])
779            outfile.createVariable('outline_segment_tags',
780                                   Character,
781                                   ('num_of_outline_segments',
782                                    'num_of_outline_segment_tag_chars'))
783            outfile.variables['outline_segment_tags'][:] = mesh['outline_segment_tags']
784    # holes
785    if (mesh['holes'].shape[0] > 0):
786        outfile.createDimension('num_of_holes', mesh['holes'].shape[0])
787        outfile.createVariable('holes', Float, ('num_of_holes',
788                                                'num_of_dimensions'))
789        outfile.variables['holes'][:] = mesh['holes']
790    # regions
791    if (mesh['regions'].shape[0] > 0):
792        outfile.createDimension('num_of_regions', mesh['regions'].shape[0])
793        outfile.createVariable('regions', Float, ('num_of_regions',
794                                                  'num_of_dimensions'))
795        outfile.createVariable('region_max_areas',
796                               Float,
797                               ('num_of_regions',))
798        outfile.variables['regions'][:] = mesh['regions']
799        outfile.variables['region_max_areas'][:] = mesh['region_max_areas']
800        if (mesh['region_tags'].shape[1] > 0):
801            outfile.createDimension('num_of_region_tag_chars',
802                                    mesh['region_tags'].shape[1])
803            outfile.createVariable('region_tags',
804                                   Character,
805                                   ('num_of_regions',
806                                    'num_of_region_tag_chars'))
807            outfile.variables['region_tags'][:] = mesh['region_tags']
808
809    # geo_reference info
810    if mesh.has_key('geo_reference') and not mesh['geo_reference'] == None:
811        mesh['geo_reference'].write_NetCDF(outfile)
812                                                 
813    outfile.close()
814
815
816   
817def _read_msh_file(file_name):
818    """
819    Read in an msh file.
820
821    """
822               
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    #No.  Where this function is used, geospatial objects should be used.
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.