""" The format for a Points dictionary is: ['pointlist'] a 2 column array describing points. 1st column x, 2nd column y. ['attributelist'], a dictionary of 1D arrays, representing attribute values at the point. The dictionary key is the attribute header. ['geo_reference'] a Geo_refernece object. Use if the point information is relative. This is optional. eg dic['pointlist'] = [[1.0,2.0],[3.0,5.0]] dic['attributelist']['elevation'] = [[7.0,5.0] The dict format for IO with mesh files is; (the triangulation) vertices: [[x1,y1],[x2,y2],...] (lists of doubles) vertex_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles) vertex_attribute_titles:[A1Title, A2Title ...] (A list of strings) segments: [[v1,v2],[v3,v4],...] (lists of integers) segment_tags : [tag,tag,...] list of strings triangles : [(v1,v2,v3), (v4,v5,v6),....] lists of points triangle_tags: [s1,s2,...] A list of strings triangle_neighbors: [[t1,t2,t3], [t4,t5,t6],..] lists of triangles (the outline) points: [[x1,y1],[x2,y2],...] (lists of doubles) point_attributes: [[a11,a12,...],[a21,a22],...] (lists of doubles) outline_segments: [[point1,point2],[p3,p4],...] (lists of integers) outline_segment_tags : [tag1,tag2,...] list of strings holes : [[x1,y1],...](List of doubles, one inside each hole region) regions : [ [x1,y1],...] (List of 4 doubles) region_tags : [tag1,tag2,...] (list of strings) region_max_areas: [ma1,ma2,...] (A list of doubles) {Convension: A -ve max area means no max area} geo_reference: a Geo_refernece object. Use if the point information is relative. This is optional. Points files are .csv for ascii and .pts for NetCDF Mesh files are .tsh for ascii and .msh for NetCDF Note: point att's have no titles - that's currently ok, since least squares adds the point info to vertices, and they have titles The format for a .tsh file is (FIXME update this) First line: <# of vertices> <# of attributes> Following lines: [attributes] One line: <# of triangles> Following lines: [attribute of region] One line: <# of segments> Following lines: [boundary tag] """ ##FIXME (DSG-DSG) Is the dict format mentioned above a list of a numeric array? # Needs to be defined from string import find, rfind import numpy as num from os.path import splitext from anuga.coordinate_transforms.geo_reference import Geo_reference, TITLE, \ TitleError from Scientific.IO.NetCDF import NetCDFFile from anuga.config import netcdf_mode_r, netcdf_mode_w, netcdf_mode_a from anuga.config import netcdf_float, netcdf_char, netcdf_int import exceptions class TitleAmountError(exceptions.Exception): pass NOMAXAREA=-999 ## # @brief Read a mesh file (.tsh or .msh) and return a dictionary. # @param ofile Path to the file to read. # @return A dictionary of data from the input file. # @note Throws IOError if can't read file. # @note This is used by pmesh. def import_mesh_file(ofile): """Read a mesh file, either .tsh or .msh Note: will throw an IOError if it can't load the file. Catch these! """ try: if ofile[-4:] == ".tsh": dict = _read_tsh_file(ofile) elif ofile[-4:] == ".msh": dict = _read_msh_file(ofile) else: msg = 'Extension .%s is unknown' % ofile[-4:] raise IOError, msg #FIXME No test for ValueError except (TitleError, SyntaxError, IndexError, ValueError): msg = 'File could not be opened' raise IOError, msg return dict ## # @brief Write a mesh file from a dictionary. # @param ofile The path of the mesh file to write. # @param mesh_dict Dictionary containing data to write to output mesh file. # @note Raises IOError if file extension is unrecognized. def export_mesh_file(ofile, mesh_dict): """Write a mesh file given a dictionary. First line: <# of vertices> <# of attributes> Following lines: [attributes] One line: <# of triangles> Following lines: [attribute of region] One line: <# of segments> Following lines: [boundary tag] """ # Ensure that if required key isn't present, we add it with [] value. # .setdefault() for dictionaries was added in python 2.0. # The only other neat way to do this is to override the dictionary # .__getitem__() method to return [] if key not found. reqd_keys = ['points', 'point_attributes', 'outline_segments', 'outline_segment_tags', 'holes', 'regions', 'region_tags', 'region_max_areas', 'vertices', 'vertex_attributes', 'vertex_attribute_titles', 'segments', 'segment_tags', 'triangles', 'triangle_tags', 'triangle_neighbors'] for k in reqd_keys: mesh_dict.setdefault(k, []) # hand-off to appropriate function if (ofile[-4:] == ".tsh"): _write_tsh_file(ofile, mesh_dict) elif (ofile[-4:] == ".msh"): _write_msh_file(ofile, mesh_dict) else: msg = 'Unknown file type %s ' % ofile raise IOError, msg ## # @brief Read a mesh file into a dictionary. # @param ofile Path of the file to read. # @return Data dictionary from the mesh (.tsh) file. def _read_tsh_file(ofile): """Read the text file format for meshes""" fd = open(ofile, 'r') dict = _read_triangulation(fd) dict_mesh = _read_outline(fd) for element in dict_mesh.keys(): dict[element] = dict_mesh[element] fd.close() return dict ## # @brief ?? # @param fd An open descriptor of the file to read. # @return A dictionary ... def _read_triangulation(fd): """Read the generated triangulation, NOT the outline.""" delimiter = " " # loading the point info line = fd.readline() fragments = line.split() if fragments == []: NumOfVertices = 0 NumOfVertAttributes = 0 else: NumOfVertices = fragments[0] NumOfVertAttributes = fragments[1] points = [] pointattributes = [] for index in range(int(NumOfVertices)): fragments = fd.readline().split() fragments.pop(0) # pop off the index # pop the x & y off so we're left with a list of attributes vert = [float(fragments.pop(0)), float(fragments.pop(0))] points.append(vert) apointattributes = [] for fragment in fragments: apointattributes.append(float(fragment)) if apointattributes != []: pointattributes.append(apointattributes) # loading the point title info line = fd.readline() vertTitle = [] for index in range(int(NumOfVertAttributes)): fragments = fd.readline().strip() vertTitle.append(fragments) # loading the triangle info line = fd.readline() fragments = line.split() NumOfTriangles = fragments[0] triangles = [] triangleattributes = [] triangleneighbors = [] for index in range(int(NumOfTriangles)): line = fd.readline() line.strip() # so we can get the region string fragments = line.split() fragments.pop(0) # pop off the index tri = [int(fragments[0]), int(fragments[1]), int(fragments[2])] triangles.append(tri) neighbors = [int(fragments[3]), int(fragments[4]), int(fragments[5])] triangleneighbors.append(neighbors) for x in range(7): # remove index [] [] line = line[find(line, delimiter):] # remove index line = line.lstrip() stringtag = line.strip() triangleattributes.append(stringtag) # loading the segment info line = fd.readline() fragments = line.split() NumOfSegments = fragments[0] segments = [] segmenttags = [] for index in range(int(NumOfSegments)): line = fd.readline() line.strip() # to get the segment string fragments = line.split() fragments.pop(0) #pop off the index seg = [int(fragments[0]), int(fragments[1])] segments.append(seg) line = line[find(line, delimiter):] # remove index line = line.lstrip() line = line[find(line, delimiter):] # remove x line = line.lstrip() line = line[find(line, delimiter):] # remove y stringtag = line.strip() segmenttags.append(stringtag) meshDict = {} meshDict['vertices'] = points if pointattributes == []: meshDict['vertex_attributes'] = None else: meshDict['vertex_attributes'] = pointattributes meshDict['triangles'] = triangles meshDict['triangle_tags'] = triangleattributes meshDict['triangle_neighbors'] = triangleneighbors meshDict['segments'] = segments meshDict['segment_tags'] = segmenttags meshDict['vertex_attribute_titles'] = vertTitle return meshDict ## # @brief # @param fd An open descriptor of the file to read. # @return A dictionary ... # @note If file has no mesh info, an empty dict will be returned. def _read_outline(fd): """Read a mesh file outline. Note, if a file has no mesh info, it can still be read - the meshdic returned will be 'empty'. """ delimiter = " " # warning: split() calls are using default whitespace # loading the point info line = fd.readline() fragments = line.split() if fragments == []: NumOfVertices = 0 NumOfVertAttributes = 0 else: NumOfVertices = fragments[0] NumOfVertAttributes = fragments[1] points = [] pointattributes = [] for index in range(int(NumOfVertices)): fragments = fd.readline().split() fragments.pop(0) # pop off the index # pop the x & y off so we're left with a list of attributes vert = [float(fragments.pop(0)), float(fragments.pop(0))] points.append(vert) apointattributes = [] for fragment in fragments: apointattributes.append(float(fragment)) pointattributes.append(apointattributes) # loading the segment info line = fd.readline() fragments = line.split() if fragments == []: NumOfSegments = 0 else: NumOfSegments = fragments[0] segments = [] segmenttags = [] for index in range(int(NumOfSegments)): line = fd.readline() fragments = line.split() fragments.pop(0) # pop off the index seg = [int(fragments[0]), int(fragments[1])] segments.append(seg) line = line[find(line, delimiter):] # remove index line = line.lstrip() line = line[find(line, delimiter):] # remove x line = line.lstrip() line = line[find(line, delimiter):] # remove y stringtag = line.strip() segmenttags.append(stringtag) # loading the hole info line = fd.readline() fragments = line.split() if fragments == []: numOfHoles = 0 else: numOfHoles = fragments[0] holes = [] for index in range(int(numOfHoles)): fragments = fd.readline().split() fragments.pop(0) #pop off the index hole = [float(fragments[0]), float(fragments[1])] holes.append(hole) # loading the region info line = fd.readline() fragments = line.split() if fragments == []: numOfRegions = 0 else: numOfRegions = fragments[0] regions = [] regionattributes = [] for index in range(int(numOfRegions)): line = fd.readline() fragments = line.split() fragments.pop(0) # pop off the index region = [float(fragments[0]), float(fragments[1])] regions.append(region) line = line[find(line, delimiter):] # remove index line = line.lstrip() line = line[find(line, delimiter):] # remove x line = line.lstrip() line = line[find(line, delimiter):] # remove y stringtag = line.strip() regionattributes.append(stringtag) regionmaxareas = [] line = fd.readline() for index in range(int(numOfRegions)): # Read in the Max area info line = fd.readline() fragments = line.split() # The try is here for format compatibility try: fragments.pop(0) # pop off the index if len(fragments) == 0: # no max area regionmaxareas.append(None) else: regionmaxareas.append(float(fragments[0])) except IndexError, e: regionmaxareas.append(None) try: geo_reference = Geo_reference(ASCIIFile=fd) except: #geo_ref not compulsory geo_reference = None meshDict = {} meshDict['points'] = points meshDict['point_attributes'] = pointattributes meshDict['outline_segments'] = segments meshDict['outline_segment_tags'] = segmenttags meshDict['holes'] = holes meshDict['regions'] = regions meshDict['region_tags'] = regionattributes meshDict['region_max_areas'] = regionmaxareas meshDict['geo_reference'] = geo_reference return meshDict ## # @brief Clean up a line with delimited fields. # @param line The string to be cleaned. # @param delimiter String used a delimiter when splitting 'line'. # @return A list of delimited fields with white-space removed from ends. # @note Will also remove any field that was originally ''. def clean_line(line, delimiter): """Clean up a line - return fields with end space removed.""" line = line.strip() numbers = line.split(delimiter) i = len(numbers) - 1 while i >= 0: if numbers[i] == '': numbers.pop(i) else: numbers[i] = numbers[i].strip() i += -1 return numbers ## # @brief Write an ASCII triangulation file. # @param fd An open descriptor to the file to write. # @param gen_dict Dictionary containing data to write. def _write_ASCII_triangulation(fd, gen_dict): vertices = gen_dict['vertices'] vertices_attributes = gen_dict['vertex_attributes'] try: vertices_attribute_titles = gen_dict['vertex_attribute_titles'] except KeyError, e: #FIXME is this the best way? if vertices_attributes == [] or vertices_attributes[0] == []: vertices_attribute_titles = [] else: raise KeyError, e triangles = gen_dict['triangles'] triangles_attributes = gen_dict['triangle_tags'] triangle_neighbors = gen_dict['triangle_neighbors'] segments = gen_dict['segments'] segment_tags = gen_dict['segment_tags'] numVert = str(len(vertices)) # Don't understand why we have to do vertices_attributes[0] is None, # but it doesn't work otherwise... if (vertices_attributes == None or numVert == "0" or len(vertices_attributes) == 0): numVertAttrib = "0" else: numVertAttrib = str(len(vertices_attributes[0])) fd.write(numVert + " " + numVertAttrib + " # <# of verts> <# of vert attributes>, next lines " " [attributes] ...Triangulation Vertices...\n") # [attributes] index = 0 for vert in vertices: attlist = "" if vertices_attributes == None or vertices_attributes == []: attlist = "" else: for att in vertices_attributes[index]: attlist = attlist + str(att) + " " attlist.strip() fd.write(str(index) + " " + str(vert[0]) + " " + str(vert[1]) + " " + attlist + "\n") index += 1 # write comments for title fd.write("# attribute column titles ...Triangulation Vertex Titles...\n") for title in vertices_attribute_titles: fd.write(title + "\n") # <# of triangles> n = len(triangles) fd.write(str(n) + " # <# of triangles>, next lines [] " "[] [attribute of region] ..." "Triangulation Triangles...\n") # # [attribute of region] for index in range(n): neighbors = "" tri = triangles[index] if triangle_neighbors == []: neighbors = "-1 -1 -1 " else: for neighbor in triangle_neighbors[index]: if neighbor: neighbors += str(neighbor) + " " else: if neighbor == 0: neighbors += "0 " else: neighbors += "-1 " # Warning even though a list is passed, only the first value # is written. There's an assumption that the list only # contains one item. This assumption is made since the # dict that's being passed around is also be used to communicate # with triangle, and it seems to have the option of returning # more than one value for triangle attributex if (triangles_attributes == None or triangles_attributes == [] or triangles_attributes[index] == ['']): att = "" else: att = str(triangles_attributes[index]) fd.write(str(index) + " " + str(tri[0]) + " " + str(tri[1]) + " " + str(tri[2]) + " " + neighbors + " " + att + "\n") # One line: <# of segments> fd.write(str(len(segments)) + " # <# of segments>, next lines " " [boundary tag] ...Triangulation Segments...\n") # Following lines: [boundary tag] for i in range(len(segments)): seg = segments[i] fd.write(str(i) + " " + str(seg[0]) + " " + str(seg[1]) + " " + str(segment_tags[i]) + "\n") ## # @brief Write triangulation and outline to a .tsh file. # @param ofile Path of the file to write. # @mesh_dict Dictionary of data to write. def _write_tsh_file(ofile, mesh_dict): fd = open(ofile, 'w') _write_ASCII_triangulation(fd, mesh_dict) _write_ASCII_outline(fd, mesh_dict) fd.close() ## # @brief Write an ASCII outline to a file. # @param fd An open descriptor of the file to write. # @param dict Dictionary holding data to write. def _write_ASCII_outline(fd, dict): points = dict['points'] point_attributes = dict['point_attributes'] segments = dict['outline_segments'] segment_tags = dict['outline_segment_tags'] holes = dict['holes'] regions = dict['regions'] region_tags = dict['region_tags'] region_max_areas = dict['region_max_areas'] num_points = str(len(points)) if (num_points == '0'): num_point_atts = '0' else: num_point_atts = str(len(point_attributes[0])) fd.write(num_points + ' ' + num_point_atts + ' # <# of verts> <# of vert attributes>, next lines ' ' [attributes] ...Mesh Vertices...\n') # [attributes] for i, point in enumerate(points): attlist = '' for att in point_attributes[i]: attlist = attlist + str(att) + ' ' attlist.strip() fd.write(str(i) + ' ' + str(point[0]) + ' ' + str(point[1]) + ' ' + attlist + '\n') # One line: <# of segments> fd.write(str(len(segments)) + ' # <# of segments>, next lines ' ' [boundary tag] ...Mesh Segments...\n') # Following lines: [boundary tag] for i,seg in enumerate(segments): fd.write(str(i) + ' ' + str(seg[0]) + ' ' + str(seg[1]) + ' ' + str(segment_tags[i]) + '\n') # One line: <# of holes> fd.write(str(len(holes)) + ' # <# of holes>, next lines ...Mesh Holes...\n') # for i,h in enumerate(holes): fd.write(str(i) + ' ' + str(h[0]) + ' ' + str(h[1]) + '\n') # One line: <# of regions> fd.write(str(len(regions)) + ' # <# of regions>, next lines ' '...Mesh Regions...\n') # for i,r in enumerate(regions): fd.write(str(i) + ' ' + str(r[0]) + ' ' + str(r[1])+ ' ' + str(region_tags[i]) + '\n') # [|''] # One line: <# of regions> fd.write(str(len(regions)) + ' # <# of regions>, next lines [Max Area] ' '...Mesh Regions...\n') for i,r in enumerate(regions): area = str(region_max_areas[i]) fd.write(str(i) + ' ' + area + '\n') # geo_reference info if dict.has_key('geo_reference') and not dict['geo_reference'] is None: dict['geo_reference'].write_ASCII(fd) ## # @brief Write a .msh file. # @param file Path to the file to write. # @param mesh Dictionary holding data to write. def _write_msh_file(file_name, mesh): """Write .msh NetCDF file WARNING: This function mangles the mesh data structure """ # FIXME(Ole and John): We ran into a problem on Bogong (64 bit) # where integers appeared as arrays. This may be similar to # problem seen by Steve in changeset:2778 where he had to wrap # them in int. Now we are trying with the native Integer format # (Int == 'l' == Int64). However, that caused casting errors, when # 64bit arrays are to be assigned to their NetCDF counterparts. It # seems that the NetCDF arrays are 32bit even though they are # created with the type Int64. Need to look at the NetCDF library # in more detail. IntType = num.int32 #IntType = Int #the triangulation mesh['vertices'] = num.array(mesh['vertices'], num.float) if mesh['vertex_attributes'] != None: mesh['vertex_attributes'] = \ num.array(mesh['vertex_attributes'], num.float) mesh['vertex_attribute_titles'] = \ num.array(mesh['vertex_attribute_titles'], num.character) mesh['segments'] = num.array(mesh['segments'], IntType) mesh['segment_tags'] = num.array(mesh['segment_tags'], num.character) mesh['triangles'] = num.array(mesh['triangles'], IntType) mesh['triangle_tags'] = num.array(mesh['triangle_tags']) mesh['triangle_neighbors'] = \ num.array(mesh['triangle_neighbors'], IntType) #the outline mesh['points'] = num.array(mesh['points'], num.float) mesh['point_attributes'] = num.array(mesh['point_attributes'], num.float) mesh['outline_segments'] = num.array(mesh['outline_segments'], IntType) mesh['outline_segment_tags'] = \ num.array(mesh['outline_segment_tags'], num.character) mesh['holes'] = num.array(mesh['holes'], num.float) mesh['regions'] = num.array(mesh['regions'], num.float) mesh['region_tags'] = num.array(mesh['region_tags'], num.character) mesh['region_max_areas'] = num.array(mesh['region_max_areas'], num.float) # NetCDF file definition try: outfile = NetCDFFile(file_name, netcdf_mode_w) except IOError: msg = 'File %s could not be created' % file_name raise Exception, msg #Create new file outfile.institution = 'Geoscience Australia' outfile.description = 'NetCDF format for compact and portable storage ' + \ 'of spatial point data' # dimension definitions - fixed outfile.createDimension('num_of_dimensions', 2) # This is 2d data outfile.createDimension('num_of_segment_ends', 2) # Segs have two points outfile.createDimension('num_of_triangle_vertices', 3) outfile.createDimension('num_of_triangle_faces', 3) outfile.createDimension('num_of_region_max_area', 1) # Create dimensions, variables and set the variables # trianglulation # vertices if (mesh['vertices'].shape[0] > 0): outfile.createDimension('num_of_vertices', mesh['vertices'].shape[0]) outfile.createVariable('vertices', netcdf_float, ('num_of_vertices', 'num_of_dimensions')) outfile.variables['vertices'][:] = mesh['vertices'] if (mesh['vertex_attributes'] != None and (mesh['vertex_attributes'].shape[0] > 0 and mesh['vertex_attributes'].shape[1] > 0)): outfile.createDimension('num_of_vertex_attributes', mesh['vertex_attributes'].shape[1]) outfile.createDimension('num_of_vertex_attribute_title_chars', mesh['vertex_attribute_titles'].shape[1]) outfile.createVariable('vertex_attributes', netcdf_float, ('num_of_vertices', 'num_of_vertex_attributes')) outfile.createVariable('vertex_attribute_titles', netcdf_char, ('num_of_vertex_attributes', 'num_of_vertex_attribute_title_chars')) outfile.variables['vertex_attributes'][:] = \ mesh['vertex_attributes'] outfile.variables['vertex_attribute_titles'][:] = \ mesh['vertex_attribute_titles'] # segments if (mesh['segments'].shape[0] > 0): outfile.createDimension('num_of_segments', mesh['segments'].shape[0]) outfile.createVariable('segments', netcdf_int, ('num_of_segments', 'num_of_segment_ends')) outfile.variables['segments'][:] = mesh['segments'] if mesh['segment_tags'].shape[1] > 0: outfile.createDimension('num_of_segment_tag_chars', mesh['segment_tags'].shape[1]) outfile.createVariable('segment_tags', netcdf_char, ('num_of_segments', 'num_of_segment_tag_chars')) outfile.variables['segment_tags'][:] = mesh['segment_tags'] # triangles if (mesh['triangles'].shape[0] > 0): outfile.createDimension('num_of_triangles', mesh['triangles'].shape[0]) outfile.createVariable('triangles', netcdf_int, ('num_of_triangles', 'num_of_triangle_vertices')) outfile.createVariable('triangle_neighbors', netcdf_int, ('num_of_triangles', 'num_of_triangle_faces')) outfile.variables['triangles'][:] = mesh['triangles'] outfile.variables['triangle_neighbors'][:] = mesh['triangle_neighbors'] if (mesh['triangle_tags'] != None and (mesh['triangle_tags'].shape[1] > 0)): outfile.createDimension('num_of_triangle_tag_chars', mesh['triangle_tags'].shape[1]) outfile.createVariable('triangle_tags', netcdf_char, ('num_of_triangles', 'num_of_triangle_tag_chars')) outfile.variables['triangle_tags'][:] = mesh['triangle_tags'] # outline # points if (mesh['points'].shape[0] > 0): outfile.createDimension('num_of_points', mesh['points'].shape[0]) outfile.createVariable('points', netcdf_float, ('num_of_points', 'num_of_dimensions')) outfile.variables['points'][:] = mesh['points'] if mesh['point_attributes'].shape[0] > 0 \ and mesh['point_attributes'].shape[1] > 0: outfile.createDimension('num_of_point_attributes', mesh['point_attributes'].shape[1]) outfile.createVariable('point_attributes', netcdf_float, ('num_of_points', 'num_of_point_attributes')) outfile.variables['point_attributes'][:] = mesh['point_attributes'] # outline_segments if mesh['outline_segments'].shape[0] > 0: outfile.createDimension('num_of_outline_segments', mesh['outline_segments'].shape[0]) outfile.createVariable('outline_segments', netcdf_int, ('num_of_outline_segments', 'num_of_segment_ends')) outfile.variables['outline_segments'][:] = mesh['outline_segments'] if (mesh['outline_segment_tags'].shape[1] > 0): outfile.createDimension('num_of_outline_segment_tag_chars', mesh['outline_segment_tags'].shape[1]) outfile.createVariable('outline_segment_tags', netcdf_char, ('num_of_outline_segments', 'num_of_outline_segment_tag_chars')) outfile.variables['outline_segment_tags'][:] = \ mesh['outline_segment_tags'] # holes if (mesh['holes'].shape[0] > 0): outfile.createDimension('num_of_holes', mesh['holes'].shape[0]) outfile.createVariable('holes', netcdf_float, ('num_of_holes', 'num_of_dimensions')) outfile.variables['holes'][:] = mesh['holes'] # regions if (mesh['regions'].shape[0] > 0): outfile.createDimension('num_of_regions', mesh['regions'].shape[0]) outfile.createVariable('regions', netcdf_float, ('num_of_regions', 'num_of_dimensions')) outfile.createVariable('region_max_areas', netcdf_float, ('num_of_regions',)) outfile.variables['regions'][:] = mesh['regions'] outfile.variables['region_max_areas'][:] = mesh['region_max_areas'] if (mesh['region_tags'].shape[1] > 0): outfile.createDimension('num_of_region_tag_chars', mesh['region_tags'].shape[1]) outfile.createVariable('region_tags', netcdf_char, ('num_of_regions', 'num_of_region_tag_chars')) outfile.variables['region_tags'][:] = mesh['region_tags'] # geo_reference info if mesh.has_key('geo_reference') and not mesh['geo_reference'] == None: mesh['geo_reference'].write_NetCDF(outfile) outfile.close() ## # @brief Read a mesh file. # @param file_name Path to the file to read. # @return A dictionary containing the .msh file data. # @note Throws IOError if file not found. def _read_msh_file(file_name): """ Read in an msh file.""" #Check contents. Get NetCDF fd = open(file_name, 'r') fd.close() # throws prints to screen if file not present fid = NetCDFFile(file_name, netcdf_mode_r) mesh = {} # Get the variables - the triangulation try: mesh['vertices'] = fid.variables['vertices'][:] except KeyError: mesh['vertices'] = num.array([], num.int) #array default# try: mesh['vertex_attributes'] = fid.variables['vertex_attributes'][:] except KeyError: mesh['vertex_attributes'] = None mesh['vertex_attribute_titles'] = [] try: titles = fid.variables['vertex_attribute_titles'][:] for i, title in enumerate(titles): mesh['vertex_attribute_titles'].append(titles[i].tostring().strip()) except KeyError: pass try: mesh['segments'] = fid.variables['segments'][:] except KeyError: mesh['segments'] = num.array([], num.int) #array default# mesh['segment_tags'] =[] try: tags = fid.variables['segment_tags'][:] for i, tag in enumerate(tags): mesh['segment_tags'].append(tags[i].tostring().strip()) except KeyError: for ob in mesh['segments']: mesh['segment_tags'].append('') try: mesh['triangles'] = fid.variables['triangles'][:] mesh['triangle_neighbors'] = fid.variables['triangle_neighbors'][:] except KeyError: mesh['triangles'] = num.array([], num.int) #array default# mesh['triangle_neighbors'] = num.array([], num.int) #array default# mesh['triangle_tags'] = [] try: tags = fid.variables['triangle_tags'][:] for i, tag in enumerate(tags): mesh['triangle_tags'].append(tags[i].tostring().strip()) except KeyError: for ob in mesh['triangles']: mesh['triangle_tags'].append('') #the outline try: mesh['points'] = fid.variables['points'][:] except KeyError: mesh['points'] = [] try: mesh['point_attributes'] = fid.variables['point_attributes'][:] except KeyError: mesh['point_attributes'] = [] for point in mesh['points']: mesh['point_attributes'].append([]) try: mesh['outline_segments'] = fid.variables['outline_segments'][:] except KeyError: mesh['outline_segments'] = num.array([], num.int) #array default# mesh['outline_segment_tags'] =[] try: tags = fid.variables['outline_segment_tags'][:] for i, tag in enumerate(tags): mesh['outline_segment_tags'].append(tags[i].tostring().strip()) except KeyError: for ob in mesh['outline_segments']: mesh['outline_segment_tags'].append('') try: mesh['holes'] = fid.variables['holes'][:] except KeyError: mesh['holes'] = num.array([], num.int) #array default# try: mesh['regions'] = fid.variables['regions'][:] except KeyError: mesh['regions'] = num.array([], num.int) #array default# mesh['region_tags'] =[] try: tags = fid.variables['region_tags'][:] for i, tag in enumerate(tags): mesh['region_tags'].append(tags[i].tostring().strip()) except KeyError: for ob in mesh['regions']: mesh['region_tags'].append('') try: mesh['region_max_areas'] = fid.variables['region_max_areas'][:] except KeyError: mesh['region_max_areas'] = num.array([], num.int) #array default# try: geo_reference = Geo_reference(NetCDFObject=fid) mesh['geo_reference'] = geo_reference except AttributeError, e: #geo_ref not compulsory mesh['geo_reference'] = None fid.close() return mesh ## # @brief Export a boundary file. # @param file_name Path to the file to write. # @param points List of point index pairs. # @paran title Title string to write into file (first line). # @param delimiter Delimiter string to write between points. # @note Used by alpha shapes def export_boundary_file(file_name, points, title, delimiter=','): """Export a boundary file. Format: First line: Title variable Following lines: [point index][delimiter][point index] file_name - the name of the new file points - List of point index pairs [[p1, p2],[p3, p4]..] title - info to write in the first line """ fd = open(file_name, 'w') fd.write(title + '\n') # [point index][delimiter][point index] for point in points: fd.write(str(point[0]) + delimiter + str(point[1]) + '\n') fd.close() ################################################################################ # IMPORT/EXPORT POINTS FILES ################################################################################ ## # @brief # @param point_atts def extent_point_atts(point_atts): """Returns 4 points representing the extent This loses attribute info. """ point_atts['pointlist'] = extent(point_atts['pointlist']) point_atts['attributelist'] = {} return point_atts ## # @brief Get an extent array for a set of points. # @param points A set of points. # @return An extent array of form [[min_x, min_y], [max_x, min_y], # [max_x, max_y], [min_x, max_y]] def extent(points): points = num.array(points, num.float) max_x = min_x = points[0][0] max_y = min_y = points[0][1] for point in points[1:]: x = point[0] if x > max_x: max_x = x if x < min_x: min_x = x y = point[1] if y > max_y: max_y = y if y < min_y: min_y = y extent = num.array([[min_x, min_y], [max_x, min_y], [max_x, max_y], [min_x, max_y]]) return extent ## # @brief Reduce a points file until number of points is less than limit. # @param infile Path to input file to thin. # @param outfile Path to output thinned file. # @param max_points Max number of points in output file. # @param verbose True if this function is to be verbose. def reduce_pts(infile, outfile, max_points, verbose = False): """Reduce a points file until less than given size. Reduces a points file by removing every second point until the # of points is less than max_points. """ # check out pts2rectangular in least squares, and the use of reduction. # Maybe it does the same sort of thing? point_atts = _read_pts_file(infile) while point_atts['pointlist'].shape[0] > max_points: if verbose: print "point_atts['pointlist'].shape[0]" point_atts = half_pts(point_atts) export_points_file(outfile, point_atts) ## # @brief # @param infile # @param max_points # @param delimiter # @param verbose True if this function is to be verbose. # @return A list of def produce_half_point_files(infile, max_points, delimiter, verbose=False): point_atts = _read_pts_file(infile) root, ext = splitext(infile) outfiles = [] if verbose: print "# of points", point_atts['pointlist'].shape[0] while point_atts['pointlist'].shape[0] > max_points: point_atts = half_pts(point_atts) if verbose: print "# of points", point_atts['pointlist'].shape[0] outfile = root + delimiter + str(point_atts['pointlist'].shape[0]) + ext outfiles.append(outfile) export_points_file(outfile, point_atts) return outfiles ## # @brief # @param point_atts Object with attribute of 'pointlist'. # @return def point_atts2array(point_atts): # convert attribute list to array of floats point_atts['pointlist'] = num.array(point_atts['pointlist'], num.float) for key in point_atts['attributelist'].keys(): point_atts['attributelist'][key] = \ num.array(point_atts['attributelist'][key], num.float) return point_atts ## # @brief ?? # @param point_atts ?? # @return ?? def half_pts(point_atts): point_atts2array(point_atts) point_atts['pointlist'] = point_atts['pointlist'][::2] for key in point_atts['attributelist'].keys(): point_atts['attributelist'][key] = point_atts['attributelist'][key][::2] return point_atts ## # @brief ?? # @param dic ?? # @return ?? def concatinate_attributelist(dic): """ giving a dic[attribute title] = attribute return list of attribute titles, array of attributes """ point_attributes = num.array([], num.float) keys = dic.keys() key = keys.pop(0) point_attributes = num.reshape(dic[key], (dic[key].shape[0], 1)) for key in keys: reshaped = num.reshape(dic[key], (dic[key].shape[0], 1)) point_attributes = num.concatenate([point_attributes, reshaped], axis=1) return dic.keys(), point_attributes ## # @brief # @param dict # @param indices_to_keep # @return # @note FIXME(dsg), turn this dict plus methods into a class? def take_points(dict, indices_to_keep): dict = point_atts2array(dict) dict['pointlist'] = num.take(dict['pointlist'], indices_to_keep) for key in dict['attributelist'].keys(): dict['attributelist'][key] = num.take(dict['attributelist'][key], indices_to_keep) return dict ## # @brief ?? # @param dict1 ?? # @param dict2 ?? # @return ?? def add_point_dictionaries(dict1, dict2): """ """ dict1 = point_atts2array(dict1) dict2 = point_atts2array(dict2) combined = {} combined['pointlist'] = num.concatenate((dict2['pointlist'], dict1['pointlist']), axis=0) atts = {} for key in dict2['attributelist'].keys(): atts[key]= num.concatenate((dict2['attributelist'][key], dict1['attributelist'][key]), axis=0) combined['attributelist'] = atts combined['geo_reference'] = dict1['geo_reference'] return combined if __name__ == "__main__": pass