source: inundation/load_mesh/loadASCII.py @ 2267

Last change on this file since 2267 was 2262, checked in by ole, 19 years ago

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