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

Last change on this file since 5079 was 4903, checked in by duncan, 17 years ago

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