source: branches/inundation-numpy-branch/load_mesh/loadASCII.py @ 7447

Last change on this file since 7447 was 2608, checked in by ole, 19 years ago

Played with custom exceptions for ANUGA
File size: 47.4 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 numpy 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 binary (.xya) or text (.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    point_atts = _read_pts_file(infile)
1060    while point_atts['pointlist'].shape[0] > max_points:
1061        if verbose: print "point_atts['pointlist'].shape[0]"
1062        point_atts = half_pts(point_atts)
1063    export_points_file(outfile, point_atts)
1064 
1065def produce_half_point_files(infile, max_points, delimiter, verbose = False):
1066    point_atts = _read_pts_file(infile)
1067    root, ext = splitext(infile)
1068    outfiles = []
1069    if verbose: print "# of points", point_atts['pointlist'].shape[0]
1070    while point_atts['pointlist'].shape[0] > max_points:
1071        point_atts = half_pts(point_atts)
1072        if verbose: print "# of points", point_atts['pointlist'].shape[0]
1073        outfile = root + delimiter + str(point_atts['pointlist'].shape[0]) + ext
1074        outfiles.append(outfile)
1075        export_points_file(outfile,
1076                  point_atts)
1077    return outfiles
1078   
1079def point_atts2array(point_atts):
1080    point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
1081   
1082    for key in point_atts['attributelist'].keys():
1083        point_atts['attributelist'][key]= array(point_atts['attributelist'][key]).astype(Float)
1084    return point_atts
1085   
1086def half_pts(point_atts):
1087    point_atts2array(point_atts)
1088    point_atts['pointlist'] = point_atts['pointlist'][::2]
1089   
1090    for key in point_atts['attributelist'].keys():
1091        point_atts['attributelist'][key]=  point_atts['attributelist'][key] [::2]
1092    return point_atts
1093
1094def concatinate_attributelist(dic):
1095    """
1096    giving a dic[attribute title] = attribute
1097    return list of attribute titles, array of attributes
1098    """
1099    point_attributes = array([]).astype(Float)
1100    keys = dic.keys()
1101    key = keys.pop(0)
1102    point_attributes = reshape(dic[key],(dic[key].shape[0],1))
1103    for key in keys:
1104        #point_attributes = concatenate([point_attributes, dic[key]], axis=1)
1105        reshaped = reshape(dic[key],(dic[key].shape[0],1))
1106        point_attributes = concatenate([point_attributes, reshaped], axis=1)
1107    return dic.keys(), point_attributes
1108
1109def _read_pts_file(file_name, verbose = False):
1110    """Read .pts NetCDF file
1111   
1112    Return a dic of array of points, and dic of array of attribute
1113
1114    eg
1115    dic['points'] = [[1.0,2.0],[3.0,5.0]]
1116    dic['attributelist']['elevation'] = [[7.0,5.0]
1117    """
1118    #FIXME: (DSG)This format has issues.
1119    # There can't be an attribute called points
1120    # consider format change
1121
1122    from Scientific.IO.NetCDF import NetCDFFile
1123
1124    if verbose: print 'Reading ', file_name
1125
1126   
1127    # see if the file is there.  Throw a QUIET IO error if it isn't
1128    fd = open(file_name,'r')
1129    fd.close()
1130
1131    #throws prints to screen if file not present
1132    fid = NetCDFFile(file_name, 'r') 
1133
1134    point_atts = {} 
1135    # Get the variables
1136    point_atts['pointlist'] = array(fid.variables['points'])
1137    keys = fid.variables.keys()
1138    if verbose: print 'Got %d variables: %s' %(len(keys), keys)
1139    try:
1140        keys.remove('points')
1141    except IOError, e:       
1142        fid.close()   
1143        msg = 'Expected keyword "points" but could not find it'
1144        raise IOError, msg
1145
1146    attributes = {}
1147    for key in keys:
1148        if verbose: print "reading attribute '%s'" %key
1149       
1150        attributes[key] = array(fid.variables[key])
1151       
1152    point_atts['attributelist'] = attributes
1153   
1154    try:
1155        geo_reference = Geo_reference(NetCDFObject=fid)
1156        point_atts['geo_reference'] = geo_reference
1157    except AttributeError, e:
1158        #geo_ref not compulsory
1159        point_atts['geo_reference'] = None
1160   
1161   
1162    fid.close()
1163   
1164    #print "point_atts",point_atts
1165    return point_atts
1166
1167def _write_pts_file(file_name, point_atts):
1168    """Write .pts NetCDF file   
1169
1170    WARNING: This function mangles the point_atts data structure
1171    """
1172    #FIXME: (DSG)This format has issues.
1173    # There can't be an attribute called points
1174    # consider format change
1175
1176    legal_keys = ['pointlist', 'attributelist', 'geo_reference']
1177    for key in point_atts.keys():
1178        msg = 'Key %s is illegal. Valid keys are %s' %(key, legal_keys) 
1179        assert key in legal_keys, msg
1180
1181    from Scientific.IO.NetCDF import NetCDFFile
1182    point_atts2array(point_atts)
1183    # NetCDF file definition
1184
1185    outfile = NetCDFFile(file_name, 'w')
1186
1187    #Create new file
1188    outfile.institution = 'Geoscience Australia'
1189    outfile.description = 'NetCDF format for compact and portable storage ' +\
1190                      'of spatial point data'
1191   
1192    # dimension definitions
1193    shape = point_atts['pointlist'].shape[0]
1194    outfile.createDimension('number_of_points', shape) 
1195    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1196   
1197    # variable definition
1198    outfile.createVariable('points', Float, ('number_of_points',
1199                                             'number_of_dimensions'))
1200
1201    #create variables 
1202    outfile.variables['points'][:] = point_atts['pointlist'] #.astype(Float32)
1203    for key in point_atts['attributelist'].keys():
1204        outfile.createVariable(key, Float, ('number_of_points',))
1205        outfile.variables[key][:] = point_atts['attributelist'][key] #.astype(Float32)
1206       
1207    if point_atts.has_key('geo_reference') and not point_atts['geo_reference'] == None:
1208        point_atts['geo_reference'].write_NetCDF(outfile)
1209       
1210    outfile.close() 
1211   
1212def _read_xya_file(fd, delimiter):
1213    #lines = fd.readlines()
1214    points = []
1215    pointattributes = []
1216    #if len(lines) <= 1:
1217    #    raise SyntaxError
1218    title = fd.readline()
1219    #title = lines.pop(0) # the first (title) line
1220    att_names = clean_line(title,delimiter)
1221
1222    att_dict = {}
1223    line = fd.readline()
1224    numbers = clean_line(line,delimiter)
1225    #for line in lines:
1226    while len(numbers) > 1:
1227        #print "line >%s" %line
1228        #numbers = clean_line(line,delimiter)
1229        #print "numbers >%s<" %numbers
1230        #if len(numbers) < 2 and numbers != []:
1231           
1232            # A line without two numbers
1233            # or a bad delimiter
1234            #FIXME dsg-dsg change error that is raised.
1235            #raise SyntaxError
1236        if numbers != []:
1237            try:
1238                x = float(numbers[0])
1239                y = float(numbers[1])
1240                points.append([x,y])
1241                numbers.pop(0)
1242                numbers.pop(0)
1243                #attributes = []
1244                #print "att_names",att_names
1245                #print "numbers",numbers
1246                if len(att_names) != len(numbers):
1247                    fd.close()
1248                    # It might not be a problem with the title
1249                    #raise TitleAmountError
1250                    raise IOError
1251                for i,num in enumerate(numbers):
1252                    num.strip()
1253                    if num != '\n' and num != '':
1254                        #attributes.append(float(num))
1255                        att_dict.setdefault(att_names[i],[]).append(float(num))
1256                   
1257            except ValueError:
1258                raise SyntaxError
1259        line = fd.readline()
1260        numbers = clean_line(line,delimiter)
1261       
1262    if line == '':
1263        # end of file
1264        geo_reference = None
1265    else:
1266        geo_reference = Geo_reference(ASCIIFile=fd,read_title=line)
1267   
1268    xya_dict = {}
1269    xya_dict['pointlist'] = array(points).astype(Float)
1270   
1271    for key in att_dict.keys():
1272        att_dict[key] = array(att_dict[key]).astype(Float)
1273    xya_dict['attributelist'] = att_dict
1274    xya_dict['geo_reference'] = geo_reference
1275    #print "xya_dict",xya_dict
1276    return xya_dict
1277
1278#FIXME(dsg), turn this dict plus methods into a class?
1279def take_points(dict,indices_to_keep):
1280    dict = point_atts2array(dict)
1281    #FIXME maybe the points data structure should become a class?
1282    dict['pointlist'] = take(dict['pointlist'],indices_to_keep)
1283
1284    for key in dict['attributelist'].keys():
1285        dict['attributelist'][key]= take(dict['attributelist'][key],
1286                                         indices_to_keep)
1287    return dict
1288   
1289def add_point_dictionaries (dict1, dict2):
1290    """
1291    """
1292    dict1 = point_atts2array(dict1)
1293    dict2 = point_atts2array(dict2)
1294   
1295    combined = {} 
1296    combined['pointlist'] = concatenate((dict2['pointlist'],
1297                                         dict1['pointlist']),axis=0)   
1298    atts = {}   
1299    for key in dict2['attributelist'].keys():
1300        atts[key]= concatenate((dict2['attributelist'][key],
1301                                dict1['attributelist'][key]), axis=0)
1302    combined['attributelist']=atts
1303    combined['geo_reference'] = dict1['geo_reference']
1304    return combined
1305                 
1306def _write_xya_file( file_name, xya_dict, delimiter = ','):
1307    """
1308    export a file, ofile, with the xya format
1309   
1310    """
1311    points = xya_dict['pointlist'] 
1312    pointattributes = xya_dict['attributelist']
1313   
1314    fd = open(file_name,'w')
1315 
1316    titlelist = ""
1317    for title in pointattributes.keys():
1318        titlelist = titlelist + title + delimiter
1319    titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
1320    fd.write(titlelist+"\n")
1321    #<vertex #> <x> <y> [attributes]
1322    for i,vert in enumerate( points):
1323       
1324        attlist = ","
1325        for att in pointattributes.keys():
1326            attlist = attlist + str(pointattributes[att][i])+ delimiter
1327        attlist = attlist[0:-len(delimiter)] # remove the last delimiter
1328        attlist.strip()
1329        fd.write( str(vert[0]) + delimiter
1330                  + str(vert[1])
1331                  + attlist + "\n")
1332   
1333    # geo_reference info
1334    if xya_dict.has_key('geo_reference') and \
1335           not xya_dict['geo_reference'] is None:
1336        xya_dict['geo_reference'].write_ASCII(fd)
1337    fd.close()
1338     
1339if __name__ == "__main__":
1340    pass
Note: See TracBrowser for help on using the repository browser.