source: anuga_core/source/anuga/load_mesh/loadASCII.py @ 4154

Last change on this file since 4154 was 4099, checked in by jakeman, 18 years ago

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