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

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

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