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

Last change on this file since 4901 was 4901, checked in by duncan, 16 years ago

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