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

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

Ongoing conversion changes.

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