source: inundation/load_mesh/loadASCII.py @ 2257

Last change on this file since 2257 was 2257, checked in by duncan, 18 years ago

comments

File size: 47.1 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    from Scientific.IO.NetCDF import NetCDFFile
1174    point_atts2array(point_atts)
1175    # NetCDF file definition
1176    outfile = NetCDFFile(file_name, 'w')
1177   
1178    #Create new file
1179    outfile.institution = 'Geoscience Australia'
1180    outfile.description = 'NetCDF format for compact and portable storage ' +\
1181                      'of spatial point data'
1182   
1183    # dimension definitions
1184    shape = point_atts['pointlist'].shape[0]
1185    outfile.createDimension('number_of_points', shape) 
1186    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1187   
1188    # variable definition
1189    outfile.createVariable('points', Float, ('number_of_points',
1190                                             'number_of_dimensions'))
1191
1192    #create variables 
1193    outfile.variables['points'][:] = point_atts['pointlist'] #.astype(Float32)
1194    for key in point_atts['attributelist'].keys():
1195        outfile.createVariable(key, Float, ('number_of_points',))
1196        outfile.variables[key][:] = point_atts['attributelist'][key] #.astype(Float32)
1197       
1198    if point_atts.has_key('geo_reference') and not point_atts['geo_reference'] == None:
1199        point_atts['geo_reference'].write_NetCDF(outfile)
1200       
1201    outfile.close() 
1202   
1203def _read_xya_file(fd, delimiter):
1204    #lines = fd.readlines()
1205    points = []
1206    pointattributes = []
1207    #if len(lines) <= 1:
1208    #    raise SyntaxError
1209    title = fd.readline()
1210    #title = lines.pop(0) # the first (title) line
1211    att_names = clean_line(title,delimiter)
1212
1213    att_dict = {}
1214    line = fd.readline()
1215    numbers = clean_line(line,delimiter)
1216    #for line in lines:
1217    while len(numbers) > 1:
1218        #print "line >%s" %line
1219        #numbers = clean_line(line,delimiter)
1220        #print "numbers >%s<" %numbers
1221        #if len(numbers) < 2 and numbers != []:
1222           
1223            # A line without two numbers
1224            # or a bad delimiter
1225            #FIXME dsg-dsg change error that is raised.
1226            #raise SyntaxError
1227        if numbers != []:
1228            try:
1229                x = float(numbers[0])
1230                y = float(numbers[1])
1231                points.append([x,y])
1232                numbers.pop(0)
1233                numbers.pop(0)
1234                #attributes = []
1235                #print "att_names",att_names
1236                #print "numbers",numbers
1237                if len(att_names) != len(numbers):
1238                    fd.close()
1239                    # It might not be a problem with the title
1240                    #raise TitleAmountError
1241                    raise IOError
1242                for i,num in enumerate(numbers):
1243                    num.strip()
1244                    if num != '\n' and num != '':
1245                        #attributes.append(float(num))
1246                        att_dict.setdefault(att_names[i],[]).append(float(num))
1247                   
1248            except ValueError:
1249                raise SyntaxError
1250        line = fd.readline()
1251        numbers = clean_line(line,delimiter)
1252       
1253    if line == '':
1254        # end of file
1255        geo_reference = None
1256    else:
1257        geo_reference = Geo_reference(ASCIIFile=fd,read_title=line)
1258   
1259    xya_dict = {}
1260    xya_dict['pointlist'] = array(points).astype(Float)
1261   
1262    for key in att_dict.keys():
1263        att_dict[key] = array(att_dict[key]).astype(Float)
1264    xya_dict['attributelist'] = att_dict
1265    xya_dict['geo_reference'] = geo_reference
1266    #print "xya_dict",xya_dict
1267    return xya_dict
1268
1269#FIXME(dsg), turn this dict plus methods into a class?
1270def take_points(dict,indices_to_keep):
1271    dict = point_atts2array(dict)
1272    #FIXME maybe the points data structure should become a class?
1273    dict['pointlist'] = take(dict['pointlist'],indices_to_keep)
1274
1275    for key in dict['attributelist'].keys():
1276        dict['attributelist'][key]= take(dict['attributelist'][key],
1277                                         indices_to_keep)
1278    return dict
1279   
1280def add_point_dictionaries (dict1, dict2):
1281    """
1282    """
1283    dict1 = point_atts2array(dict1)
1284    dict2 = point_atts2array(dict2)
1285   
1286    combined = {} 
1287    combined['pointlist'] = concatenate((dict2['pointlist'],
1288                                         dict1['pointlist']),axis=0)   
1289    atts = {}   
1290    for key in dict2['attributelist'].keys():
1291        atts[key]= concatenate((dict2['attributelist'][key],
1292                                dict1['attributelist'][key]), axis=0)
1293    combined['attributelist']=atts
1294    combined['geo_reference'] = dict1['geo_reference']
1295    return combined
1296                 
1297def _write_xya_file( file_name, xya_dict, delimiter = ','):
1298    """
1299    export a file, ofile, with the xya format
1300   
1301    """
1302    points = xya_dict['pointlist'] 
1303    pointattributes = xya_dict['attributelist']
1304   
1305    fd = open(file_name,'w')
1306 
1307    titlelist = ""
1308    for title in pointattributes.keys():
1309        titlelist = titlelist + title + delimiter
1310    titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
1311    fd.write(titlelist+"\n")
1312    #<vertex #> <x> <y> [attributes]
1313    for i,vert in enumerate( points):
1314       
1315        attlist = ","
1316        for att in pointattributes.keys():
1317            attlist = attlist + str(pointattributes[att][i])+ delimiter
1318        attlist = attlist[0:-len(delimiter)] # remove the last delimiter
1319        attlist.strip()
1320        fd.write( str(vert[0]) + delimiter
1321                  + str(vert[1])
1322                  + attlist + "\n")
1323   
1324    # geo_reference info
1325    if xya_dict.has_key('geo_reference') and \
1326           not xya_dict['geo_reference'] is None:
1327        xya_dict['geo_reference'].write_ASCII(fd)
1328    fd.close()
1329     
1330if __name__ == "__main__":
1331    pass
Note: See TracBrowser for help on using the repository browser.