source: branches/source_numpy_conversion/anuga/load_mesh/loadASCII.py @ 6768

Last change on this file since 6768 was 5946, checked in by rwilson, 16 years ago

More NumPy? changes.
File size: 41.2 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
59import numpy
60from os.path import splitext
61
62from anuga.config import Float, Character, Int
63from anuga.coordinate_transforms.geo_reference import Geo_reference,TITLE, TitleError
64
65from Scientific.IO.NetCDF import NetCDFFile
66   
67import exceptions
68class TitleAmountError(exceptions.Exception): pass
69
70NOMAXAREA=-999
71
72# This is used by pmesh
73def import_mesh_file(ofile):
74    """
75    read a mesh file, either .tsh or .msh
76   
77    Note: will throw an IOError if it can't load the file.
78    Catch these!
79    """
80    try:
81        if ofile[-4:]== ".tsh":
82            dict = _read_tsh_file(ofile)
83        elif ofile[-4:]== ".msh":
84            dict = _read_msh_file(ofile)
85        else:     
86            msg = 'Extension .%s is unknown' %ofile[-4:]
87            raise IOError, msg
88    except (TitleError,SyntaxError,IndexError, ValueError): #FIXME No test for ValueError
89            msg = 'File could not be opened'
90            raise IOError, msg 
91    return dict
92
93
94def export_mesh_file(ofile,mesh_dict):
95    """
96    write a file, ofile, with the format
97   
98    First line:  <# of vertices> <# of attributes>
99    Following lines:  <vertex #> <x> <y> [attributes]
100    One line:  <# of triangles>
101    Following lines:  <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
102    One line:  <# of segments>
103    Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
104    """
105    #FIXME DSG-anyone: automate
106    if not mesh_dict.has_key('points'):
107        mesh_dict['points'] = []
108    if not mesh_dict.has_key('point_attributes'):
109        mesh_dict['point_attributes'] = []
110    if not mesh_dict.has_key('outline_segments'):
111        mesh_dict['outline_segments'] = []
112    if not mesh_dict.has_key('outline_segment_tags'):
113        mesh_dict['outline_segment_tags'] = []
114    if not mesh_dict.has_key('holes'):
115        mesh_dict['holes'] = []
116    if not mesh_dict.has_key('regions'):
117        mesh_dict['regions'] = []
118       
119    if not mesh_dict.has_key('region_tags'):
120        mesh_dict['region_tags'] = []
121    if not mesh_dict.has_key('region_max_areas'):
122        mesh_dict['region_max_areas'] = []
123    if not mesh_dict.has_key('vertices'):
124        mesh_dict['vertices'] = []
125    if not mesh_dict.has_key('vertex_attributes'):
126        mesh_dict['vertex_attributes'] = []
127    if not mesh_dict.has_key('vertex_attribute_titles'):
128        mesh_dict['vertex_attribute_titles'] = []
129    if not mesh_dict.has_key('segments'):
130        mesh_dict['segments'] = []
131    if not mesh_dict.has_key('segment_tags'):
132        mesh_dict['segment_tags'] = []
133    if not mesh_dict.has_key('triangles'):
134        mesh_dict['triangles'] = []
135    if not mesh_dict.has_key('triangle_tags'):
136        mesh_dict['triangle_tags'] = []
137    if not mesh_dict.has_key('triangle_neighbors'):
138        mesh_dict['triangle_neighbors'] = []
139    #print "DSG************"
140    #print "mesh_dict",mesh_dict
141    #print "DSG************"
142   
143    if (ofile[-4:] == ".tsh"):
144        _write_tsh_file(ofile,mesh_dict)
145    elif (ofile[-4:] == ".msh"):
146        _write_msh_file(ofile, mesh_dict)
147    else:
148        msg = 'Unknown file type %s ' %ofile
149        raise IOError, msg
150
151def _read_tsh_file(ofile):
152    """
153    Read the text file format for meshes
154    """
155    fd = open(ofile,'r')
156    dict = _read_triangulation(fd)
157    dict_mesh = _read_outline(fd)
158    for element in dict_mesh.keys():
159        dict[element] = dict_mesh[element]
160    fd.close()
161    return dict
162
163   
164def _read_triangulation(fd):
165    """
166    Read the generated triangulation, NOT the outline.
167    """
168    delimiter = " "
169    ######### loading the point info
170    line = fd.readline()
171    #print line
172    fragments = line.split()
173    #for fragment in fragments:
174    #print fragment
175    if fragments ==[]:
176        NumOfVertices = 0
177        NumOfVertAttributes = 0
178    else:
179        NumOfVertices = fragments[0]
180        NumOfVertAttributes = fragments[1]
181    points = []
182    pointattributes = []
183    for index in range(int(NumOfVertices)):       
184        #print index
185        fragments = fd.readline().split()
186        #print fragments
187        fragments.pop(0) #pop off the index
188       
189        # pop the x y off so we're left with a list of attributes       
190        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
191        points.append(vert)
192        apointattributes  = []
193        #print fragments
194        for fragment in fragments:
195            apointattributes.append(float(fragment))
196        if apointattributes != []:
197            pointattributes.append(apointattributes)
198       
199    ######### loading the point title info
200    line = fd.readline()
201    #print "point title comments",line
202    vertTitle = []
203    for index in range(int(NumOfVertAttributes)):       
204        #print index
205        fragments = fd.readline().strip() 
206        vertTitle.append(fragments)
207        #print vertTitle
208       
209    ######### loading the triangle info
210    line = fd.readline()
211    #print "triangle comments",line
212    fragments = line.split()
213    #for fragment in fragments:
214    #    print fragment
215    NumOfTriangles = fragments[0]
216    triangles = []
217    triangleattributes = []
218    triangleneighbors = []
219    for index in range(int(NumOfTriangles)):       
220        #print index
221        line = fd.readline()
222        line.strip() # so we can get the region string
223        fragments = line.split()
224        #print "triangle info", fragments
225        fragments.pop(0) #pop off the index
226       
227        tri = [int(fragments[0]),int(fragments[1]),int(fragments[2])]
228        triangles.append(tri)
229        neighbors = [int(fragments[3]),int(fragments[4]),int(fragments[5])]
230        triangleneighbors.append(neighbors)
231        for x in range(7): # remove index [<vertex #>] [<neigbouring tri #>]
232            line = line[find(line,delimiter):] # remove index
233            line = line.lstrip()
234        stringtag = line.strip()
235        triangleattributes.append(stringtag)
236       
237    ######### loading the segment info
238    line = fd.readline()
239    #print "seg comment line",line
240    fragments = line.split()
241    #for fragment in fragments:
242    #    print fragment
243    NumOfSegments = fragments[0]
244    segments = []
245    segmenttags = []
246    for index in range(int(NumOfSegments)):       
247        #print index
248        line = fd.readline()
249        line.strip() # to get the segment string
250        fragments = line.split()
251        #print fragments
252        fragments.pop(0) #pop off the index
253        seg = [int(fragments[0]),int(fragments[1])]
254        segments.append(seg)
255        line = line[find(line,delimiter):] # remove index
256        line = line.lstrip()
257        line = line[find(line,delimiter):] # remove x
258        line = line.lstrip()
259        line = line[find(line,delimiter):] # remove y
260        stringtag = line.strip()
261        segmenttags.append(stringtag)
262
263           
264    meshDict = {}
265    meshDict['vertices'] = points
266    if pointattributes == []:
267        meshDict['vertex_attributes'] = None
268    else:
269        meshDict['vertex_attributes'] = pointattributes
270    meshDict['triangles'] = triangles
271    meshDict['triangle_tags'] = triangleattributes
272    meshDict['triangle_neighbors'] = triangleneighbors
273    meshDict['segments'] = segments
274    meshDict['segment_tags'] = segmenttags
275    meshDict['vertex_attribute_titles'] = vertTitle
276   
277    return meshDict
278   
279def _read_outline(fd):
280    """
281    Note, if a file has no mesh info, it can still be read - the meshdic
282    returned will be 'empty'.
283    """
284    delimiter = " " # warning: split() calls are using default whitespace
285   
286    ######### loading the point info
287    line = fd.readline()
288    #print 'read_outline - line',line
289    fragments = line.split()
290    #for fragment in fragments:
291    if fragments ==[]:
292        NumOfVertices = 0
293        NumOfVertAttributes = 0
294    else:
295        NumOfVertices = fragments[0]
296        NumOfVertAttributes = fragments[1]
297    points = []
298    pointattributes = []
299    for index in range(int(NumOfVertices)):       
300        #print index
301        fragments = fd.readline().split() 
302        #print fragments
303        fragments.pop(0) #pop off the index
304        # pop the x y off so we're left with a list of attributes
305        vert = [float(fragments.pop(0)),float(fragments.pop(0))]
306        points.append(vert)
307        apointattributes  = []
308        #print fragments
309        for fragment in fragments:
310            apointattributes.append(float(fragment))
311        pointattributes.append(apointattributes)
312       
313
314    ######### loading the segment info
315    line = fd.readline()
316    fragments = line.split()
317    #for fragment in fragments:
318    #    print fragment
319    if fragments ==[]:
320        NumOfSegments = 0
321    else:
322        NumOfSegments = fragments[0]
323    segments = []
324    segmenttags = []
325    for index in range(int(NumOfSegments)): 
326        #print index
327        line = fd.readline()
328        fragments = line.split()
329        #print fragments
330        fragments.pop(0) #pop off the index
331       
332        seg = [int(fragments[0]),int(fragments[1])]
333        segments.append(seg)
334        line = line[find(line,delimiter):] # remove index
335        line = line.lstrip()
336        line = line[find(line,delimiter):] # remove x
337        line = line.lstrip()
338        line = line[find(line,delimiter):] # remove y
339        stringtag = line.strip()
340        segmenttags.append(stringtag) 
341
342    ######### loading the hole info
343    line = fd.readline()
344    #print line
345    fragments = line.split()
346    #for fragment in fragments:
347    #    print fragment
348    if fragments ==[]:
349        numOfHoles = 0
350    else:
351        numOfHoles = fragments[0]
352    holes = []
353    for index in range(int(numOfHoles)):       
354        #print index
355        fragments = fd.readline().split()
356        #print fragments
357        fragments.pop(0) #pop off the index
358        hole = [float(fragments[0]),float(fragments[1])]
359        holes.append(hole)
360
361   
362    ######### loading the region info
363    line = fd.readline()
364    #print line
365    fragments = line.split()
366    #for fragment in fragments:
367    #    print fragment
368    if fragments ==[]:
369        numOfRegions = 0
370    else:
371        numOfRegions = fragments[0]
372    regions = []
373    regionattributes = []
374    for index in range(int(numOfRegions)):
375        line = fd.readline()
376        fragments = line.split()
377        #print fragments
378        fragments.pop(0) #pop off the index
379        region = [float(fragments[0]),float(fragments[1])]
380        regions.append(region)
381
382        line = line[find(line,delimiter):] # remove index
383        line = line.lstrip()
384        line = line[find(line,delimiter):] # remove x
385        line = line.lstrip()
386        line = line[find(line,delimiter):] # remove y
387        stringtag = line.strip()
388        regionattributes.append(stringtag)
389    regionmaxareas = []
390    line = fd.readline()
391    for index in range(int(numOfRegions)): # Read in the Max area info
392        line = fd.readline()
393        fragments = line.split()
394        #print fragments
395        # The try is here for format compatibility
396        try:
397            fragments.pop(0) #pop off the index
398            if len(fragments) == 0: #no max area
399                regionmaxareas.append(None)
400            else:
401                regionmaxareas.append(float(fragments[0]))
402        except IndexError, e:
403            regionmaxareas.append(None)
404
405    try:
406        geo_reference = Geo_reference(ASCIIFile=fd)
407    except:
408        #geo_ref not compulsory
409        geo_reference = None
410       
411    meshDict = {}
412    meshDict['points'] = points
413    meshDict['point_attributes'] = pointattributes
414    meshDict['outline_segments'] = segments
415    meshDict['outline_segment_tags'] = segmenttags
416    meshDict['holes'] = holes
417    meshDict['regions'] = regions
418    meshDict['region_tags'] = regionattributes
419    meshDict['region_max_areas'] = regionmaxareas
420    meshDict['geo_reference'] = geo_reference
421   
422    #print "meshDict",meshDict
423    return meshDict
424
425def clean_line(line,delimiter):     
426    """Remove whitespace
427    """
428    #print ">%s" %line
429    line = line.strip()
430    #print "stripped>%s" %line
431    numbers = line.split(delimiter)
432    i = len(numbers) - 1
433    while i >= 0:
434        if numbers[i] == '':
435            numbers.pop(i)
436        else:
437            numbers[i] = numbers[i].strip()
438       
439        i += -1
440    #for num in numbers:
441    #    print "num>%s<" %num
442    return numbers
443
444def _write_ASCII_triangulation(fd,
445                               gen_dict):
446    vertices = gen_dict['vertices']
447    vertices_attributes = gen_dict['vertex_attributes']
448    try:
449        vertices_attribute_titles = gen_dict['vertex_attribute_titles']   
450    except KeyError, e:
451        #FIXME is this the best way?
452        if vertices_attributes == [] or vertices_attributes[0] == []:
453             vertices_attribute_titles = []
454        else:
455            raise KeyError, e
456       
457    triangles = gen_dict['triangles']
458    triangles_attributes = gen_dict['triangle_tags']
459    triangle_neighbors = gen_dict['triangle_neighbors']
460       
461    segments = gen_dict['segments']
462    segment_tags = gen_dict['segment_tags']
463     
464    numVert = str(len(vertices))
465    #print "load vertices_attributes", vertices_attributes
466    # Don't understand why we have to do vertices_attributes[0] is None,
467    # but it doesn't work otherwise...
468    if vertices_attributes == None or \
469           (numVert == "0") or \
470           len(vertices_attributes) == 0:
471        numVertAttrib = "0"
472    else:
473        numVertAttrib = str(len(vertices_attributes[0]))
474    fd.write(numVert + " " + numVertAttrib + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Triangulation Vertices..." + "\n")
475
476    #<vertex #> <x> <y> [attributes]   
477    index = 0
478    for vert in vertices:
479        attlist = ""
480       
481        if vertices_attributes == None or \
482               vertices_attributes == []:
483            attlist = ""
484        else:
485            for att in vertices_attributes[index]:
486                attlist = attlist + str(att)+" "
487        attlist.strip()
488        fd.write(str(index) + " "
489                 + str(vert[0]) + " "
490                 + str(vert[1]) + " "
491                 + attlist + "\n")
492        index += 1
493
494    # write comments for title
495    fd.write("# attribute column titles ...Triangulation Vertex Titles..." + "\n")
496    for title in vertices_attribute_titles:
497        fd.write(title + "\n")
498       
499    #<# of triangles>
500    n = len(triangles)
501    fd.write(str(n) + " # <# of triangles>, next lines <triangle #> [<vertex #>] [<neigbouring triangle #>] [attribute of region] ...Triangulation Triangles..." + "\n")
502       
503    # <triangle #> <vertex #>  <vertex #> <vertex #> <neigbouring triangle #> <neigbouring triangle #> <neigbouring triangle #> [attribute of region]
504    for index in range(n):
505        neighbors = ""
506        tri = triangles[index]
507        if triangle_neighbors == []:
508            neighbors = "-1 -1 -1 "
509        else:
510            for neighbor in triangle_neighbors[index]:
511                if neighbor:
512                    neighbors += str(neighbor) + " "
513                else:
514                    if neighbor == 0:
515                        neighbors +=  "0 "
516                    else:
517                        neighbors +=  "-1 "
518        #Warning even though a list is past, only the first value
519        #is written.  There's an assumption that the list only
520        # contains one item. This assumption is made since the
521        #dict that's being passed around is also be used to communicate
522        # with triangle, and it seems to have the option of returning
523        # more than one value for triangle attributex
524        if triangles_attributes == None or \
525               triangles_attributes == [] or \
526               triangles_attributes[index] == ['']:
527            att = ""
528        else:
529            att = str(triangles_attributes[index])
530        fd.write(str(index) + " "
531                 + str(tri[0]) + " " 
532                 + str(tri[1]) + " " 
533                 + str(tri[2]) + " " 
534                 + neighbors + " "
535                 + att + "\n")
536           
537    #One line:  <# of segments>
538    fd.write(str(len(segments)) + 
539             " # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Triangulation Segments..." + "\n")
540       
541    #Following lines:  <segment #> <vertex #>  <vertex #> [boundary tag]
542    for i in range(len(segments)):
543        seg = segments[i]
544        fd.write(str(i) + " "
545                 + str(seg[0]) + " " 
546                 + str(seg[1]) + " " 
547                 + str(segment_tags[i]) + "\n")
548
549
550def _write_tsh_file(ofile,mesh_dict):
551            fd = open(ofile,'w')
552            _write_ASCII_triangulation(fd,mesh_dict)
553            _write_ASCII_outline(fd,mesh_dict)   
554            fd.close()
555
556def _write_ASCII_outline(fd,
557                     dict):
558    points = dict['points']
559    point_attributes = dict['point_attributes']
560    segments = dict['outline_segments']
561    segment_tags = dict['outline_segment_tags']
562    holes = dict['holes']
563    regions = dict['regions']
564    region_tags = dict['region_tags']
565    region_max_areas = dict['region_max_areas']
566       
567    num_points = str(len(points))
568    if (num_points == "0"):
569        num_point_atts = "0"
570    else:
571        num_point_atts = str(len(point_attributes[0]))
572    fd.write(num_points + " " + num_point_atts + " # <# of verts> <# of vert attributes>, next lines <vertex #> <x> <y> [attributes] ...Mesh Vertices..." + "\n")
573       
574    # <x> <y> [attributes]
575    for i,point in enumerate(points):
576        attlist = ""
577        for att in point_attributes[i]:
578            attlist = attlist + str(att)+" "
579        attlist.strip()
580        fd.write(str(i) + " "
581                 +str(point[0]) + " "
582                 + str(point[1]) + " "
583                 + attlist + "\n")
584           
585    #One line:  <# of segments>
586    fd.write(str(len(segments)) + 
587             " # <# of segments>, next lines <segment #> <vertex #>  <vertex #> [boundary tag] ...Mesh Segments..." + "\n")
588       
589    #Following lines:  <vertex #>  <vertex #> [boundary tag]
590    for i,seg in enumerate(segments):
591        fd.write(str(i) + " "
592                 + str(seg[0]) + " " 
593                 + str(seg[1]) + " " 
594                 + str(segment_tags[i]) + "\n")
595
596    #One line:  <# of holes>
597    fd.write(str(len(holes)) + 
598             " # <# of holes>, next lines <Hole #> <x> <y> ...Mesh Holes..." + "\n")   
599    # <x> <y>
600    for i,h in enumerate(holes):
601        fd.write(str(i) + " "
602                 + str(h[0]) + " "
603                 + str(h[1]) + "\n")
604       
605    #One line:  <# of regions>
606    fd.write(str(len(regions)) + 
607             " # <# of regions>, next lines <Region #> <x> <y> <tag>...Mesh Regions..." + "\n")   
608    # <index> <x> <y> <tag>
609    #print "regions",regions
610    for i,r in enumerate(regions):
611        #print "r",r
612        fd.write(str(i) + " " 
613                 + str(r[0]) + " "
614                 + str(r[1])+ " "
615                 + str(region_tags[i]) + "\n")
616       
617    # <index> [<MaxArea>|""]
618       
619    #One line:  <# of regions>
620    fd.write(str(len(regions)) + 
621             " # <# of regions>, next lines <Region #> [Max Area] ...Mesh Regions..." + "\n") 
622    for i,r in enumerate(regions):
623        area = str(region_max_areas[i])
624           
625        fd.write(str(i) + " " + area + "\n")
626
627    # geo_reference info
628    if dict.has_key('geo_reference') and not dict['geo_reference'] is None:
629        dict['geo_reference'].write_ASCII(fd)
630
631def _write_msh_file(file_name, mesh):
632    """Write .msh NetCDF file   
633
634    WARNING: This function mangles the mesh data structure
635    """
636 
637    # FIXME(Ole and John): We ran into a problem on Bogong (64 bit)
638    # where integers appeared as arrays.  This may be similar to
639    # problem seen by Steve in changeset:2778 where he had to wrap
640    # them in int.  Now we are trying with the native Integer format
641    # (Int == 'l' == Int64). However, that caused casting errors, when
642    # 64bit arrays are to be assigned to their NetCDF counterparts. It
643    # seems that the NetCDF arrays are 32bit even though they are
644    # created with the type Int64. Need to look at the NetCDF library
645    # in more detail.
646   
647    IntType = numpy.int32
648    #IntType = Int
649   
650    #
651    #the triangulation
652    mesh['vertices'] = numpy.array(mesh['vertices']).astype(numpy.float)
653    if mesh['vertex_attributes'] != None:
654        mesh['vertex_attributes'] = \
655               numpy.array(mesh['vertex_attributes']).astype(numpy.float)
656    mesh['vertex_attribute_titles'] = numpy.array(mesh['vertex_attribute_titles']).astype(numpy.character) 
657    mesh['segments'] = numpy.array(mesh['segments']).astype(IntType)
658    mesh['segment_tags'] = numpy.array(mesh['segment_tags']).astype(numpy.character)
659    mesh['triangles'] = numpy.array(mesh['triangles']).astype(IntType)
660    mesh['triangle_tags'] = numpy.array(mesh['triangle_tags']) #.astype(numpy.character)
661    mesh['triangle_neighbors'] = numpy.array(mesh['triangle_neighbors']).astype(IntType) 
662   
663    #the outline
664    mesh['points'] = numpy.array(mesh['points']).astype(numpy.float)
665    mesh['point_attributes'] = numpy.array(mesh['point_attributes']).astype(numpy.float)
666    mesh['outline_segments'] = numpy.array(mesh['outline_segments']).astype(IntType)
667    mesh['outline_segment_tags'] = numpy.array(mesh['outline_segment_tags']).astype(numpy.character)
668    mesh['holes'] = numpy.array(mesh['holes']).astype(numpy.float)
669    mesh['regions'] = numpy.array(mesh['regions']).astype(numpy.float)
670    mesh['region_tags'] = numpy.array(mesh['region_tags']).astype(numpy.character)
671    mesh['region_max_areas'] = numpy.array(mesh['region_max_areas']).astype(numpy.float)
672   
673    #mesh = mesh_dict2array(mesh)
674    #print "new_mesh",new_mesh
675    #print "mesh",mesh
676   
677    # NetCDF file definition
678    try:
679        outfile = NetCDFFile(file_name, 'w')
680    except IOError:
681        msg = 'File %s could not be created' %file_name
682        raise msg
683       
684    #Create new file
685    outfile.institution = 'Geoscience Australia'
686    outfile.description = 'NetCDF format for compact and portable storage ' +\
687                      'of spatial point data'
688   
689    # dimension definitions - fixed
690    outfile.createDimension('num_of_dimensions', 2) #This is 2d data
691    outfile.createDimension('num_of_segment_ends', 2) #Segs have two points
692    outfile.createDimension('num_of_triangle_vertices', 3)
693    outfile.createDimension('num_of_triangle_faces', 3)
694    outfile.createDimension('num_of_region_max_area', 1)
695
696    # Create dimensions, variables and set the variables
697   
698    # trianglulation
699    # vertices
700    if (mesh['vertices'].shape[0] > 0):
701        outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0])
702        outfile.createVariable('vertices', Float, ('num_of_vertices',
703                                                   'num_of_dimensions'))
704        outfile.variables['vertices'][:] = mesh['vertices']
705        if mesh['vertex_attributes']  != None and \
706               (mesh['vertex_attributes'].shape[0] > 0 and mesh['vertex_attributes'].shape[1] > 0):
707            outfile.createDimension('num_of_vertex_attributes',
708                                    mesh['vertex_attributes'].shape[1])
709            outfile.createDimension('num_of_vertex_attribute_title_chars',
710                                    mesh['vertex_attribute_titles'].shape[1])
711            outfile.createVariable('vertex_attributes',
712                                   Float,
713                                   ('num_of_vertices',
714                                    'num_of_vertex_attributes'))
715            outfile.createVariable('vertex_attribute_titles',
716                                   Character,
717                                   ( 'num_of_vertex_attributes',
718                                     'num_of_vertex_attribute_title_chars' ))
719            outfile.variables['vertex_attributes'][:] = \
720                                                      mesh['vertex_attributes']
721            outfile.variables['vertex_attribute_titles'][:] = \
722                                     mesh['vertex_attribute_titles']
723    # segments
724    if (mesh['segments'].shape[0] > 0):       
725        outfile.createDimension('num_of_segments',
726                                mesh['segments'].shape[0])
727        outfile.createVariable('segments',
728                               Int,
729                               ('num_of_segments', 'num_of_segment_ends'))
730        outfile.variables['segments'][:] = mesh['segments']
731        if (mesh['segment_tags'].shape[1] > 0):
732            outfile.createDimension('num_of_segment_tag_chars',
733                                    mesh['segment_tags'].shape[1])
734            outfile.createVariable('segment_tags',
735                                   Character,
736                                   ('num_of_segments',
737                                    'num_of_segment_tag_chars'))
738            outfile.variables['segment_tags'][:] = mesh['segment_tags']
739    # triangles   
740    if (mesh['triangles'].shape[0] > 0):
741        outfile.createDimension('num_of_triangles',
742                                mesh['triangles'].shape[0])
743        outfile.createVariable('triangles',
744                               Int,
745                               ('num_of_triangles',
746                                'num_of_triangle_vertices'))
747        outfile.createVariable('triangle_neighbors',
748                               Int,
749                               ('num_of_triangles',
750                                'num_of_triangle_faces'))
751        outfile.variables['triangles'][:] = mesh['triangles']
752        outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors']
753        if mesh['triangle_tags'] != None and \
754               (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                               Int,
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    mesh = {}
846    # Get the variables
847    # the triangulation
848    try:
849        mesh['vertices'] = fid.variables['vertices'][:]
850    except KeyError:
851        mesh['vertices'] = numpy.array([])
852    try:
853        mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:]
854    except KeyError:
855        mesh['vertex_attributes'] = None
856        #for ob in mesh['vertices']:
857            #mesh['vertex_attributes'].append([])
858   
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'] = numpy.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'] = numpy.array([])
883        mesh['triangle_neighbors'] = numpy.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'] = numpy.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'] = numpy.array([])
920    try:
921        mesh['regions'] = fid.variables['regions'][:]
922    except KeyError:
923        mesh['regions'] = numpy.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'] = numpy.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 = numpy.array(points).astype(numpy.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 = numpy.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'] = numpy.array(point_atts['pointlist']).astype(numpy.float)
1031   
1032    for key in point_atts['attributelist'].keys():
1033        point_atts['attributelist'][key]= numpy.array(point_atts['attributelist'][key]).astype(numpy.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 = numpy.array([]).astype(numpy.float)
1050    keys = dic.keys()
1051    key = keys.pop(0)
1052    point_attributes = numpy.reshape(dic[key],(dic[key].shape[0],1))
1053    for key in keys:
1054        #point_attributes = numpy.concatenate([point_attributes, dic[key]], axis=1)
1055        reshaped = numpy.reshape(dic[key],(dic[key].shape[0],1))
1056        point_attributes = numpy.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'] = numpy.take(dict['pointlist'],indices_to_keep)
1064
1065    for key in dict['attributelist'].keys():
1066        dict['attributelist'][key]= numpy.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'] = numpy.concatenate((dict2['pointlist'],
1078                                         dict1['pointlist']),axis=0)   
1079    atts = {}   
1080    for key in dict2['attributelist'].keys():
1081        atts[key]= numpy.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.