source: inundation/ga/storm_surge/pmesh/load_mesh/loadASCII.py @ 1668

Last change on this file since 1668 was 1657, checked in by ole, 20 years ago

Comments and work on benchmark 2 (lwru2.py)
New unit test of least squares populating domain
File size: 46.5 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    outfile = NetCDFFile(file_name, 'w')
664   
665    #Create new file
666    outfile.institution = 'Geoscience Australia'
667    outfile.description = 'NetCDF format for compact and portable storage ' +\
668                      'of spatial point data'
669   
670    # dimension definitions - fixed
671    outfile.createDimension('num_of_dimensions', 2) #This is 2d data
672    outfile.createDimension('num_of_segment_ends', 2) #Segs have two points
673    outfile.createDimension('num_of_triangle_vertices', 3)
674    outfile.createDimension('num_of_triangle_faces', 3)
675    outfile.createDimension('num_of_region_max_area', 1)
676
677    # Create dimensions, variables and set the variables
678   
679    # trianglulation
680    # vertices
681    if (mesh['vertices'].shape[0] > 0):
682        outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
683        outfile.createVariable('vertices', Float, ('num_of_vertices',
684                                                   'num_of_dimensions'))
685        outfile.variables['vertices'][:] = mesh['vertices']
686        if (mesh['vertex_attributes'].shape[0] > 0 and mesh['vertex_attributes'].shape[1] > 0):
687            outfile.createDimension('num_of_vertex_attributes',
688                                    mesh['vertex_attributes'].shape[1])
689            outfile.createDimension('num_of_vertex_attribute_title_chars',
690                                    mesh['vertex_attribute_titles'].shape[1])
691            outfile.createVariable('vertex_attributes',
692                                   Float,
693                                   ('num_of_vertices',
694                                    'num_of_vertex_attributes'))   
695            outfile.createVariable('vertex_attribute_titles',
696                                   Character,
697                                   ( 'num_of_vertex_attributes',
698                                     'num_of_vertex_attribute_title_chars' ))
699            outfile.variables['vertex_attributes'][:] = \
700                                                      mesh['vertex_attributes']
701            outfile.variables['vertex_attribute_titles'][:] = \
702                                     mesh['vertex_attribute_titles']
703    # segments
704    if (mesh['segments'].shape[0] > 0):       
705        outfile.createDimension('num_of_segments',
706                                mesh['segments'].shape[0])
707        outfile.createVariable('segments',
708                               IntType,
709                               ('num_of_segments', 'num_of_segment_ends'))
710        outfile.variables['segments'][:] = mesh['segments']
711        if (mesh['segment_tags'].shape[1] > 0):
712            outfile.createDimension('num_of_segment_tag_chars',
713                                    mesh['segment_tags'].shape[1])
714            outfile.createVariable('segment_tags',
715                                   Character,
716                                   ('num_of_segments',
717                                    'num_of_segment_tag_chars'))
718            outfile.variables['segment_tags'][:] = mesh['segment_tags']
719    # triangles   
720    if (mesh['triangles'].shape[0] > 0):
721        outfile.createDimension('num_of_triangles',
722                                mesh['triangles'].shape[0])
723        outfile.createVariable('triangles',
724                               IntType,
725                               ('num_of_triangles',
726                                'num_of_triangle_vertices'))
727        outfile.createVariable('triangle_neighbors',
728                               IntType,
729                               ('num_of_triangles',
730                                'num_of_triangle_faces'))
731        outfile.variables['triangles'][:] = mesh['triangles']
732        outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
733        if (mesh['triangle_tags'].shape[1] > 0):
734            outfile.createDimension('num_of_triangle_tag_chars',
735                                    mesh['triangle_tags'].shape[1]) 
736            outfile.createVariable('triangle_tags',
737                                   Character,
738                                   ('num_of_triangles',
739                                    'num_of_triangle_tag_chars'))
740            outfile.variables['triangle_tags'][:] = mesh['triangle_tags']
741   
742
743    # outline
744    # points
745    if (mesh['points'].shape[0] > 0):
746        outfile.createDimension('num_of_points', mesh['points'].shape[0])
747        outfile.createVariable('points', Float, ('num_of_points',
748                                                 'num_of_dimensions'))
749        outfile.variables['points'][:] = mesh['points'] 
750        if (mesh['point_attributes'].shape[0] > 0  and mesh['point_attributes'].shape[1] > 0):
751            outfile.createDimension('num_of_point_attributes',
752                                    mesh['point_attributes'].shape[1])
753            outfile.createVariable('point_attributes',
754                                   Float,
755                                   ('num_of_points',
756                                    'num_of_point_attributes'))
757            outfile.variables['point_attributes'][:] = mesh['point_attributes']
758    # outline_segments
759    if (mesh['outline_segments'].shape[0] > 0):
760        outfile.createDimension('num_of_outline_segments',
761                                mesh['outline_segments'].shape[0])
762        outfile.createVariable('outline_segments',
763                               IntType,
764                               ('num_of_outline_segments',
765                                'num_of_segment_ends'))
766        outfile.variables['outline_segments'][:] = mesh['outline_segments']
767        if (mesh['outline_segment_tags'].shape[1] > 0):
768            outfile.createDimension('num_of_outline_segment_tag_chars',
769                                    mesh['outline_segment_tags'].shape[1])
770            outfile.createVariable('outline_segment_tags',
771                                   Character,
772                                   ('num_of_outline_segments',
773                                    'num_of_outline_segment_tag_chars'))
774            outfile.variables['outline_segment_tags'][:] = mesh['outline_segment_tags']
775    # holes
776    if (mesh['holes'].shape[0] > 0):
777        outfile.createDimension('num_of_holes', mesh['holes'].shape[0])
778        outfile.createVariable('holes', Float, ('num_of_holes',
779                                                'num_of_dimensions'))
780        outfile.variables['holes'][:] = mesh['holes']
781    # regions
782    if (mesh['regions'].shape[0] > 0):
783        outfile.createDimension('num_of_regions', mesh['regions'].shape[0])
784        outfile.createVariable('regions', Float, ('num_of_regions',
785                                                  'num_of_dimensions'))
786        outfile.createVariable('region_max_areas',
787                               Float,
788                               ('num_of_regions',))
789        outfile.variables['regions'][:] = mesh['regions']
790        outfile.variables['region_max_areas'][:] = mesh['region_max_areas']
791        if (mesh['region_tags'].shape[1] > 0):
792            outfile.createDimension('num_of_region_tag_chars',
793                                    mesh['region_tags'].shape[1])
794            outfile.createVariable('region_tags',
795                                   Character,
796                                   ('num_of_regions',
797                                    'num_of_region_tag_chars'))
798            outfile.variables['region_tags'][:] = mesh['region_tags']
799
800    # geo_reference info
801    if mesh.has_key('geo_reference') and not mesh['geo_reference'] == None:
802        mesh['geo_reference'].write_NetCDF(outfile)
803                                                 
804    outfile.close()
805
806
807   
808def _read_msh_file(file_name):
809    """
810    Read in an msh file.
811
812    """
813       
814    from Scientific.IO.NetCDF import NetCDFFile         
815           
816    #Check contents
817    #Get NetCDF
818   
819    # see if the file is there.  Throw a QUIET IO error if it isn't
820    fd = open(file_name,'r')
821    fd.close()
822
823    #throws prints to screen if file not present
824    fid = NetCDFFile(file_name, 'r') 
825
826    mesh = {}
827    # Get the variables
828    # the triangulation
829    try:
830        mesh['vertices'] = fid.variables['vertices'][:]
831    except KeyError:
832        mesh['vertices'] = array([])
833    try:
834        mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:]
835    except KeyError:
836        mesh['vertex_attributes'] = []
837        for ob in mesh['vertices']:
838            mesh['vertex_attributes'].append([])
839    mesh['vertex_attribute_titles'] = []
840    try:
841        titles = fid.variables['vertex_attribute_titles'][:]
842        for i, title in enumerate(titles):
843            mesh['vertex_attribute_titles'].append(titles[i].tostring().strip())
844    except KeyError:
845        pass 
846    try:
847        mesh['segments'] = fid.variables['segments'][:]
848    except KeyError:
849        mesh['segments'] = array([])
850    mesh['segment_tags'] =[]
851    try:
852        tags = fid.variables['segment_tags'][:]
853        for i, tag in enumerate(tags):
854            mesh['segment_tags'].append(tags[i].tostring().strip())
855    except KeyError:
856        for ob in mesh['segments']:
857            mesh['segment_tags'].append('')
858    try:
859        mesh['triangles'] = fid.variables['triangles'][:]
860        mesh['triangle_neighbors'] = fid.variables['triangle_neighbors'][:]
861    except KeyError:
862        mesh['triangles'] = array([])
863        mesh['triangle_neighbors'] = array([])
864    mesh['triangle_tags'] =[]
865    try:
866        tags = fid.variables['triangle_tags'][:]
867        for i, tag in enumerate(tags):
868            mesh['triangle_tags'].append(tags[i].tostring().strip())
869    except KeyError:
870        for ob in mesh['triangles']:
871            mesh['triangle_tags'].append('')
872   
873    #the outline
874    try:
875        mesh['points'] = fid.variables['points'][:]
876    except KeyError:
877        mesh['points'] = []
878    try:
879        mesh['point_attributes'] = fid.variables['point_attributes'][:]
880    except KeyError:
881        mesh['point_attributes'] = []
882        for point in mesh['points']:
883            mesh['point_attributes'].append([])
884    try:
885        mesh['outline_segments'] = fid.variables['outline_segments'][:]
886    except KeyError:
887        mesh['outline_segments'] = array([])
888    mesh['outline_segment_tags'] =[]
889    try:
890        tags = fid.variables['outline_segment_tags'][:]
891        for i, tag in enumerate(tags):
892            mesh['outline_segment_tags'].append(tags[i].tostring().strip())
893    except KeyError:
894        for ob in mesh['outline_segments']:
895            mesh['outline_segment_tags'].append('')
896    try:
897        mesh['holes'] = fid.variables['holes'][:]
898    except KeyError:
899        mesh['holes'] = array([])
900    try:
901        mesh['regions'] = fid.variables['regions'][:]
902    except KeyError:
903        mesh['regions'] = array([])
904    mesh['region_tags'] =[]
905    try:
906        tags = fid.variables['region_tags'][:]
907        for i, tag in enumerate(tags):
908            mesh['region_tags'].append(tags[i].tostring().strip())
909    except KeyError:
910        for ob in mesh['regions']:
911            mesh['region_tags'].append('')
912    try:
913        mesh['region_max_areas'] = fid.variables['region_max_areas'][:]
914    except KeyError:
915        mesh['region_max_areas'] = array([])
916    #mesh[''] = fid.variables[''][:]
917     
918    try:
919        geo_reference = Geo_reference(NetCDFObject=fid)
920        mesh['geo_reference'] = geo_reference
921    except AttributeError, e:
922        #geo_ref not compulsory
923        mesh['geo_reference'] = None
924   
925    fid.close()
926     
927    return mesh
928           
929
930## used by alpha shapes
931   
932def export_boundary_file( file_name, points, title, delimiter = ','):
933    """
934    export a file, ofile, with the format
935   
936    First line: Title variable
937    Following lines:  [point index][delimiter][point index]
938
939    file_name - the name of the new file
940    points - List of point index pairs [[p1, p2],[p3, p4]..]
941    title - info to write in the first line
942    """
943   
944    fd = open(file_name,'w')
945   
946    fd.write(title+"\n")
947    #[point index][delimiter][point index]
948    for point in points:
949        fd.write( str(point[0]) + delimiter
950                  + str(point[1]) + "\n")
951    fd.close()
952
953
954###
955#  IMPORT/EXPORT XYA FILES
956###
957
958def export_points_file(ofile,point_dict):
959    """
960    write a points file, ofile, as a binary (.msh) or text (.tsh) file
961    """
962    #this was done for all keys in the mesh file.
963    #if not mesh_dict.has_key('points'):
964    #    mesh_dict['points'] = []
965    if (ofile[-4:] == ".xya"):
966        _write_xya_file(ofile, point_dict)
967    elif (ofile[-4:] == ".pts"):
968        _write_pts_file(ofile, point_dict)
969    else:
970        msg = 'Unknown file type %s ' %ofile
971        raise IOError, msg
972               
973def import_points_file(ofile, delimiter = None, verbose = False):
974    """
975    load an .xya or .pts file
976
977    Note: will throw an IOError if it can't load the file.
978    Catch these!
979    """
980   
981    if ofile[-4:]== ".xya":
982        try:
983            if delimiter == None:
984                try:
985                    fd = open(ofile)
986                    points_dict = _read_xya_file(fd, ',')
987                except SyntaxError:
988                    fd.close()
989                    fd = open(ofile)
990                    points_dict = _read_xya_file(fd, ' ')
991            else:
992                fd = open(ofile)
993                points_dict = _read_xya_file(fd, delimiter)
994            fd.close()
995        except (IndexError,ValueError,SyntaxError):
996            fd.close()   
997            msg = 'Could not open file %s ' %ofile
998            raise IOError, msg
999        except IOError:
1000            # Catch this to add an error message
1001            msg = 'Could not open file %s ' %ofile
1002            raise IOError, msg
1003       
1004    elif ofile[-4:]== ".pts":
1005        try:
1006            points_dict = _read_pts_file(ofile, verbose)       
1007        except IOError, e:   
1008            msg = 'Could not open file %s ' %ofile
1009            raise IOError, msg
1010    else:     
1011        msg = 'Extension %s is unknown' %ofile[-4:]
1012        raise IOError, msg
1013    return points_dict
1014
1015def extent_point_atts(point_atts):
1016    """
1017    Returns 4 points representing the extent
1018    This losses attribute info.
1019    """
1020    point_atts['pointlist'] = extent(point_atts['pointlist'])
1021    point_atts['attributelist'] = {}
1022    return point_atts
1023   
1024def extent(points):
1025    points = array(points).astype(Float)
1026    max_x = min_x = points[0][0]
1027    max_y = min_y = points[0][1]
1028    for point in points[1:]:
1029        x = point[0] 
1030        if x > max_x: max_x = x
1031        if x < min_x: min_x = x           
1032        y = point[1] 
1033        if y > max_y: max_y = y
1034        if y < min_y: min_y = y
1035    extent = array([[min_x,min_y],[max_x,min_y],[max_x,max_y],[min_x,max_y]])
1036    #print "extent",extent
1037    return extent
1038   
1039def reduce_pts(infile, outfile, max_points, verbose = False):
1040    point_atts = _read_pts_file(infile)
1041    while point_atts['pointlist'].shape[0] > max_points:
1042        if verbose: print "point_atts['pointlist'].shape[0]"
1043        point_atts = half_pts(point_atts)
1044    export_points_file(outfile, point_atts)
1045 
1046def produce_half_point_files(infile, max_points, delimiter, verbose = False):
1047    point_atts = _read_pts_file(infile)
1048    root, ext = splitext(infile)
1049    outfiles = []
1050    if verbose: print "# of points", point_atts['pointlist'].shape[0]
1051    while point_atts['pointlist'].shape[0] > max_points:
1052        point_atts = half_pts(point_atts)
1053        if verbose: print "# of points", point_atts['pointlist'].shape[0]
1054        outfile = root + delimiter + str(point_atts['pointlist'].shape[0]) + ext
1055        outfiles.append(outfile)
1056        export_points_file(outfile,
1057                  point_atts)
1058    return outfiles
1059   
1060def point_atts2array(point_atts):
1061    point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
1062   
1063    for key in point_atts['attributelist'].keys():
1064        point_atts['attributelist'][key]= array(point_atts['attributelist'][key]).astype(Float)
1065    return point_atts
1066   
1067def half_pts(point_atts):
1068    point_atts2array(point_atts)
1069    point_atts['pointlist'] = point_atts['pointlist'][::2]
1070   
1071    for key in point_atts['attributelist'].keys():
1072        point_atts['attributelist'][key]=  point_atts['attributelist'][key] [::2]
1073    return point_atts
1074
1075def concatinate_attributelist(dic):
1076    """
1077    giving a dic[attribute title] = attribute
1078    return list of attribute titles, array of attributes
1079    """
1080    point_attributes = array([]).astype(Float)
1081    keys = dic.keys()
1082    key = keys.pop(0)
1083    point_attributes = reshape(dic[key],(dic[key].shape[0],1))
1084    for key in keys:
1085        #point_attributes = concatenate([point_attributes, dic[key]], axis=1)
1086        reshaped = reshape(dic[key],(dic[key].shape[0],1))
1087        point_attributes = concatenate([point_attributes, reshaped], axis=1)
1088    return dic.keys(), point_attributes
1089
1090def _read_pts_file(file_name, verbose = False):
1091    """Read .pts NetCDF file
1092   
1093    Return a dic of array of points, and dic of array of attribute
1094
1095    eg
1096    dic['points'] = [[1.0,2.0],[3.0,5.0]]
1097    dic['attributelist']['elevation'] = [[7.0,5.0]
1098    """
1099    #FIXME: (DSG)This format has issues.
1100    # There can't be an attribute called points
1101    # consider format change
1102
1103    from Scientific.IO.NetCDF import NetCDFFile
1104
1105    if verbose: print 'Reading ', file_name
1106
1107   
1108    # see if the file is there.  Throw a QUIET IO error if it isn't
1109    fd = open(file_name,'r')
1110    fd.close()
1111
1112    #throws prints to screen if file not present
1113    fid = NetCDFFile(file_name, 'r') 
1114
1115    point_atts = {} 
1116    # Get the variables
1117    point_atts['pointlist'] = array(fid.variables['points'])
1118    keys = fid.variables.keys()
1119    if verbose: print 'Got %d variables: %s' %(len(keys), keys)
1120    try:
1121        keys.remove('points')
1122    except IOError, e:       
1123        fid.close()   
1124        msg = 'Expected keyword "points" but could not find it'
1125        raise IOError, msg
1126
1127    attributes = {}
1128    for key in keys:
1129        if verbose: print "reading attribute '%s'" %key
1130       
1131        attributes[key] = array(fid.variables[key])
1132       
1133    point_atts['attributelist'] = attributes
1134   
1135    try:
1136        geo_reference = Geo_reference(NetCDFObject=fid)
1137        point_atts['geo_reference'] = geo_reference
1138    except AttributeError, e:
1139        #geo_ref not compulsory
1140        point_atts['geo_reference'] = None
1141   
1142   
1143    fid.close()
1144   
1145    #print "point_atts",point_atts
1146    return point_atts
1147
1148def _write_pts_file(file_name, point_atts):
1149    """Write .pts NetCDF file   
1150
1151    WARNING: This function mangles the point_atts data structure
1152    """
1153    #FIXME: (DSG)This format has issues.
1154    # There can't be an attribute called points
1155    # consider format change
1156   
1157    from Scientific.IO.NetCDF import NetCDFFile
1158    point_atts2array(point_atts)
1159    # NetCDF file definition
1160    outfile = NetCDFFile(file_name, 'w')
1161   
1162    #Create new file
1163    outfile.institution = 'Geoscience Australia'
1164    outfile.description = 'NetCDF format for compact and portable storage ' +\
1165                      'of spatial point data'
1166   
1167    # dimension definitions
1168    shape = point_atts['pointlist'].shape[0]
1169    outfile.createDimension('number_of_points', shape) 
1170    outfile.createDimension('number_of_dimensions', 2) #This is 2d data
1171   
1172    # variable definition
1173    outfile.createVariable('points', Float, ('number_of_points',
1174                                             'number_of_dimensions'))
1175
1176    #create variables 
1177    outfile.variables['points'][:] = point_atts['pointlist'] #.astype(Float32)
1178    for key in point_atts['attributelist'].keys():
1179        outfile.createVariable(key, Float, ('number_of_points',))
1180        outfile.variables[key][:] = point_atts['attributelist'][key] #.astype(Float32)
1181       
1182    if point_atts.has_key('geo_reference') and not point_atts['geo_reference'] == None:
1183        point_atts['geo_reference'].write_NetCDF(outfile)
1184       
1185    outfile.close() 
1186   
1187def _read_xya_file(fd, delimiter):
1188    #lines = fd.readlines()
1189    points = []
1190    pointattributes = []
1191    #if len(lines) <= 1:
1192    #    raise SyntaxError
1193    title = fd.readline()
1194    #title = lines.pop(0) # the first (title) line
1195    att_names = clean_line(title,delimiter)
1196
1197    att_dict = {}
1198    line = fd.readline()
1199    numbers = clean_line(line,delimiter)
1200    #for line in lines:
1201    while len(numbers) > 1:
1202        #print "line >%s" %line
1203        #numbers = clean_line(line,delimiter)
1204        #print "numbers >%s<" %numbers
1205        #if len(numbers) < 2 and numbers != []:
1206           
1207            # A line without two numbers
1208            # or a bad delimiter
1209            #FIXME dsg-dsg change error that is raised.
1210            #raise SyntaxError
1211        if numbers != []:
1212            try:
1213                x = float(numbers[0])
1214                y = float(numbers[1])
1215                points.append([x,y])
1216                numbers.pop(0)
1217                numbers.pop(0)
1218                #attributes = []
1219                #print "att_names",att_names
1220                #print "numbers",numbers
1221                if len(att_names) != len(numbers):
1222                    fd.close()
1223                    # It might not be a problem with the title
1224                    #raise TitleAmountError
1225                    raise IOError
1226                for i,num in enumerate(numbers):
1227                    num.strip()
1228                    if num != '\n' and num != '':
1229                        #attributes.append(float(num))
1230                        att_dict.setdefault(att_names[i],[]).append(float(num))
1231                   
1232            except ValueError:
1233                raise SyntaxError
1234        line = fd.readline()
1235        numbers = clean_line(line,delimiter)
1236       
1237    if line == '':
1238        # end of file
1239        geo_reference = None
1240    else:
1241        geo_reference = Geo_reference(ASCIIFile=fd,read_title=line)
1242   
1243    xya_dict = {}
1244    xya_dict['pointlist'] = array(points).astype(Float)
1245   
1246    for key in att_dict.keys():
1247        att_dict[key] = array(att_dict[key]).astype(Float)
1248    xya_dict['attributelist'] = att_dict
1249    xya_dict['geo_reference'] = geo_reference
1250    #print "xya_dict",xya_dict
1251    return xya_dict
1252
1253#FIXME(dsg), turn this dict plus methods into a class?
1254def take_points(dict,indices_to_keep):
1255    dict = point_atts2array(dict)
1256    #FIXME maybe the points data structure should become a class?
1257    dict['pointlist'] = take(dict['pointlist'],indices_to_keep)
1258
1259    for key in dict['attributelist'].keys():
1260        dict['attributelist'][key]= take(dict['attributelist'][key],
1261                                         indices_to_keep)
1262    return dict
1263   
1264def add_point_dictionaries (dict1, dict2):
1265    """
1266    """
1267    dict1 = point_atts2array(dict1)
1268    dict2 = point_atts2array(dict2)
1269   
1270    combined = {} 
1271    combined['pointlist'] = concatenate((dict2['pointlist'],
1272                                         dict1['pointlist']),axis=0)   
1273    atts = {}   
1274    for key in dict2['attributelist'].keys():
1275        atts[key]= concatenate((dict2['attributelist'][key],
1276                                dict1['attributelist'][key]), axis=0)
1277    combined['attributelist']=atts
1278    combined['geo_reference'] = dict1['geo_reference']
1279    return combined
1280                 
1281def _write_xya_file( file_name, xya_dict, delimiter = ','):
1282    """
1283    export a file, ofile, with the xya format
1284   
1285    """
1286    points = xya_dict['pointlist'] 
1287    pointattributes = xya_dict['attributelist']
1288   
1289    fd = open(file_name,'w')
1290 
1291    titlelist = ""
1292    for title in pointattributes.keys():
1293        titlelist = titlelist + title + delimiter
1294    titlelist = titlelist[0:-len(delimiter)] # remove the last delimiter
1295    fd.write(titlelist+"\n")
1296    #<vertex #> <x> <y> [attributes]
1297    for i,vert in enumerate( points):
1298       
1299        attlist = ","
1300        for att in pointattributes.keys():
1301            attlist = attlist + str(pointattributes[att][i])+ delimiter
1302        attlist = attlist[0:-len(delimiter)] # remove the last delimiter
1303        attlist.strip()
1304        fd.write( str(vert[0]) + delimiter
1305                  + str(vert[1])
1306                  + attlist + "\n")
1307   
1308    # geo_reference info
1309    if xya_dict.has_key('geo_reference') and \
1310           not xya_dict['geo_reference'] is None:
1311        xya_dict['geo_reference'].write_ASCII(fd)
1312    fd.close()
1313     
1314if __name__ == "__main__":
1315    pass
Note: See TracBrowser for help on using the repository browser.