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

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

generalising so info without attributes can be used
File size: 39.7 KB
Line 
1"""
2The format for a .xya file is:
31st line:     [attribute names]
4other lines:  x y [attributes]
5
6for example:
7elevation, friction
80.6, 0.7, 4.9, 0.3
91.9, 2.8, 5, 0.3
102.7, 2.4, 5.2, 0.3
11
12The first two columns are always implicitly assumed to be x, y coordinates.
13Use the same delimiter for the attribute names and the data
14
15An xya file can optionally end with
16#geo reference
1756
18466600.0
198644444.0
20
21When the 1st # is the zone,
222nd # the xllcorner and
233rd # the yllcorner
24The format for a Points dictionary is:
25
26  ['pointlist'] a 2 column array describing points. 1st column x, 2nd column y.
27  ['attributelist'], a dictionary of 1D arrays, representing attribute values
28  at the point.  The dictionary key is the attribute header.
29  ['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       
492        if vertices_attributes == []:
493            attlist = ""
494        else:
495            for att in vertices_attributes[index]:
496                attlist = attlist + str(att)+" "
497        attlist.strip()
498        fd.write(str(index) + " "
499                 + str(vert[0]) + " "
500                 + str(vert[1]) + " "
501                 + attlist + "\n")
502        index += 1
503
504    # write comments for title
505    fd.write("# attribute column titles ...Triangulation Vertex Titles..." + "\n")
506    for title in vertices_attribute_titles:
507        fd.write(title + "\n")
508       
509    #<# of triangles>
510    n = len(triangles)
511    fd.write(str(n) + " # <# of triangles>, next lines <triangle #> [<vertex #>] [<neigbouring triangle #>] [attribute of region] ...Triangulation Triangles..." + "\n")
512       
513    # <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
514    for index in range(n):
515        neighbors = ""
516        tri = triangles[index]
517        if triangle_neighbors == []:
518            neighbors = "-1 -1 -1 "
519        else:
520            for neighbor in triangle_neighbors[index]:
521                if neighbor:
522                    neighbors += str(neighbor) + " "
523                else:
524                    if neighbor == 0:
525                        neighbors +=  "0 "
526                    else:
527                        neighbors +=  "-1 "
528        #Warning even though a list is past, only the first value
529        #is written.  There's an assumption that the list only
530        # contains one item. This assumption is made since the
531        #dict that's being passed around is also be used to communicate
532        # with triangle, and it seems to have the option of returning
533        # more than one value for triangle attributex
534        if triangles_attributes == [] or triangles_attributes[index] == ['']:
535            att = ""
536        else:
537            att = str(triangles_attributes[index])
538        fd.write(str(index) + " "
539                 + str(tri[0]) + " " 
540                 + str(tri[1]) + " " 
541                 + str(tri[2]) + " " 
542                 + neighbors + " "
543                 + att + "\n")
544           
545    #One line:  <# of segments>
546    fd.write(str(len(segments)) + 
547             " # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Triangulation Segments..." + "\n")
548       
549    #Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
550    for i in range(len(segments)):
551        seg = segments[i]
552        fd.write(str(i) + " "
553                 + str(seg[0]) + " " 
554                 + str(seg[1]) + " " 
555                 + str(segment_tags[i]) + "\n")
556
557
558def _write_tsh_file(ofile,mesh_dict):
559            fd = open(ofile,'w')
560            _write_ASCII_triangulation(fd,mesh_dict)
561            _write_ASCII_outline(fd,mesh_dict)   
562            fd.close()
563
564def _write_ASCII_outline(fd,
565                     dict):
566    points = dict['points']
567    point_attributes = dict['point_attributes']
568    segments = dict['outline_segments']
569    segment_tags = dict['outline_segment_tags']
570    holes = dict['holes']
571    regions = dict['regions']
572    region_tags = dict['region_tags']
573    region_max_areas = dict['region_max_areas']
574       
575    num_points = str(len(points))
576    if (num_points == "0"):
577        num_point_atts = "0"
578    else:
579        num_point_atts = str(len(point_attributes[0]))
580    fd.write(num_points + " " + num_point_atts + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Mesh Vertices..." + "\n")
581       
582    # <x> <y> [attributes]
583    for i,point in enumerate(points):
584        attlist = ""
585        for att in point_attributes[i]:
586            attlist = attlist + str(att)+" "
587        attlist.strip()
588        fd.write(str(i) + " "
589                 +str(point[0]) + " "
590                 + str(point[1]) + " "
591                 + attlist + "\n")
592           
593    #One line:  <# of segments>
594    fd.write(str(len(segments)) + 
595             " # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Mesh Segments..." + "\n")
596       
597    #Following lines:  <vertex #>  <vertex #> [boundary tag]
598    for i,seg in enumerate(segments):
599        fd.write(str(i) + " "
600                 + str(seg[0]) + " " 
601                 + str(seg[1]) + " " 
602                 + str(segment_tags[i]) + "\n")
603
604    #One line:  <# of holes>
605    fd.write(str(len(holes)) + 
606             " # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes..." + "\n")   
607    # <x> <y>
608    for i,h in enumerate(holes):
609        fd.write(str(i) + " "
610                 + str(h[0]) + " "
611                 + str(h[1]) + "\n")
612       
613    #One line:  <# of regions>
614    fd.write(str(len(regions)) + 
615             " # <# of regions>, next lines <Region #> <x> <y> <tag>...Mesh Regions..." + "\n")   
616    # <index> <x> <y> <tag>
617    #print "regions",regions
618    for i,r in enumerate(regions):
619        #print "r",r
620        fd.write(str(i) + " " 
621                 + str(r[0]) + " "
622                 + str(r[1])+ " "
623                 + str(region_tags[i]) + "\n")
624       
625    # <index> [<MaxArea>|""]
626       
627    #One line:  <# of regions>
628    fd.write(str(len(regions)) + 
629             " # <# of regions>, next lines <Region #> [Max Area] ...Mesh Regions..." + "\n") 
630    for i,r in enumerate(regions):
631        area = str(region_max_areas[i])
632           
633        fd.write(str(i) + " " + area + "\n")
634
635    # geo_reference info
636    if dict.has_key('geo_reference') and not dict['geo_reference'] is None:
637        dict['geo_reference'].write_ASCII(fd)
638
639def _write_msh_file(file_name, mesh):
640    """Write .msh NetCDF file   
641
642    WARNING: This function mangles the mesh data structure
643    """
644 
645    # FIXME(Ole and John): We ran into a problem on Bogong (64 bit) where integers appeared as arrays.
646    # This may be similar to problem seen by Steve in changeset:2778 where he had to wrap them in int.
647    # Now we are trying with the native Integer format (Int == 'l' == Int64). However, that caused casting errors, when 64bit arrays
648    # are to be assigned to their NetCDF counterparts. It seems that the NetCDF arrays are 32bit even though they are created with
649    # the type Int64. Need to look at the NetCDF library in more detail.
650   
651    IntType = Int32
652    #IntType = Int
653   
654    #
655    #the triangulation
656    mesh['vertices'] = array(mesh['vertices']).astype(Float)
657    mesh['vertex_attributes'] = array(mesh['vertex_attributes']).astype(Float)
658    mesh['vertex_attribute_titles'] = array(mesh['vertex_attribute_titles']).astype(Character) 
659    mesh['segments'] = array(mesh['segments']).astype(IntType)
660    mesh['segment_tags'] = array(mesh['segment_tags']).astype(Character)
661    mesh['triangles'] = array(mesh['triangles']).astype(IntType)
662    mesh['triangle_tags'] = array(mesh['triangle_tags']) #.astype(Character)
663    mesh['triangle_neighbors'] = array(mesh['triangle_neighbors']).astype(IntType) 
664   
665    #the outline
666    mesh['points'] = array(mesh['points']).astype(Float)
667    mesh['point_attributes'] = array(mesh['point_attributes']).astype(Float)
668    mesh['outline_segments'] = array(mesh['outline_segments']).astype(IntType)
669    mesh['outline_segment_tags'] = array(mesh['outline_segment_tags']).astype(Character)
670    mesh['holes'] = array(mesh['holes']).astype(Float)
671    mesh['regions'] = array(mesh['regions']).astype(Float)
672    mesh['region_tags'] = array(mesh['region_tags']).astype(Character)
673    mesh['region_max_areas'] = array(mesh['region_max_areas']).astype(Float)
674   
675    #mesh = mesh_dict2array(mesh)
676    #print "new_mesh",new_mesh
677    #print "mesh",mesh
678   
679    # NetCDF file definition
680    try:
681        outfile = NetCDFFile(file_name, 'w')
682    except IOError:
683        msg = 'File %s could not be created' %file_name
684        raise msg
685       
686    #Create new file
687    outfile.institution = 'Geoscience Australia'
688    outfile.description = 'NetCDF format for compact and portable storage ' +\
689                      'of spatial point data'
690   
691    # dimension definitions - fixed
692    outfile.createDimension('num_of_dimensions', 2) #This is 2d data
693    outfile.createDimension('num_of_segment_ends', 2) #Segs have two points
694    outfile.createDimension('num_of_triangle_vertices', 3)
695    outfile.createDimension('num_of_triangle_faces', 3)
696    outfile.createDimension('num_of_region_max_area', 1)
697
698    # Create dimensions, variables and set the variables
699   
700    # trianglulation
701    # vertices
702    if (mesh['vertices'].shape[0] > 0):
703        outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
704        outfile.createVariable('vertices', Float, ('num_of_vertices',
705                                                   'num_of_dimensions'))
706        outfile.variables['vertices'][:] = mesh['vertices']
707        if (mesh['vertex_attributes'].shape[0] > 0 and mesh['vertex_attributes'].shape[1] > 0):
708            outfile.createDimension('num_of_vertex_attributes',
709                                    mesh['vertex_attributes'].shape[1])
710            outfile.createDimension('num_of_vertex_attribute_title_chars',
711                                    mesh['vertex_attribute_titles'].shape[1])
712            outfile.createVariable('vertex_attributes',
713                                   Float,
714                                   ('num_of_vertices',
715                                    'num_of_vertex_attributes'))   
716            outfile.createVariable('vertex_attribute_titles',
717                                   Character,
718                                   ( 'num_of_vertex_attributes',
719                                     'num_of_vertex_attribute_title_chars' ))
720            outfile.variables['vertex_attributes'][:] = \
721                                                      mesh['vertex_attributes']
722            outfile.variables['vertex_attribute_titles'][:] = \
723                                     mesh['vertex_attribute_titles']
724    # segments
725    if (mesh['segments'].shape[0] > 0):       
726        outfile.createDimension('num_of_segments',
727                                mesh['segments'].shape[0])
728        outfile.createVariable('segments',
729                               IntType,
730                               ('num_of_segments', 'num_of_segment_ends'))
731        outfile.variables['segments'][:] = mesh['segments']
732        if (mesh['segment_tags'].shape[1] > 0):
733            outfile.createDimension('num_of_segment_tag_chars',
734                                    mesh['segment_tags'].shape[1])
735            outfile.createVariable('segment_tags',
736                                   Character,
737                                   ('num_of_segments',
738                                    'num_of_segment_tag_chars'))
739            outfile.variables['segment_tags'][:] = mesh['segment_tags']
740    # triangles   
741    if (mesh['triangles'].shape[0] > 0):
742        outfile.createDimension('num_of_triangles',
743                                mesh['triangles'].shape[0])
744        outfile.createVariable('triangles',
745                               IntType,
746                               ('num_of_triangles',
747                                'num_of_triangle_vertices'))
748        outfile.createVariable('triangle_neighbors',
749                               IntType,
750                               ('num_of_triangles',
751                                'num_of_triangle_faces'))
752        outfile.variables['triangles'][:] = mesh['triangles']
753        outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
754        if (mesh['triangle_tags'].shape[1] > 0):
755            outfile.createDimension('num_of_triangle_tag_chars',
756                                    mesh['triangle_tags'].shape[1]) 
757            outfile.createVariable('triangle_tags',
758                                   Character,
759                                   ('num_of_triangles',
760                                    'num_of_triangle_tag_chars'))
761            outfile.variables['triangle_tags'][:] = mesh['triangle_tags']
762   
763
764    # outline
765    # points
766    if (mesh['points'].shape[0] > 0):
767        outfile.createDimension('num_of_points', mesh['points'].shape[0])
768        outfile.createVariable('points', Float, ('num_of_points',
769                                                 'num_of_dimensions'))
770        outfile.variables['points'][:] = mesh['points'] 
771        if (mesh['point_attributes'].shape[0] > 0  and mesh['point_attributes'].shape[1] > 0):
772            outfile.createDimension('num_of_point_attributes',
773                                    mesh['point_attributes'].shape[1])
774            outfile.createVariable('point_attributes',
775                                   Float,
776                                   ('num_of_points',
777                                    'num_of_point_attributes'))
778            outfile.variables['point_attributes'][:] = mesh['point_attributes']
779    # outline_segments
780    if (mesh['outline_segments'].shape[0] > 0):
781        outfile.createDimension('num_of_outline_segments',
782                                mesh['outline_segments'].shape[0])
783        outfile.createVariable('outline_segments',
784                               IntType,
785                               ('num_of_outline_segments',
786                                'num_of_segment_ends'))
787        outfile.variables['outline_segments'][:] = mesh['outline_segments']
788        if (mesh['outline_segment_tags'].shape[1] > 0):
789            outfile.createDimension('num_of_outline_segment_tag_chars',
790                                    mesh['outline_segment_tags'].shape[1])
791            outfile.createVariable('outline_segment_tags',
792                                   Character,
793                                   ('num_of_outline_segments',
794                                    'num_of_outline_segment_tag_chars'))
795            outfile.variables['outline_segment_tags'][:] = mesh['outline_segment_tags']
796    # holes
797    if (mesh['holes'].shape[0] > 0):
798        outfile.createDimension('num_of_holes', mesh['holes'].shape[0])
799        outfile.createVariable('holes', Float, ('num_of_holes',
800                                                'num_of_dimensions'))
801        outfile.variables['holes'][:] = mesh['holes']
802    # regions
803    if (mesh['regions'].shape[0] > 0):
804        outfile.createDimension('num_of_regions', mesh['regions'].shape[0])
805        outfile.createVariable('regions', Float, ('num_of_regions',
806                                                  'num_of_dimensions'))
807        outfile.createVariable('region_max_areas',
808                               Float,
809                               ('num_of_regions',))
810        outfile.variables['regions'][:] = mesh['regions']
811        outfile.variables['region_max_areas'][:] = mesh['region_max_areas']
812        if (mesh['region_tags'].shape[1] > 0):
813            outfile.createDimension('num_of_region_tag_chars',
814                                    mesh['region_tags'].shape[1])
815            outfile.createVariable('region_tags',
816                                   Character,
817                                   ('num_of_regions',
818                                    'num_of_region_tag_chars'))
819            outfile.variables['region_tags'][:] = mesh['region_tags']
820
821    # geo_reference info
822    if mesh.has_key('geo_reference') and not mesh['geo_reference'] == None:
823        mesh['geo_reference'].write_NetCDF(outfile)
824                                                 
825    outfile.close()
826
827
828   
829def _read_msh_file(file_name):
830    """
831    Read in an msh file.
832
833    """
834               
835           
836    #Check contents
837    #Get NetCDF
838   
839    # see if the file is there.  Throw a QUIET IO error if it isn't
840    fd = open(file_name,'r')
841    fd.close()
842
843    #throws prints to screen if file not present
844    fid = NetCDFFile(file_name, 'r') 
845
846    mesh = {}
847    # Get the variables
848    # the triangulation
849    try:
850        mesh['vertices'] = fid.variables['vertices'][:]
851    except KeyError:
852        mesh['vertices'] = array([])
853    try:
854        mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:]
855    except KeyError:
856        mesh['vertex_attributes'] = []
857        for ob in mesh['vertices']:
858            mesh['vertex_attributes'].append([])
859    mesh['vertex_attribute_titles'] = []
860    try:
861        titles = fid.variables['vertex_attribute_titles'][:]
862        for i, title in enumerate(titles):
863            mesh['vertex_attribute_titles'].append(titles[i].tostring().strip())
864    except KeyError:
865        pass 
866    try:
867        mesh['segments'] = fid.variables['segments'][:]
868    except KeyError:
869        mesh['segments'] = array([])
870    mesh['segment_tags'] =[]
871    try:
872        tags = fid.variables['segment_tags'][:]
873        for i, tag in enumerate(tags):
874            mesh['segment_tags'].append(tags[i].tostring().strip())
875    except KeyError:
876        for ob in mesh['segments']:
877            mesh['segment_tags'].append('')
878    try:
879        mesh['triangles'] = fid.variables['triangles'][:]
880        mesh['triangle_neighbors'] = fid.variables['triangle_neighbors'][:]
881    except KeyError:
882        mesh['triangles'] = array([])
883        mesh['triangle_neighbors'] = array([])
884    mesh['triangle_tags'] =[]
885    try:
886        tags = fid.variables['triangle_tags'][:]
887        for i, tag in enumerate(tags):
888            mesh['triangle_tags'].append(tags[i].tostring().strip())
889    except KeyError:
890        for ob in mesh['triangles']:
891            mesh['triangle_tags'].append('')
892   
893    #the outline
894    try:
895        mesh['points'] = fid.variables['points'][:]
896    except KeyError:
897        mesh['points'] = []
898    try:
899        mesh['point_attributes'] = fid.variables['point_attributes'][:]
900    except KeyError:
901        mesh['point_attributes'] = []
902        for point in mesh['points']:
903            mesh['point_attributes'].append([])
904    try:
905        mesh['outline_segments'] = fid.variables['outline_segments'][:]
906    except KeyError:
907        mesh['outline_segments'] = array([])
908    mesh['outline_segment_tags'] =[]
909    try:
910        tags = fid.variables['outline_segment_tags'][:]
911        for i, tag in enumerate(tags):
912            mesh['outline_segment_tags'].append(tags[i].tostring().strip())
913    except KeyError:
914        for ob in mesh['outline_segments']:
915            mesh['outline_segment_tags'].append('')
916    try:
917        mesh['holes'] = fid.variables['holes'][:]
918    except KeyError:
919        mesh['holes'] = array([])
920    try:
921        mesh['regions'] = fid.variables['regions'][:]
922    except KeyError:
923        mesh['regions'] = array([])
924    mesh['region_tags'] =[]
925    try:
926        tags = fid.variables['region_tags'][:]
927        for i, tag in enumerate(tags):
928            mesh['region_tags'].append(tags[i].tostring().strip())
929    except KeyError:
930        for ob in mesh['regions']:
931            mesh['region_tags'].append('')
932    try:
933        mesh['region_max_areas'] = fid.variables['region_max_areas'][:]
934    except KeyError:
935        mesh['region_max_areas'] = array([])
936    #mesh[''] = fid.variables[''][:]
937     
938    try:
939        geo_reference = Geo_reference(NetCDFObject=fid)
940        mesh['geo_reference'] = geo_reference
941    except AttributeError, e:
942        #geo_ref not compulsory
943        mesh['geo_reference'] = None
944   
945    fid.close()
946     
947    return mesh
948           
949
950## used by alpha shapes
951   
952def export_boundary_file( file_name, points, title, delimiter = ','):
953    """
954    export a file, ofile, with the format
955   
956    First line: Title variable
957    Following lines:  [point index][delimiter][point index]
958
959    file_name - the name of the new file
960    points - List of point index pairs [[p1, p2],[p3, p4]..]
961    title - info to write in the first line
962    """
963   
964    fd = open(file_name,'w')
965   
966    fd.write(title+"\n")
967    #[point index][delimiter][point index]
968    for point in points:
969        fd.write( str(point[0]) + delimiter
970                  + str(point[1]) + "\n")
971    fd.close()
972
973
974###
975#  IMPORT/EXPORT POINTS FILES
976###
977
978def extent_point_atts(point_atts):
979    """
980    Returns 4 points representing the extent
981    This losses attribute info.
982    """
983    point_atts['pointlist'] = extent(point_atts['pointlist'])
984    point_atts['attributelist'] = {}
985    return point_atts
986   
987def extent(points):
988    points = array(points).astype(Float)
989    max_x = min_x = points[0][0]
990    max_y = min_y = points[0][1]
991    for point in points[1:]:
992        x = point[0] 
993        if x > max_x: max_x = x
994        if x < min_x: min_x = x           
995        y = point[1] 
996        if y > max_y: max_y = y
997        if y < min_y: min_y = y
998    extent = array([[min_x,min_y],[max_x,min_y],[max_x,max_y],[min_x,max_y]])
999    #print "extent",extent
1000    return extent
1001   
1002def reduce_pts(infile, outfile, max_points, verbose = False):
1003    """
1004    reduces a points file by removong every second point until the # of points
1005    is less than max_points.
1006    """
1007    # check out pts2rectangular in least squares, and the use of reduction.
1008    # Maybe it does the same sort of thing?
1009    point_atts = _read_pts_file(infile)
1010    while point_atts['pointlist'].shape[0] > max_points:
1011        if verbose: print "point_atts['pointlist'].shape[0]"
1012        point_atts = half_pts(point_atts)
1013    export_points_file(outfile, point_atts)
1014 
1015def produce_half_point_files(infile, max_points, delimiter, verbose = False):
1016    point_atts = _read_pts_file(infile)
1017    root, ext = splitext(infile)
1018    outfiles = []
1019    if verbose: print "# of points", point_atts['pointlist'].shape[0]
1020    while point_atts['pointlist'].shape[0] > max_points:
1021        point_atts = half_pts(point_atts)
1022        if verbose: print "# of points", point_atts['pointlist'].shape[0]
1023        outfile = root + delimiter + str(point_atts['pointlist'].shape[0]) + ext
1024        outfiles.append(outfile)
1025        export_points_file(outfile,
1026                  point_atts)
1027    return outfiles
1028   
1029def point_atts2array(point_atts):
1030    point_atts['pointlist'] = array(point_atts['pointlist']).astype(Float)
1031   
1032    for key in point_atts['attributelist'].keys():
1033        point_atts['attributelist'][key]= array(point_atts['attributelist'][key]).astype(Float)
1034    return point_atts
1035   
1036def half_pts(point_atts):
1037    point_atts2array(point_atts)
1038    point_atts['pointlist'] = point_atts['pointlist'][::2]
1039   
1040    for key in point_atts['attributelist'].keys():
1041        point_atts['attributelist'][key]=  point_atts['attributelist'][key] [::2]
1042    return point_atts
1043
1044def concatinate_attributelist(dic):
1045    """
1046    giving a dic[attribute title] = attribute
1047    return list of attribute titles, array of attributes
1048    """
1049    point_attributes = array([]).astype(Float)
1050    keys = dic.keys()
1051    key = keys.pop(0)
1052    point_attributes = reshape(dic[key],(dic[key].shape[0],1))
1053    for key in keys:
1054        #point_attributes = concatenate([point_attributes, dic[key]], axis=1)
1055        reshaped = reshape(dic[key],(dic[key].shape[0],1))
1056        point_attributes = concatenate([point_attributes, reshaped], axis=1)
1057    return dic.keys(), point_attributes
1058
1059#FIXME(dsg), turn this dict plus methods into a class?
1060def take_points(dict,indices_to_keep):
1061    dict = point_atts2array(dict)
1062    #FIXME maybe the points data structure should become a class?
1063    dict['pointlist'] = take(dict['pointlist'],indices_to_keep)
1064
1065    for key in dict['attributelist'].keys():
1066        dict['attributelist'][key]= take(dict['attributelist'][key],
1067                                         indices_to_keep)
1068    return dict
1069   
1070def add_point_dictionaries (dict1, dict2):
1071    """
1072    """
1073    dict1 = point_atts2array(dict1)
1074    dict2 = point_atts2array(dict2)
1075   
1076    combined = {} 
1077    combined['pointlist'] = concatenate((dict2['pointlist'],
1078                                         dict1['pointlist']),axis=0)   
1079    atts = {}   
1080    for key in dict2['attributelist'].keys():
1081        atts[key]= concatenate((dict2['attributelist'][key],
1082                                dict1['attributelist'][key]), axis=0)
1083    combined['attributelist']=atts
1084    combined['geo_reference'] = dict1['geo_reference']
1085    return combined
1086     
1087if __name__ == "__main__":
1088    pass
Note: See TracBrowser for help on using the repository browser.