source: inundation/load_mesh/loadASCII.py @ 2942

Last change on this file since 2942 was 2941, checked in by duncan, 19 years ago

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