source: inundation/load_mesh/loadASCII.py @ 2253

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

creating flatter structure

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