source: branches/numpy/anuga/load_mesh/loadASCII.py @ 7176

Last change on this file since 7176 was 7176, checked in by rwilson, 15 years ago

Back-merge from Numeric to numpy.

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