source: anuga_core/source/anuga/load_mesh/loadASCII.py @ 6074

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

Looked at NetCDF code, also did a pep8.

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