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

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

numpy changes

File size: 39.8 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[0] > 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:
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:
836        pass
837
838    try:
839        mesh['triangles'] = fid.variables['triangles'][:]
840        mesh['triangle_neighbors'] = fid.variables['triangle_neighbors'][:]
841    except KeyError:
842        mesh['triangles'] = num.array([], num.int)              #array default#
843        mesh['triangle_neighbors'] = num.array([], num.int)     #array default#
844
845    mesh['triangle_tags'] = []
846    try:
847        tags = fid.variables['triangle_tags'][:]
848        mesh['triangle_tags'] = [x.tostring().strip() for x in tags]
849    except KeyError:
850        for ob in mesh['triangles']:
851            mesh['triangle_tags'].append('')
852
853    #the outline
854    try:
855        mesh['points'] = fid.variables['points'][:]
856    except KeyError:
857        mesh['points'] = []
858
859    try:
860        mesh['point_attributes'] = fid.variables['point_attributes'][:]
861    except KeyError:
862        mesh['point_attributes'] = []
863        for point in mesh['points']:
864            mesh['point_attributes'].append([])
865
866    try:
867        mesh['outline_segments'] = fid.variables['outline_segments'][:]
868    except KeyError:
869        mesh['outline_segments'] = num.array([], num.int)      #array default#
870
871    mesh['outline_segment_tags'] =[]
872    try:
873        tags = fid.variables['outline_segment_tags'][:]
874        for i, tag in enumerate(tags):
875            mesh['outline_segment_tags'].append(tags[i].tostring().strip())
876    except KeyError:
877        for ob in mesh['outline_segments']:
878            mesh['outline_segment_tags'].append('')
879
880    try:
881        mesh['holes'] = fid.variables['holes'][:]
882    except KeyError:
883        mesh['holes'] = num.array([], num.int)          #array default#
884
885    try:
886        mesh['regions'] = fid.variables['regions'][:]
887    except KeyError:
888        mesh['regions'] = num.array([], num.int)        #array default#
889
890    mesh['region_tags'] =[]
891    try:
892        tags = fid.variables['region_tags'][:]
893        for i, tag in enumerate(tags):
894            mesh['region_tags'].append(tags[i].tostring().strip())
895    except KeyError:
896        for ob in mesh['regions']:
897            mesh['region_tags'].append('')
898
899    try:
900        mesh['region_max_areas'] = fid.variables['region_max_areas'][:]
901    except KeyError:
902        mesh['region_max_areas'] = num.array([], num.int)      #array default#
903
904    try:
905        geo_reference = Geo_reference(NetCDFObject=fid)
906        mesh['geo_reference'] = geo_reference
907    except AttributeError, e:
908        #geo_ref not compulsory
909        mesh['geo_reference'] = None
910
911    fid.close()
912
913    return mesh
914
915
916##
917# @brief Export a boundary file.
918# @param file_name Path to the file to write.
919# @param points List of point index pairs.
920# @paran title Title string to write into file (first line).
921# @param delimiter Delimiter string to write between points.
922# @note Used by alpha shapes
923def export_boundary_file(file_name, points, title, delimiter=','):
924    """Export a boundary file.
925
926    Format:
927    First line: Title variable
928    Following lines:  [point index][delimiter][point index]
929
930    file_name - the name of the new file
931    points - List of point index pairs [[p1, p2],[p3, p4]..]
932    title - info to write in the first line
933    """
934
935    fd = open(file_name, 'w')
936
937    fd.write(title + '\n')
938
939    # [point index][delimiter][point index]
940    for point in points:
941        fd.write(str(point[0]) + delimiter + str(point[1]) + '\n')
942
943    fd.close()
944
945################################################################################
946#  IMPORT/EXPORT POINTS FILES
947################################################################################
948
949##
950# @brief
951# @param point_atts
952def extent_point_atts(point_atts):
953    """Returns 4 points representing the extent
954    This loses attribute info.
955    """
956
957    point_atts['pointlist'] = extent(point_atts['pointlist'])
958    point_atts['attributelist'] = {}
959
960    return point_atts
961
962
963##
964# @brief Get an extent array for a set of points.
965# @param points A set of points.
966# @return An extent array of form [[min_x, min_y], [max_x, min_y],
967#                                  [max_x, max_y], [min_x, max_y]]
968def extent(points):
969    points = num.array(points, num.float)
970
971    max_x = min_x = points[0][0]
972    max_y = min_y = points[0][1]
973
974    for point in points[1:]:
975        x = point[0]
976        if x > max_x:
977            max_x = x
978        if x < min_x:
979            min_x = x
980
981        y = point[1]
982        if y > max_y:
983            max_y = y
984        if y < min_y:
985            min_y = y
986
987    extent = num.array([[min_x, min_y],
988                        [max_x, min_y],
989                        [max_x, max_y],
990                        [min_x, max_y]])
991
992    return extent
993
994
995##
996# @brief Reduce a points file until number of points is less than limit.
997# @param infile Path to input file to thin.
998# @param outfile Path to output thinned file.
999# @param max_points Max number of points in output file.
1000# @param verbose True if this function is to be verbose.
1001def reduce_pts(infile, outfile, max_points, verbose = False):
1002    """Reduce a points file until less than given size.
1003
1004    Reduces a points file by removing every second point until the # of points
1005    is less than max_points.
1006    """
1007
1008    # check out pts2rectangular in least squares, and the use of reduction.
1009    # Maybe it does the same sort of thing?
1010    point_atts = _read_pts_file(infile)
1011
1012    while point_atts['pointlist'].shape[0] > max_points:
1013        if verbose: print "point_atts['pointlist'].shape[0]"
1014        point_atts = half_pts(point_atts)
1015
1016    export_points_file(outfile, point_atts)
1017
1018
1019##
1020# @brief
1021# @param infile
1022# @param max_points
1023# @param delimiter
1024# @param verbose True if this function is to be verbose.
1025# @return A list of
1026def produce_half_point_files(infile, max_points, delimiter, verbose=False):
1027    point_atts = _read_pts_file(infile)
1028    root, ext = splitext(infile)
1029    outfiles = []
1030
1031    if verbose: print "# of points", point_atts['pointlist'].shape[0]
1032
1033    while point_atts['pointlist'].shape[0] > max_points:
1034        point_atts = half_pts(point_atts)
1035
1036        if verbose: print "# of points", point_atts['pointlist'].shape[0]
1037
1038        outfile = root + delimiter + str(point_atts['pointlist'].shape[0]) + ext
1039        outfiles.append(outfile)
1040        export_points_file(outfile, point_atts)
1041
1042    return outfiles
1043
1044
1045##
1046# @brief
1047# @param point_atts Object with attribute of 'pointlist'.
1048# @return
1049def point_atts2array(point_atts):
1050    # convert attribute list to array of floats
1051    point_atts['pointlist'] = num.array(point_atts['pointlist'], num.float)
1052
1053    for key in point_atts['attributelist'].keys():
1054        point_atts['attributelist'][key] = \
1055            num.array(point_atts['attributelist'][key], num.float)
1056
1057    return point_atts
1058
1059
1060##
1061# @brief ??
1062# @param point_atts ??
1063# @return ??
1064def half_pts(point_atts):
1065    point_atts2array(point_atts)
1066    point_atts['pointlist'] = point_atts['pointlist'][::2]
1067
1068    for key in point_atts['attributelist'].keys():
1069        point_atts['attributelist'][key] = point_atts['attributelist'][key][::2]
1070
1071    return point_atts
1072
1073##
1074# @brief ??
1075# @param dic ??
1076# @return ??
1077def concatinate_attributelist(dic):
1078    """
1079    giving a dic[attribute title] = attribute
1080    return list of attribute titles, array of attributes
1081    """
1082
1083    point_attributes = num.array([], num.float)
1084    keys = dic.keys()
1085    key = keys.pop(0)
1086    point_attributes = num.reshape(dic[key], (dic[key].shape[0], 1))
1087    for key in keys:
1088        reshaped = num.reshape(dic[key], (dic[key].shape[0], 1))
1089        point_attributes = num.concatenate([point_attributes, reshaped], axis=1)
1090
1091    return dic.keys(), point_attributes
1092
1093
1094##
1095# @brief
1096# @param dict
1097# @param indices_to_keep
1098# @return
1099# @note FIXME(dsg), turn this dict plus methods into a class?
1100def take_points(dict, indices_to_keep):
1101    dict = point_atts2array(dict)
1102    dict['pointlist'] = num.take(dict['pointlist'], indices_to_keep, axis=0)
1103
1104    for key in dict['attributelist'].keys():
1105        dict['attributelist'][key] = num.take(dict['attributelist'][key],
1106                                              indices_to_keep, axis=0)
1107
1108    return dict
1109
1110##
1111# @brief ??
1112# @param dict1 ??
1113# @param dict2 ??
1114# @return ??
1115def add_point_dictionaries(dict1, dict2):
1116    """
1117    """
1118
1119    dict1 = point_atts2array(dict1)
1120    dict2 = point_atts2array(dict2)
1121
1122    combined = {}
1123    combined['pointlist'] = num.concatenate((dict2['pointlist'],
1124                                             dict1['pointlist']), axis=0)
1125
1126    atts = {}
1127    for key in dict2['attributelist'].keys():
1128        atts[key]= num.concatenate((dict2['attributelist'][key],
1129                                    dict1['attributelist'][key]), axis=0)
1130    combined['attributelist'] = atts
1131    combined['geo_reference'] = dict1['geo_reference']
1132
1133    return combined
1134
1135
1136if __name__ == "__main__":
1137    pass
Note: See TracBrowser for help on using the repository browser.