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

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

Initial commit of numpy changes. Still a long way to go.

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