Source code for discretize.base.mesh_io

import os
import json
import numpy as np

from discretize.utils import mkvc
from discretize.base.base_mesh import BaseMesh
from discretize.utils.code_utils import deprecate_method
import warnings

try:
    from ..mixins import InterfaceTensorread_vtk
except ImportError:
    InterfaceTensorread_vtk = object


[docs]def load_mesh(file_name): """ Open a json file and load the mesh into the target class As long as there are no namespace conflicts, the target __class__ will be stored on the properties.HasProperties registry and may be fetched from there. :param str file_name: name of file to read in """ with open(file_name, "r") as outfile: jsondict = json.load(outfile) data = BaseMesh.deserialize(jsondict, trusted=True) return data
[docs]class TensorMeshIO(InterfaceTensorread_vtk): @classmethod def _readUBC_3DMesh(TensorMesh, file_name): """Read UBC GIF 3D tensor mesh and generate same dimension TensorMesh. Input: :param string file_name: path to the UBC GIF mesh file Output: :rtype: TensorMesh :return: The tensor mesh for the file_name. """ # Interal function to read cell size lines for the UBC mesh files. def readCellLine(line): line_list = [] for seg in line.split(): if "*" in seg: sp = seg.split("*") seg_arr = np.ones((int(sp[0]),)) * float(sp[1]) else: seg_arr = np.array([float(seg)], float) line_list.append(seg_arr) return np.concatenate(line_list) # Read the file as line strings, remove lines with comment = ! msh = np.genfromtxt(file_name, delimiter="\n", dtype=np.str, comments="!") # Fist line is the size of the model sizeM = np.array(msh[0].split(), dtype=float) # Second line is the South-West-Top corner coordinates. origin = np.array(msh[1].split(), dtype=float) # Read the cell sizes h1 = readCellLine(msh[2]) h2 = readCellLine(msh[3]) h3temp = readCellLine(msh[4]) # Invert the indexing of the vector to start from the bottom. h3 = h3temp[::-1] # Adjust the reference point to the bottom south west corner origin[2] = origin[2] - np.sum(h3) # Make the mesh tensMsh = TensorMesh([h1, h2, h3], origin=origin) return tensMsh @classmethod def _readUBC_2DMesh(TensorMesh, file_name): """Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg Input: :param string file_name: path to the UBC GIF mesh file Output: :rtype: TensorMesh :return: SimPEG TensorMesh 2D object """ fopen = open(file_name, "r") # Read down the file and unpack dx vector def unpackdx(fid, nrows): for ii in range(nrows): line = fid.readline() var = np.array(line.split(), dtype=float) if ii == 0: x0 = var[0] xvec = np.ones(int(var[2])) * (var[1] - var[0]) / int(var[2]) xend = var[1] else: xvec = np.hstack( (xvec, np.ones(int(var[1])) * (var[0] - xend) / int(var[1])) ) xend = var[0] return x0, xvec # Start with dx block # First line specifies the number of rows for x-cells line = fopen.readline() # Strip comments lines while line.startswith("!"): line = fopen.readline() nl = np.array(line.split(), dtype=int) [x0, dx] = unpackdx(fopen, nl[0]) # Move down the file until reaching the z-block line = fopen.readline() if not line: line = fopen.readline() # End with dz block # First line specifies the number of rows for z-cells line = fopen.readline() nl = np.array(line.split(), dtype=int) [z0, dz] = unpackdx(fopen, nl[0]) # Flip z0 to be the bottom of the mesh for SimPEG z0 = -(z0 + sum(dz)) dz = dz[::-1] # Make the mesh tensMsh = TensorMesh([dx, dz], origin=(x0, z0)) fopen.close() return tensMsh
[docs] @classmethod def read_UBC(TensorMesh, file_name, directory=""): """Wrapper to Read UBC GIF 2D and 3D tensor mesh and generate same dimension TensorMesh. Input: :param str file_name: path to the UBC GIF mesh file or just its name if directory is specified :param str directory: directory where the UBC GIF file lives Output: :rtype: TensorMesh :return: The tensor mesh for the file_name. """ # Check the expected mesh dimensions fname = os.path.join(directory, file_name) # Read the file as line strings, remove lines with comment = ! msh = np.genfromtxt( fname, delimiter="\n", dtype=np.str, comments="!", max_rows=1 ) # Fist line is the size of the model sizeM = np.array(msh.ravel()[0].split(), dtype=float) # Check if the mesh is a UBC 2D mesh if sizeM.shape[0] == 1: Tnsmsh = TensorMesh._readUBC_2DMesh(fname) # Check if the mesh is a UBC 3D mesh elif sizeM.shape[0] == 3: Tnsmsh = TensorMesh._readUBC_3DMesh(fname) else: raise Exception("File format not recognized") return Tnsmsh
def _readModelUBC_2D(mesh, file_name): """ Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg Input: :param string file_name: path to the UBC GIF 2D model file Output: :rtype: numpy.ndarray :return: model with TensorMesh ordered """ # Open fileand skip header... assume that we know the mesh already obsfile = np.genfromtxt(file_name, delimiter=" \n", dtype=np.str, comments="!") dim = tuple(np.array(obsfile[0].split(), dtype=int)) if mesh.shape_cells != dim: raise Exception("Dimension of the model and mesh mismatch") model = [] for line in obsfile[1:]: model.extend([float(val) for val in line.split()]) model = np.asarray(model) if not len(model) == mesh.nC: raise Exception( """Something is not right, expected size is {:d} but unwrap vector is size {:d}""".format( mesh.nC, len(model) ) ) return model.reshape(mesh.vnC, order="F")[:, ::-1].reshape(-1, order="F") def _readModelUBC_3D(mesh, file_name): """Read UBC 3DTensor mesh model and generate 3D Tensor mesh model Input: :param string file_name: path to the UBC GIF mesh file to read Output: :rtype: numpy.ndarray :return: model with TensorMesh ordered """ f = open(file_name, "r") model = np.array(list(map(float, f.readlines()))) f.close() nCx, nCy, nCz = mesh.shape_cells model = np.reshape(model, (nCz, nCx, nCy), order="F") model = model[::-1, :, :] model = np.transpose(model, (1, 2, 0)) model = mkvc(model) return model
[docs] def read_model_UBC(mesh, file_name, directory=""): """Read UBC 2D or 3D Tensor mesh model and generate Tensor mesh model Input: :param str file_name: path to the UBC GIF mesh file to read or just its name if directory is specified :param str directory: directory where the UBC GIF file lives Output: :rtype: numpy.ndarray :return: model with TensorMesh ordered """ fname = os.path.join(directory, file_name) if mesh.dim == 3: model = mesh._readModelUBC_3D(fname) elif mesh.dim == 2: model = mesh._readModelUBC_2D(fname) else: raise Exception("mesh must be a Tensor Mesh 2D or 3D") return model
[docs] def write_model_UBC(mesh, file_name, model, directory=""): """Writes a model associated with a TensorMesh to a UBC-GIF format model file. Input: :param str file_name: File to write to or just its name if directory is specified :param str directory: directory where the UBC GIF file lives :param numpy.ndarray model: The model """ fname = os.path.join(directory, file_name) if mesh.dim == 3: # Reshape model to a matrix modelMat = mesh.reshape(model, "CC", "CC", "M") # Transpose the axes modelMatT = modelMat.transpose((2, 0, 1)) # Flip z to positive down modelMatTR = mkvc(modelMatT[::-1, :, :]) np.savetxt(fname, modelMatTR.ravel()) elif mesh.dim == 2: modelMat = mesh.reshape(model, "CC", "CC", "M").T[::-1] f = open(fname, "w") f.write("{:d} {:d}\n".format(*mesh.shape_cells)) f.close() f = open(fname, "ab") np.savetxt(f, modelMat) f.close() else: raise Exception("mesh must be a Tensor Mesh 2D or 3D")
def _writeUBC_3DMesh(mesh, file_name, comment_lines=""): """Writes a TensorMesh to a UBC-GIF format mesh file. Input: :param string file_name: File to write to :param dict models: A dictionary of the models """ if not mesh.dim == 3: raise Exception("Mesh must be 3D") s = comment_lines s += "{0:d} {1:d} {2:d}\n".format(*tuple(mesh.vnC)) # Have to it in the same operation or use mesh.origin.copy(), # otherwise the mesh.origin is updated. origin = mesh.origin + np.array([0, 0, mesh.h[2].sum()]) nCx, nCy, nCz = mesh.shape_cells s += "{0:.6f} {1:.6f} {2:.6f}\n".format(*tuple(origin)) s += ("%.6f " * nCx + "\n") % tuple(mesh.h[0]) s += ("%.6f " * nCy + "\n") % tuple(mesh.h[1]) s += ("%.6f " * nCz + "\n") % tuple(mesh.h[2][::-1]) f = open(file_name, "w") f.write(s) f.close() def _writeUBC_2DMesh(mesh, file_name, comment_lines=""): """Writes a TensorMesh to a UBC-GIF format mesh file. Input: :param string file_name: File to write to :param dict models: A dictionary of the models """ if not mesh.dim == 2: raise Exception("Mesh must be 2D") def writeF(fx, outStr=""): # Init i = 0 origin = True x0 = fx[i] f = fx[i] number_segment = 0 auxStr = "" while True: i = i + 1 if i >= fx.size: break dx = -f + fx[i] f = fx[i] n = 1 for j in range(i + 1, fx.size): if -f + fx[j] == dx: n += 1 i += 1 f = fx[j] else: break number_segment += 1 if origin: auxStr += "{:.10f} {:.10f} {:d} \n".format(x0, f, n) origin = False else: auxStr += "{:.10f} {:d} \n".format(f, n) auxStr = "{:d}\n".format(number_segment) + auxStr outStr += auxStr return outStr # Grab face coordinates fx = mesh.nodes_x fz = -mesh.nodes_y[::-1] # Create the string outStr = comment_lines outStr = writeF(fx, outStr=outStr) outStr += "\n" outStr = writeF(fz, outStr=outStr) # Write file f = open(file_name, "w") f.write(outStr) f.close()
[docs] def write_UBC(mesh, file_name, models=None, directory="", comment_lines=""): """Writes a TensorMesh to a UBC-GIF format mesh file. Input: :param str file_name: File to write to :param str directory: directory where to save model :param dict models: A dictionary of the models :param str comment_lines: comment lines preceded with '!' to add """ fname = os.path.join(directory, file_name) if mesh.dim == 3: mesh._writeUBC_3DMesh(fname, comment_lines=comment_lines) elif mesh.dim == 2: mesh._writeUBC_2DMesh(fname, comment_lines=comment_lines) else: raise Exception("mesh must be a Tensor Mesh 2D or 3D") if models is None: return if not isinstance(models, dict): raise TypeError("models must be a dict") for key in models: if not isinstance(key, str): raise TypeError( "The dict key must be a string representing the file name" ) mesh.write_model_UBC(key, models[key], directory=directory)
# DEPRECATED
[docs] @classmethod def readUBC(TensorMesh, file_name, directory=""): warnings.warn( "TensorMesh.readUBC has been deprecated and will be removed in" "discretize 1.0.0. please use TensorMesh.read_UBC", DeprecationWarning, ) return TensorMesh.read_UBC(file_name, directory)
readModelUBC = deprecate_method( "read_model_UBC", "readModelUBC", removal_version="1.0.0" ) writeUBC = deprecate_method("write_UBC", "writeUBC", removal_version="1.0.0") writeModelUBC = deprecate_method( "write_model_UBC", "writeModelUBC", removal_version="1.0.0" )
[docs]class TreeMeshIO(object):
[docs] @classmethod def read_UBC(TreeMesh, meshFile, directory=""): """Read UBC 3D OcTree mesh file Input: :param str meshFile: path to the UBC GIF OcTree mesh file to read :rtype: discretize.TreeMesh :return: The octree mesh """ fname = os.path.join(directory, meshFile) fileLines = np.genfromtxt(fname, dtype=str, delimiter="\n", comments="!") nCunderMesh = np.array(fileLines[0].split("!")[0].split(), dtype=int) tswCorn = np.array(fileLines[1].split("!")[0].split(), dtype=float) smallCell = np.array(fileLines[2].split("!")[0].split(), dtype=float) # Read the index array indArr = np.genfromtxt( (line.encode("utf8") for line in fileLines[4::]), dtype=np.int ) nCunderMesh = nCunderMesh[:len(tswCorn)] # remove information related to core hs = [np.ones(nr) * sz for nr, sz in zip(nCunderMesh, smallCell)] origin = tswCorn origin[-1] -= np.sum(hs[-1]) ls = np.log2(nCunderMesh).astype(int) # if all ls are equal if min(ls) == max(ls): max_level = ls[0] else: max_level = min(ls) + 1 mesh = TreeMesh(hs, origin=origin) levels = indArr[:, -1] indArr = indArr[:, :-1] indArr -= 1 # shift by 1.... indArr = 2 * indArr + levels[:, None] # get cell center index indArr[:, -1] = 2 * nCunderMesh[-1] - indArr[:, -1] # switch direction of iz levels = max_level - np.log2(levels) # calculate level mesh.__setstate__((indArr, levels)) return mesh
[docs] def read_model_UBC(mesh, file_name): """Read UBC OcTree model and get vector :param string file_name: path to the UBC GIF model file to read :rtype: numpy.ndarray :return: OcTree model """ if type(file_name) is list: out = {} for f in file_name: out[f] = mesh.read_model_UBC(f) return out modArr = np.loadtxt(file_name) ubc_order = mesh._ubc_order # order_ubc will re-order from treemesh ordering to UBC ordering # need the opposite operation un_order = np.empty_like(ubc_order) un_order[ubc_order] = np.arange(len(ubc_order)) model = modArr[un_order].copy() # ensure a contiguous array return model
[docs] def write_UBC(mesh, file_name, models=None, directory=""): """Write UBC ocTree mesh and model files from a octree mesh and model. :param string file_name: File to write to :param dict models: Models in a dict, where each key is the file_name :param str directory: directory where to save model(s) """ uniform_hs = np.array([np.allclose(h, h[0]) for h in mesh.h]) if np.any(~uniform_hs): raise Exception("UBC form does not support variable cell widths") nCunderMesh = np.array([h.size for h in mesh.h], dtype=np.int64) tswCorn = mesh.origin.copy() tswCorn[-1] += np.sum(mesh.h[-1]) smallCell = np.array([h[0] for h in mesh.h]) nrCells = mesh.nC indArr, levels = mesh._ubc_indArr ubc_order = mesh._ubc_order indArr = indArr[ubc_order] levels = levels[ubc_order] # Write the UBC octree mesh file head = " ".join([f"{int(n)}" for n in nCunderMesh]) + " \n" head += " ".join([f"{v:.4f}" for v in tswCorn]) + " \n" head += " ".join([f"{v:.3f}" for v in smallCell]) + " \n" head += f"{int(nrCells)}" np.savetxt(file_name, np.c_[indArr, levels], fmt="%i", header=head, comments="") # Print the models if models is None: return if not isinstance(models, dict): raise TypeError("models must be a dict") for key in models: if not isinstance(key, str): raise TypeError( "The dict key must be a string representing the file name" ) mesh.write_model_UBC(key, models[key], directory=directory)
[docs] def write_model_UBC(mesh, file_name, model, directory=""): """Writes a model associated with a TreeMesh to a UBC-GIF format model file. Input: :param str file_name: File to write to or just its name if directory is specified :param str directory: directory where the UBC GIF file lives :param numpy.ndarray model: The model """ if type(file_name) is list: for f, m in zip(file_name, model): mesh.write_model_UBC(f, m) else: ubc_order = mesh._ubc_order fname = os.path.join(directory, file_name) m = model[ubc_order] np.savetxt(fname, m)
# DEPRECATED
[docs] @classmethod def readUBC(TreeMesh, file_name, directory=""): warnings.warn( "TensorMesh.readUBC has been deprecated and will be removed in" "discretize 1.0.0. please use TensorMesh.read_UBC", DeprecationWarning, ) return TreeMesh.read_UBC(file_name, directory)
readModelUBC = deprecate_method( "read_model_UBC", "readModelUBC", removal_version="1.0.0" ) writeUBC = deprecate_method("write_UBC", "writeUBC", removal_version="1.0.0") writeModelUBC = deprecate_method( "write_model_UBC", "writeModelUBC", removal_version="1.0.0" )