source: branches/Numeric_anuga_source/anuga/load_mesh/loadASCII.py @ 9065

Last change on this file since 9065 was 6166, checked in by rwilson, 16 years ago

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