Source code for fatiando.mesher

"""
Generate and operate on various kinds of meshes and geometric elements

**Geometric elements**

* :class:`~fatiando.mesher.Polygon`
* :class:`~fatiando.mesher.Square`
* :class:`~fatiando.mesher.Prism`
* :class:`~fatiando.mesher.PolygonalPrism`
* :class:`~fatiando.mesher.Sphere`
* :class:`~fatiando.mesher.Tesseroid`

**Meshes**

* :class:`~fatiando.mesher.SquareMesh`
* :class:`~fatiando.mesher.PrismMesh`
* :class:`~fatiando.mesher.PrismRelief`
* :class:`~fatiando.mesher.TesseroidMesh`
* :class:`~fatiando.mesher.PointGrid`

**Utility functions**

* :func:`~fatiando.mesher.extract`: Extract the values of a physical
  property from the cells in a list
* :func:`~fatiando.mesher.vfilter`: Remove cells whose physical property
  value falls outside a given range
* :func:`~fatiando.mesher.vremove`: Remove the cells with a given physical
  property value

----

"""
from __future__ import division
import numpy
import scipy.special
import scipy.interpolate
import copy as cp

from . import gridder
from . import utils


[docs]class GeometricElement(object): """ Base class for all geometric elements. """ def __init__(self, props): self.props = {} if props is not None: for p in props: self.props[p] = props[p]
[docs] def addprop(self, prop, value): """ Add a physical property to this geometric element. If it already has the property, the given value will overwrite the existing one. Parameters: * prop : str Name of the physical property. * value : float The value of this physical property. """ self.props[prop] = value
[docs] def copy(self): """ Return a deep copy of the current instance.""" return cp.deepcopy(self)
[docs]class Polygon(GeometricElement): """ Create a polygon object. .. note:: Most applications require the vertices to be **clockwise**! Parameters: * vertices : list of lists List of [x, y] pairs with the coordinates of the vertices. * props : dict Physical properties assigned to the polygon. Ex: ``props={'density':10, 'susceptibility':10000}`` Examples:: >>> poly = Polygon([[0, 0], [1, 4], [2, 5]], {'density': 500}) >>> poly.props {'density': 500} >>> poly.nverts 3 >>> poly.vertices array([[0, 0], [1, 4], [2, 5]]) >>> poly.x array([0, 1, 2]) >>> poly.y array([0, 4, 5]) """ def __init__(self, vertices, props=None): super(Polygon, self).__init__(props) self._vertices = numpy.asarray(vertices) @property def vertices(self): return self._vertices @property def nverts(self): return len(self.vertices) @property def x(self): return self.vertices[:, 0] @property def y(self): return self.vertices[:, 1]
[docs]class Square(Polygon): """ Create a square object. Parameters: * bounds : list = [x1, x2, y1, y2] Coordinates of the top right and bottom left corners of the square * props : dict Physical properties assigned to the square. Ex: ``props={'density':10, 'slowness':10000}`` Example:: >>> sq = Square([0, 1, 2, 4], {'density': 750}) >>> sq.bounds [0, 1, 2, 4] >>> sq.x1 0 >>> sq.x2 1 >>> sq.props {'density': 750} >>> sq.addprop('magnetization', 100) >>> sq.props['magnetization'] 100 A square can be used as a :class:`~fatiando.mesher.Polygon`:: >>> sq.vertices array([[0, 2], [1, 2], [1, 4], [0, 4]]) >>> sq.x array([0, 1, 1, 0]) >>> sq.y array([2, 2, 4, 4]) >>> sq.nverts 4 """ def __init__(self, bounds, props=None): super(Square, self).__init__(None, props) self.x1, self.x2, self.y1, self.y2 = bounds @property def bounds(self): """ The x, y boundaries of the square as [xmin, xmax, ymin, ymax] """ return [self.x1, self.x2, self.y1, self.y2] @property def vertices(self): """ The vertices of the square. """ verts = numpy.array( [[self.x1, self.y1], [self.x2, self.y1], [self.x2, self.y2], [self.x1, self.y2]]) return verts def __str__(self): """Return a string representation of the square.""" names = [('x1', self.x1), ('x2', self.x2), ('y1', self.y1), ('y2', self.y2)] names.extend((p, self.props[p]) for p in sorted(self.props)) return ' | '.join('%s:%g' % (n, v) for n, v in names)
[docs]class SquareMesh(object): """ Generate a 2D regular mesh of squares. For all purposes, :class:`~fatiando.mesher.SquareMesh` can be used as a list of :class:`~fatiando.mesher.Square`. The order of the squares in the list is: x directions varies first, then y. Parameters: * bounds : list = [x1, x2, y1, y2] Boundaries of the mesh * shape : tuple = (ny, nx) Number of squares in the y and x dimension, respectively * props : dict Physical properties of each square in the mesh. Each key should be the name of a physical property. The corresponding value should be a list with the values of that particular property on each square of the mesh. Examples: >>> mesh = SquareMesh((0, 4, 0, 6), (2, 2)) >>> for s in mesh: ... print s x1:0 | x2:2 | y1:0 | y2:3 x1:2 | x2:4 | y1:0 | y2:3 x1:0 | x2:2 | y1:3 | y2:6 x1:2 | x2:4 | y1:3 | y2:6 >>> print mesh[1] x1:2 | x2:4 | y1:0 | y2:3 >>> print mesh[-1] x1:2 | x2:4 | y1:3 | y2:6 With physical properties:: >>> mesh = SquareMesh((0, 4, 0, 6), (2, 1), {'slowness':[3.4, 8.6]}) >>> for s in mesh: ... print s x1:0 | x2:4 | y1:0 | y2:3 | slowness:3.4 x1:0 | x2:4 | y1:3 | y2:6 | slowness:8.6 Or:: >>> mesh = SquareMesh((0, 4, 0, 6), (2, 1)) >>> mesh.addprop('slowness', [3.4, 8.6]) >>> for s in mesh: ... print s x1:0 | x2:4 | y1:0 | y2:3 | slowness:3.4 x1:0 | x2:4 | y1:3 | y2:6 | slowness:8.6 """ def __init__(self, bounds, shape, props=None): object.__init__(self) ny, nx = shape size = int(nx * ny) x1, x2, y1, y2 = bounds dx = (x2 - x1)/nx dy = (y2 - y1)/ny self.bounds = bounds self.shape = tuple(int(i) for i in shape) self.size = size self.dims = (dx, dy) # props has to be None, not {} by default because {} would be permanent # for all instaces of the class (like a class variable) and changes # to one instace would lead to changes in another (and a huge mess) if props is None: self.props = {} else: self.props = props # The index of the current square in an iteration. Needed when mesh is # used as an iterator self.i = 0 # List of masked squares. Will return None if trying to access them self.mask = [] def __len__(self): return self.size def __getitem__(self, index): # To walk backwards in the list if index < 0: index = self.size + index if index in self.mask: return None ny, nx = self.shape j = index//nx i = index - j*nx x1 = self.bounds[0] + self.dims[0] * i x2 = x1 + self.dims[0] y1 = self.bounds[2] + self.dims[1] * j y2 = y1 + self.dims[1] props = dict([p, self.props[p][index]] for p in self.props) return Square((x1, x2, y1, y2), props=props) def __iter__(self): self.i = 0 return self def next(self): if self.i >= self.size: raise StopIteration square = self.__getitem__(self.i) self.i += 1 return square
[docs] def addprop(self, prop, values): """ Add physical property values to the cells in the mesh. Different physical properties of the mesh are stored in a dictionary. Parameters: * prop : str Name of the physical property * values : list or array The value of this physical property in each square of the mesh. For the ordering of squares in the mesh see :class:`~fatiando.mesher.SquareMesh` """ self.props[prop] = values
[docs] def img2prop(self, fname, vmin, vmax, prop): """ Load the physical property value from an image file. The image is converted to gray scale and the gray intensity of each pixel is used to set the value of the physical property of the cells in the mesh. Gray intensity values are scaled to the range ``[vmin, vmax]``. If the shape of image (number of pixels in y and x) is different from the shape of the mesh, the image will be interpolated to match the shape of the mesh. Parameters: * fname : str Name of the image file * vmax, vmin : float Range of physical property values (used to convert the gray scale to physical property values) * prop : str Name of the physical property """ self.props[prop] = utils.fromimage(fname, ranges=[vmin, vmax], shape=self.shape)[::-1, :].ravel()
[docs] def get_xs(self): """ Get a list of the x coordinates of the corners of the cells in the mesh. If the mesh has nx cells, get_xs() will return nx + 1 values. """ dx, dy = self.dims x1, x2, y1, y2 = self.bounds ny, nx = self.shape xs = numpy.arange(x1, x2 + dx, dx, 'f') if len(xs) == nx + 2: return xs[0:-1] elif len(xs) == nx: xs = xs.tolist() xs.append(x2) return numpy.array(xs) else: return xs
[docs] def get_ys(self): """ Get a list of the y coordinates of the corners of the cells in the mesh. If the mesh has ny cells, get_ys() will return ny + 1 values. """ dx, dy = self.dims x1, x2, y1, y2 = self.bounds ny, nx = self.shape ys = numpy.arange(y1, y2, dy, 'f') if len(ys) == ny + 2: return ys[0:-1] elif len(ys) == ny: ys = ys.tolist() ys.append(y2) return numpy.array(ys) else: return ys
[docs] def copy(self): """ Return a deep copy of the current instance.""" return cp.deepcopy(self)
[docs]class Prism(GeometricElement): """ Create a 3D right rectangular prism. .. note:: The coordinate system used is x -> North, y -> East and z -> Down Parameters: * x1, x2 : float South and north borders of the prism * y1, y2 : float West and east borders of the prism * z1, z2 : float Top and bottom of the prism * props : dict Physical properties assigned to the prism. Ex: ``props={'density':10, 'magnetization':10000}`` Examples: >>> from fatiando.mesher import Prism >>> p = Prism(1, 2, 3, 4, 5, 6, {'density':200}) >>> p.props['density'] 200 >>> print p.get_bounds() [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] >>> print p x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6 | density:200 >>> p = Prism(1, 2, 3, 4, 5, 6) >>> print p x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6 >>> p.addprop('density', 2670) >>> print p x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6 | density:2670 """ def __init__(self, x1, x2, y1, y2, z1, z2, props=None): GeometricElement.__init__(self, props) self.x1 = float(x1) self.x2 = float(x2) self.y1 = float(y1) self.y2 = float(y2) self.z1 = float(z1) self.z2 = float(z2) def __str__(self): """Return a string representation of the prism.""" names = [('x1', self.x1), ('x2', self.x2), ('y1', self.y1), ('y2', self.y2), ('z1', self.z1), ('z2', self.z2)] names.extend((p, self.props[p]) for p in sorted(self.props)) return ' | '.join('%s:%g' % (n, v) for n, v in names)
[docs] def get_bounds(self): """ Get the bounding box of the prism (i.e., the borders of the prism). Returns: * bounds : list ``[x1, x2, y1, y2, z1, z2]``, the bounds of the prism Examples: >>> p = Prism(1, 2, 3, 4, 5, 6) >>> print p.get_bounds() [1.0, 2.0, 3.0, 4.0, 5.0, 6.0] """ return [self.x1, self.x2, self.y1, self.y2, self.z1, self.z2]
[docs] def center(self): """ Return the coordinates of the center of the prism. Returns: * coords : list = [xc, yc, zc] Coordinates of the center Example: >>> prism = Prism(1, 2, 1, 3, 0, 2) >>> print prism.center() [ 1.5 2. 1. ] """ xc = 0.5 * (self.x1 + self.x2) yc = 0.5 * (self.y1 + self.y2) zc = 0.5 * (self.z1 + self.z2) return numpy.array([xc, yc, zc])
[docs]class Tesseroid(GeometricElement): """ Create a tesseroid (spherical prism). Parameters: * w, e : float West and east borders of the tesseroid in decimal degrees * s, n : float South and north borders of the tesseroid in decimal degrees * top, bottom : float Bottom and top of the tesseroid with respect to the mean earth radius in meters. Ex: if the top is 100 meters above the mean earth radius, ``top=100``, if 100 meters below ``top=-100``. * props : dict Physical properties assigned to the tesseroid. Ex: ``props={'density':10, 'magnetization':10000}`` Examples: >>> from fatiando.mesher import Tesseroid >>> t = Tesseroid(1, 2, 3, 4, 6, 5, {'density':200}) >>> t.props['density'] 200 >>> print t.get_bounds() [1.0, 2.0, 3.0, 4.0, 6.0, 5.0] >>> print t w:1 | e:2 | s:3 | n:4 | top:6 | bottom:5 | density:200 >>> t = Tesseroid(1, 2, 3, 4, 6, 5) >>> print t w:1 | e:2 | s:3 | n:4 | top:6 | bottom:5 >>> t.addprop('density', 2670) >>> print t w:1 | e:2 | s:3 | n:4 | top:6 | bottom:5 | density:2670 """ def __init__(self, w, e, s, n, top, bottom, props=None): GeometricElement.__init__(self, props) self.w = float(w) self.e = float(e) self.s = float(s) self.n = float(n) self.bottom = float(bottom) self.top = float(top) def __str__(self): """Return a string representation of the tesseroid.""" names = [('w', self.w), ('e', self.e), ('s', self.s), ('n', self.n), ('top', self.top), ('bottom', self.bottom)] names.extend((p, self.props[p]) for p in sorted(self.props)) return ' | '.join('%s:%g' % (n, v) for n, v in names)
[docs] def get_bounds(self): """ Get the bounding box of the tesseroid (i.e., the borders). Returns: * bounds : list ``[w, e, s, n, top, bottom]``, the bounds of the tesseroid Examples: >>> t = Tesseroid(1, 2, 3, 4, 6, 5) >>> print t.get_bounds() [1.0, 2.0, 3.0, 4.0, 6.0, 5.0] """ return [self.w, self.e, self.s, self.n, self.top, self.bottom]
[docs] def half(self, lon=True, lat=True, r=True): """ Divide the tesseroid in 2 halfs for each dimension (total 8) The smaller tesseroids will share the large one's props. Parameters: * lon, lat, r : True or False Dimensions along which the tesseroid will be split in half. Returns: * tesseroids : list A list of maximum 8 tesseroids that make up the larger one. Examples:: >>> tess = Tesseroid(-10, 10, -20, 20, 0, -40, {'density':2}) >>> split = tess.half() >>> print len(split) 8 >>> for t in split: ... print t w:-10 | e:0 | s:-20 | n:0 | top:-20 | bottom:-40 | density:2 w:-10 | e:0 | s:-20 | n:0 | top:0 | bottom:-20 | density:2 w:-10 | e:0 | s:0 | n:20 | top:-20 | bottom:-40 | density:2 w:-10 | e:0 | s:0 | n:20 | top:0 | bottom:-20 | density:2 w:0 | e:10 | s:-20 | n:0 | top:-20 | bottom:-40 | density:2 w:0 | e:10 | s:-20 | n:0 | top:0 | bottom:-20 | density:2 w:0 | e:10 | s:0 | n:20 | top:-20 | bottom:-40 | density:2 w:0 | e:10 | s:0 | n:20 | top:0 | bottom:-20 | density:2 >>> tess = Tesseroid(-15, 15, -20, 20, 0, -40) >>> split = tess.half(lat=False) >>> print len(split) 4 >>> for t in split: ... print t w:-15 | e:0 | s:-20 | n:20 | top:-20 | bottom:-40 w:-15 | e:0 | s:-20 | n:20 | top:0 | bottom:-20 w:0 | e:15 | s:-20 | n:20 | top:-20 | bottom:-40 w:0 | e:15 | s:-20 | n:20 | top:0 | bottom:-20 """ dlon = 0.5 * (self.e - self.w) dlat = 0.5 * (self.n - self.s) dh = 0.5 * (self.top - self.bottom) wests = [self.w, self.w + dlon] souths = [self.s, self.s + dlat] bottoms = [self.bottom, self.bottom + dh] if not lon: dlon *= 2 wests.pop() if not lat: dlat *= 2 souths.pop() if not r: dh *= 2 bottoms.pop() split = [ Tesseroid(i, i + dlon, j, j + dlat, k + dh, k, props=self.props) for i in wests for j in souths for k in bottoms] return split
[docs] def split(self, nlon, nlat, nh): """ Split the tesseroid into smaller ones. The smaller tesseroids will share the large one's props. Parameters: * nlon, nlat, nh : int The number of sections to split in the longitudinal, latitudinal, and vertical dimensions Returns: * tesseroids : list A list of nlon*nlat*nh tesseroids that make up the larger one. Examples:: >>> tess = Tesseroid(-10, 10, -20, 20, 0, -40, {'density':2}) >>> split = tess.split(1, 2, 2) >>> print len(split) 4 >>> for t in split: ... print t w:-10 | e:10 | s:-20 | n:0 | top:-20 | bottom:-40 | density:2 w:-10 | e:10 | s:-20 | n:0 | top:0 | bottom:-20 | density:2 w:-10 | e:10 | s:0 | n:20 | top:-20 | bottom:-40 | density:2 w:-10 | e:10 | s:0 | n:20 | top:0 | bottom:-20 | density:2 >>> tess = Tesseroid(-15, 15, -20, 20, 0, -40) >>> split = tess.split(3, 1, 1) >>> print len(split) 3 >>> for t in split: ... print t w:-15 | e:-5 | s:-20 | n:20 | top:0 | bottom:-40 w:-5 | e:5 | s:-20 | n:20 | top:0 | bottom:-40 w:5 | e:15 | s:-20 | n:20 | top:0 | bottom:-40 """ wests = numpy.linspace(self.w, self.e, nlon + 1) souths = numpy.linspace(self.s, self.n, nlat + 1) bottoms = numpy.linspace(self.bottom, self.top, nh + 1) dlon = wests[1] - wests[0] dlat = souths[1] - souths[0] dh = bottoms[1] - bottoms[0] tesseroids = [ Tesseroid(i, i + dlon, j, j + dlat, k + dh, k, props=self.props) for i in wests[:-1] for j in souths[:-1] for k in bottoms[:-1]] return tesseroids
[docs]class Sphere(GeometricElement): """ Create a sphere. .. note:: The coordinate system used is x -> North, y -> East and z -> Down Parameters: * x, y, z : float The coordinates of the center of the sphere * radius : float The radius of the sphere * props : dict Physical properties assigned to the prism. Ex: ``props={'density':10, 'magnetization':10000}`` Examples: >>> s = Sphere(1, 2, 3, 10, {'magnetization':200}) >>> s.props['magnetization'] 200 >>> s.addprop('density', 20) >>> print s.props['density'] 20 >>> print s x:1 | y:2 | z:3 | radius:10 | density:20 | magnetization:200 >>> s = Sphere(1, 2, 3, 4) >>> print s x:1 | y:2 | z:3 | radius:4 >>> s.addprop('density', 2670) >>> print s x:1 | y:2 | z:3 | radius:4 | density:2670 """ def __init__(self, x, y, z, radius, props=None): GeometricElement.__init__(self, props) self.x = float(x) self.y = float(y) self.z = float(z) self.radius = float(radius) self.center = numpy.array([x, y, z]) def __str__(self): """Return a string representation of the sphere.""" names = [('x', self.x), ('y', self.y), ('z', self.z), ('radius', self.radius)] names.extend((p, self.props[p]) for p in sorted(self.props)) return ' | '.join('%s:%g' % (n, v) for n, v in names)
[docs]class PolygonalPrism(GeometricElement): """ Create a 3D prism with polygonal crossection. .. note:: The coordinate system used is x -> North, y -> East and z -> Down .. note:: *vertices* must be **CLOCKWISE** or will give inverse result. Parameters: * vertices : list of lists Coordinates of the vertices. A list of ``[x, y]`` pairs. * z1, z2 : float Top and bottom of the prism * props : dict Physical properties assigned to the prism. Ex: ``props={'density':10, 'magnetization':10000}`` Examples: >>> verts = [[1, 1], [1, 2], [2, 2], [2, 1]] >>> p = PolygonalPrism(verts, 0, 3, props={'temperature':25}) >>> p.props['temperature'] 25 >>> print p.x [ 1. 1. 2. 2.] >>> print p.y [ 1. 2. 2. 1.] >>> print p.z1, p.z2 0.0 3.0 >>> p.addprop('density', 2670) >>> print p.props['density'] 2670 """ def __init__(self, vertices, z1, z2, props=None): GeometricElement.__init__(self, props) self.x = numpy.fromiter((v[0] for v in vertices), dtype=numpy.float) self.y = numpy.fromiter((v[1] for v in vertices), dtype=numpy.float) self.z1 = float(z1) self.z2 = float(z2) self.nverts = len(vertices)
[docs] def topolygon(self): """ Get the polygon describing the prism viewed from above. Returns: * polygon : :func:`fatiando.mesher.Polygon` The polygon Example: >>> verts = [[1, 1], [1, 2], [2, 2], [2, 1]] >>> p = PolygonalPrism(verts, 0, 100) >>> poly = p.topolygon() >>> print poly.x [ 1. 1. 2. 2.] >>> print poly.y [ 1. 2. 2. 1.] """ verts = numpy.transpose([self.x, self.y]) return Polygon(verts, self.props)
[docs]class PointGrid(object): """ Create a regular grid of 3D point sources (spheres of unit volume). Use this as a 1D list of :class:`~fatiando.mesher.Sphere`. Grid points are ordered like a C matrix, first each row in a column, then change columns. In this case, the x direction (North-South) are the rows and y (East-West) are the columns. Parameters: * area : list = [x1, x2, y1, y2] The area where the grid will be spread out * z : float or 1d-array The z coordinates of each point in the grid (remember, z is positive downward). * shape : tuple = (nx, ny) The number of points in the x and y directions * props : dict Physical properties of each point in the grid. Each key should be the name of a physical property. The corresponding value should be a list with the values of that particular property for each point in the grid. Examples:: >>> g = PointGrid([0, 10, 2, 6], 200, (2, 3)) >>> g.shape (2, 3) >>> g.size 6 >>> g[0].center array([ 0., 2., 200.]) >>> g[-1].center array([ 10., 6., 200.]) >>> for p in g: ... p.center array([ 0., 2., 200.]) array([ 0., 4., 200.]) array([ 0., 6., 200.]) array([ 10., 2., 200.]) array([ 10., 4., 200.]) array([ 10., 6., 200.]) >>> g.x.reshape(g.shape) array([[ 0., 0., 0.], [ 10., 10., 10.]]) >>> g.y.reshape(g.shape) array([[ 2., 4., 6.], [ 2., 4., 6.]]) >>> g.dx, g.dy (10.0, 2.0) """ def __init__(self, area, z, shape, props=None): object.__init__(self) self.area = area self.shape = shape if props is None: self.props = {} else: self.props = props nx, ny = shape self.size = nx*ny self.z = numpy.zeros(self.size) + z self.radius = scipy.special.cbrt(3. / (4. * numpy.pi)) self.x, self.y = gridder.regular(area, shape) # The spacing between points self.dx, self.dy = gridder.spacing(area, shape) def __len__(self): return self.size def __getitem__(self, index): if not isinstance(index, int): raise IndexError('Invalid index type. Should be int.') if index >= self.size or index < -self.size: raise IndexError('Grid index out of range.') # To walk backwards in the list if index < 0: index = self.size + index props = dict([p, self.props[p][index]] for p in self.props) sphere = Sphere(self.x[index], self.y[index], self.z[index], self.radius, props=props) return sphere def __iter__(self): self.i = 0 return self def next(self): if self.i >= self.size: raise StopIteration sphere = self.__getitem__(self.i) self.i += 1 return sphere
[docs] def addprop(self, prop, values): """ Add physical property values to the points in the grid. Different physical properties of the grid are stored in a dictionary. Parameters: * prop : str Name of the physical property. * values : list or array Value of this physical property in each point of the grid """ self.props[prop] = values
[docs] def split(self, shape): """ Divide the grid into subgrids. .. note:: Remember that x is the North-South direction and y is East-West. Parameters: * shape : tuple = (nx, ny) Number of subgrids along the x and y directions, respectively. Returns: * subgrids : list List of :class:`~fatiando.mesher.PointGrid` Examples:: >>> import numpy >>> z = numpy.linspace(0, 1100, 12) >>> g = PointGrid((0, 3, 0, 2), z, (4, 3)) >>> g.addprop('bla', [1, 2, 3, ... 4, 5, 6, ... 7, 8, 9, ... 10, 11, 12]) >>> grids = g.split((2, 3)) >>> for s in grids: ... s.props['bla'] array([1, 4]) array([2, 5]) array([3, 6]) array([ 7, 10]) array([ 8, 11]) array([ 9, 12]) >>> for s in grids: ... s.x array([ 0., 1.]) array([ 0., 1.]) array([ 0., 1.]) array([ 2., 3.]) array([ 2., 3.]) array([ 2., 3.]) >>> for s in grids: ... s.y array([ 0., 0.]) array([ 1., 1.]) array([ 2., 2.]) array([ 0., 0.]) array([ 1., 1.]) array([ 2., 2.]) >>> for s in grids: ... s.z array([ 0., 300.]) array([ 100., 400.]) array([ 200., 500.]) array([ 600., 900.]) array([ 700., 1000.]) array([ 800., 1100.]) """ nx, ny = shape totalx, totaly = self.shape if totalx % nx != 0 or totaly % ny != 0: raise ValueError( 'Cannot split! nx and ny must be divisible by grid shape') x1, x2, y1, y2 = self.area xs = numpy.linspace(x1, x2, totalx) ys = numpy.linspace(y1, y2, totaly) mx, my = (totalx//nx, totaly//ny) dx, dy = self.dx*(mx - 1), self.dy*(my - 1) subs = [] for i, xstart in enumerate(xs[::mx]): for j, ystart in enumerate(ys[::my]): area = [xstart, xstart + dx, ystart, ystart + dy] props = {} for p in self.props: pmatrix = numpy.reshape(self.props[p], self.shape) props[p] = pmatrix[i*mx:(i + 1)*mx, j*my:(j + 1)*my].ravel() zmatrix = numpy.reshape(self.z, self.shape) zs = zmatrix[i*mx:(i + 1)*mx, j*my:(j + 1)*my].ravel() subs.append(PointGrid(area, zs, (mx, my), props)) return subs
[docs] def copy(self): """ Return a deep copy of the current instance.""" return cp.deepcopy(self)
[docs]class PrismRelief(object): """ Generate a 3D model of a relief (topography) using prisms. Use to generate: * topographic model * basin model * Moho model * etc PrismRelief can used as list of prisms. It acts as an iteratior (so you can loop over prisms). It also has a ``__getitem__`` method to access individual elements in the mesh. In practice, PrismRelief should be able to be passed to any function that asks for a list of prisms, like :func:`fatiando.gravmag.prism.gz`. Parameters: * ref : float Reference level. Prisms will have: * bottom on zref and top on z if z > zref; * bottom on z and top on zref otherwise. * dims : tuple = (dy, dx) Dimensions of the prisms in the y and x directions * nodes : list of lists = [x, y, z] Coordinates of the center of the top face of each prism.x, y, and z are lists with the x, y and z coordinates on a regular grid. """ def __init__(self, ref, dims, nodes): object.__init__(self) x, y, z = nodes if len(x) != len(y) != len(z): raise ValueError( "nodes has x, y, z coordinate arrays of different lengths") self.x, self.y, self.z = x, y, z self.size = len(x) self.ref = ref self.dy, self.dx = dims self.props = {} # The index of the current prism in an iteration. Needed when mesh is # used as an iterator self.i = 0 def __len__(self): return self.size def __iter__(self): self.i = 0 return self def __getitem__(self, index): # To walk backwards in the list if index < 0: index = self.size + index xc, yc, zc = self.x[index], self.y[index], self.z[index] x1 = xc - 0.5 * self.dx x2 = xc + 0.5 * self.dx y1 = yc - 0.5 * self.dy y2 = yc + 0.5 * self.dy if zc <= self.ref: z1 = zc z2 = self.ref else: z1 = self.ref z2 = zc props = dict([p, self.props[p][index]] for p in self.props) return Prism(x1, x2, y1, y2, z1, z2, props=props) def next(self): if self.i >= self.size: raise StopIteration prism = self.__getitem__(self.i) self.i += 1 return prism
[docs] def addprop(self, prop, values): """ Add physical property values to the prisms. .. warning:: If the z value of any point in the relief is below the reference level, its corresponding prism will have the physical property value with oposite sign than was assigned to it. Parameters: * prop : str Name of the physical property. * values : list List or array with the value of this physical property in each prism of the relief. """ def correct(v, i): if self.z[i] > self.ref: return -v return v self.props[prop] = [correct(v, i) for i, v in enumerate(values)]
[docs] def copy(self): """ Return a deep copy of the current instance.""" return cp.deepcopy(self)
[docs]class PrismMesh(object): """ Generate a 3D regular mesh of right rectangular prisms. Prisms are ordered as follows: first layers (z coordinate), then EW rows (y) and finaly x coordinate (NS). .. note:: Remember that the coordinate system is x->North, y->East and z->Down Ex: in a mesh with shape ``(3,3,3)`` the 15th element (index 14) has z index 1 (second layer), y index 1 (second row), and x index 2 (third element in the column). :class:`~fatiando.mesher.PrismMesh` can used as list of prisms. It acts as an iteratior (so you can loop over prisms). It also has a ``__getitem__`` method to access individual elements in the mesh. In practice, :class:`~fatiando.mesher.PrismMesh` should be able to be passed to any function that asks for a list of prisms, like :func:`fatiando.gravmag.prism.gz`. To make the mesh incorporate a topography, use :meth:`~fatiando.mesher.PrismMesh.carvetopo` Parameters: * bounds : list = [xmin, xmax, ymin, ymax, zmin, zmax] Boundaries of the mesh. * shape : tuple = (nz, ny, nx) Number of prisms in the x, y, and z directions. * props : dict Physical properties of each prism in the mesh. Each key should be the name of a physical property. The corresponding value should be a list with the values of that particular property on each prism of the mesh. Examples: >>> from fatiando.mesher import PrismMesh >>> mesh = PrismMesh((0, 1, 0, 2, 0, 3), (1, 2, 2)) >>> for p in mesh: ... print p x1:0 | x2:0.5 | y1:0 | y2:1 | z1:0 | z2:3 x1:0.5 | x2:1 | y1:0 | y2:1 | z1:0 | z2:3 x1:0 | x2:0.5 | y1:1 | y2:2 | z1:0 | z2:3 x1:0.5 | x2:1 | y1:1 | y2:2 | z1:0 | z2:3 >>> print mesh[0] x1:0 | x2:0.5 | y1:0 | y2:1 | z1:0 | z2:3 >>> print mesh[-1] x1:0.5 | x2:1 | y1:1 | y2:2 | z1:0 | z2:3 One with physical properties:: >>> props = {'density':[2670.0, 1000.0]} >>> mesh = PrismMesh((0, 2, 0, 4, 0, 3), (1, 1, 2), props=props) >>> for p in mesh: ... print p x1:0 | x2:1 | y1:0 | y2:4 | z1:0 | z2:3 | density:2670 x1:1 | x2:2 | y1:0 | y2:4 | z1:0 | z2:3 | density:1000 or equivalently:: >>> mesh = PrismMesh((0, 2, 0, 4, 0, 3), (1, 1, 2)) >>> mesh.addprop('density', [200, -1000.0]) >>> for p in mesh: ... print p x1:0 | x2:1 | y1:0 | y2:4 | z1:0 | z2:3 | density:200 x1:1 | x2:2 | y1:0 | y2:4 | z1:0 | z2:3 | density:-1000 You can use :meth:`~fatiando.mesher.PrismMesh.get_xs` (and similar methods for y and z) to get the x coordinates of the prisms in the mesh:: >>> mesh = PrismMesh((0, 2, 0, 4, 0, 3), (1, 1, 2)) >>> print mesh.get_xs() [ 0. 1. 2.] >>> print mesh.get_ys() [ 0. 4.] >>> print mesh.get_zs() [ 0. 3.] The ``shape`` of the mesh must be integer! >>> mesh = PrismMesh((0, 2, 0, 4, 0, 3), (1, 1, 2.5)) Traceback (most recent call last): ... AttributeError: Invalid mesh shape (1, 1, 2.5). shape must be integers """ celltype = Prism def __init__(self, bounds, shape, props=None): object.__init__(self) nz, ny, nx = shape if (not isinstance(nx, int) or not isinstance(ny, int) or not isinstance(nz, int)): raise AttributeError( 'Invalid mesh shape {}. shape must be integers'.format( str(shape))) size = int(nx * ny * nz) x1, x2, y1, y2, z1, z2 = bounds dx = (x2 - x1)/nx dy = (y2 - y1)/ny dz = (z2 - z1)/nz self.shape = tuple(int(i) for i in shape) self.size = size self.dims = (dx, dy, dz) self.bounds = bounds if props is None: self.props = {} else: self.props = props # The index of the current prism in an iteration. Needed when mesh is # used as an iterator self.i = 0 # List of masked prisms. Will return None if trying to access them self.mask = [] # Wether or not to change heights to z coordinate self.zdown = True def __len__(self): return self.size def __getitem__(self, index): if index >= self.size or index < -self.size: raise IndexError('mesh index out of range') # To walk backwards in the list if index < 0: index = self.size + index if index in self.mask: return None nz, ny, nx = self.shape k = index//(nx*ny) j = (index - k*(nx*ny))//nx i = (index - k*(nx*ny) - j*nx) x1 = self.bounds[0] + self.dims[0] * i x2 = x1 + self.dims[0] y1 = self.bounds[2] + self.dims[1] * j y2 = y1 + self.dims[1] z1 = self.bounds[4] + self.dims[2] * k z2 = z1 + self.dims[2] props = dict([p, self.props[p][index]] for p in self.props) return self.celltype(x1, x2, y1, y2, z1, z2, props=props) def __iter__(self): self.i = 0 return self def next(self): if self.i >= self.size: raise StopIteration prism = self.__getitem__(self.i) self.i += 1 return prism
[docs] def addprop(self, prop, values): """ Add physical property values to the cells in the mesh. Different physical properties of the mesh are stored in a dictionary. Parameters: * prop : str Name of the physical property. * values : list or array Value of this physical property in each prism of the mesh. For the ordering of prisms in the mesh see :class:`~fatiando.mesher.PrismMesh` """ self.props[prop] = values
[docs] def carvetopo(self, x, y, height): """ Mask (remove) prisms from the mesh that are above the topography. Accessing the ith prism will return None if it was masked (above the topography). Also mask prisms outside of the topography grid provided. The topography height information does not need to be on a regular grid, it will be interpolated. Parameters: * x, y : lists x and y coordinates of the grid points * height : list or array Array with the height of the topography """ nz, ny, nx = self.shape x1, x2, y1, y2, z1, z2 = self.bounds dx, dy, dz = self.dims # The coordinates of the centers of the cells xc = numpy.arange(x1, x2, dx) + 0.5 * dx # Sometimes arange returns more due to rounding if len(xc) > nx: xc = xc[:-1] yc = numpy.arange(y1, y2, dy) + 0.5 * dy if len(yc) > ny: yc = yc[:-1] zc = numpy.arange(z1, z2, dz) + 0.5 * dz if len(zc) > nz: zc = zc[:-1] XC, YC = numpy.meshgrid(xc, yc) topo = scipy.interpolate.griddata((x, y), height, (XC, YC), method='cubic').ravel() if self.zdown: # -1 if to transform height into z coordinate topo = -1 * topo # griddata returns a masked array. If the interpolated point is out of # of the data range, mask will be True. Use this to remove all cells # below a masked topo point (ie, one with no height information) if numpy.ma.isMA(topo): topo_mask = topo.mask else: topo_mask = [False for i in xrange(len(topo))] c = 0 for cellz in zc: for h, masked in zip(topo, topo_mask): if (masked or (cellz < h and self.zdown) or (cellz > h and not self.zdown)): self.mask.append(c) c += 1
[docs] def get_xs(self): """ Return an array with the x coordinates of the prisms in mesh. """ x1, x2, y1, y2, z1, z2 = self.bounds dx, dy, dz = self.dims nz, ny, nx = self.shape xs = numpy.arange(x1, x2 + dx, dx) if xs.size > nx + 1: return xs[:-1] return xs
[docs] def get_ys(self): """ Return an array with the y coordinates of the prisms in mesh. """ x1, x2, y1, y2, z1, z2 = self.bounds dx, dy, dz = self.dims nz, ny, nx = self.shape ys = numpy.arange(y1, y2 + dy, dy) if ys.size > ny + 1: return ys[:-1] return ys
[docs] def get_zs(self): """ Return an array with the z coordinates of the prisms in mesh. """ x1, x2, y1, y2, z1, z2 = self.bounds dx, dy, dz = self.dims nz, ny, nx = self.shape zs = numpy.arange(z1, z2 + dz, dz) if zs.size > nz + 1: return zs[:-1] return zs
[docs] def get_layer(self, i): """ Return the set of prisms corresponding to the ith layer of the mesh. Parameters: * i : int The index of the layer Returns: * prisms : list of :class:`~fatiando.mesher.Prism` The prisms in the ith layer Examples:: >>> mesh = PrismMesh((0, 2, 0, 2, 0, 2), (2, 2, 2)) >>> layer = mesh.get_layer(0) >>> for p in layer: ... print p x1:0 | x2:1 | y1:0 | y2:1 | z1:0 | z2:1 x1:1 | x2:2 | y1:0 | y2:1 | z1:0 | z2:1 x1:0 | x2:1 | y1:1 | y2:2 | z1:0 | z2:1 x1:1 | x2:2 | y1:1 | y2:2 | z1:0 | z2:1 >>> layer = mesh.get_layer(1) >>> for p in layer: ... print p x1:0 | x2:1 | y1:0 | y2:1 | z1:1 | z2:2 x1:1 | x2:2 | y1:0 | y2:1 | z1:1 | z2:2 x1:0 | x2:1 | y1:1 | y2:2 | z1:1 | z2:2 x1:1 | x2:2 | y1:1 | y2:2 | z1:1 | z2:2 """ nz, ny, nx = self.shape if i >= nz or i < 0: raise IndexError('Layer index %d is out of range.' % (i)) start = i * nx * ny end = (i + 1) * nx * ny layer = [self.__getitem__(p) for p in xrange(start, end)] return layer
[docs] def layers(self): """ Returns an iterator over the layers of the mesh. Examples:: >>> mesh = PrismMesh((0, 2, 0, 2, 0, 2), (2, 2, 2)) >>> for layer in mesh.layers(): ... for p in layer: ... print p x1:0 | x2:1 | y1:0 | y2:1 | z1:0 | z2:1 x1:1 | x2:2 | y1:0 | y2:1 | z1:0 | z2:1 x1:0 | x2:1 | y1:1 | y2:2 | z1:0 | z2:1 x1:1 | x2:2 | y1:1 | y2:2 | z1:0 | z2:1 x1:0 | x2:1 | y1:0 | y2:1 | z1:1 | z2:2 x1:1 | x2:2 | y1:0 | y2:1 | z1:1 | z2:2 x1:0 | x2:1 | y1:1 | y2:2 | z1:1 | z2:2 x1:1 | x2:2 | y1:1 | y2:2 | z1:1 | z2:2 """ nz, ny, nx = self.shape for i in xrange(nz): yield self.get_layer(i)
[docs] def dump(self, meshfile, propfile, prop): r""" Dump the mesh to a file in the format required by UBC-GIF program MeshTools3D. Parameters: * meshfile : str or file Output file to save the mesh. Can be a file name or an open file. * propfile : str or file Output file to save the physical properties *prop*. Can be a file name or an open file. * prop : str The name of the physical property in the mesh that will be saved to *propfile*. .. note:: Uses -10000000 as the dummy value for plotting topography Examples: >>> from StringIO import StringIO >>> meshfile = StringIO() >>> densfile = StringIO() >>> mesh = PrismMesh((0, 10, 0, 20, 0, 5), (1, 2, 2)) >>> mesh.addprop('density', [1, 2, 3, 4]) >>> mesh.dump(meshfile, densfile, 'density') >>> print meshfile.getvalue().strip() 2 2 1 0 0 0 2*10 2*5 1*5 >>> print densfile.getvalue().strip() 1.0000 3.0000 2.0000 4.0000 """ if prop not in self.props: raise ValueError("mesh doesn't have a '%s' property." % (prop)) isstr = False if isinstance(meshfile, str): isstr = True meshfile = open(meshfile, 'w') nz, ny, nx = self.shape x1, x2, y1, y2, z1, z2 = self.bounds dx, dy, dz = self.dims meshfile.writelines([ "%d %d %d\n" % (ny, nx, nz), "%g %g %g\n" % (y1, x1, -z1), "%d*%g\n" % (ny, dy), "%d*%g\n" % (nx, dx), "%d*%g" % (nz, dz)]) if isstr: meshfile.close() values = numpy.fromiter(self.props[prop], dtype=numpy.float) # Replace the masked cells with a dummy value values[self.mask] = -10000000 reordered = numpy.ravel(numpy.reshape(values, self.shape), order='F') numpy.savetxt(propfile, reordered, fmt='%.4f')
[docs] def copy(self): """ Return a deep copy of the current instance.""" return cp.deepcopy(self)
[docs]class TesseroidMesh(PrismMesh): """ Generate a 3D regular mesh of tesseroids. Tesseroids are ordered as follows: first layers (height coordinate), then N-S rows and finaly E-W. Ex: in a mesh with shape ``(3,3,3)`` the 15th element (index 14) has height index 1 (second layer), y index 1 (second row), and x index 2 ( third element in the column). This class can used as list of tesseroids. It acts as an iteratior (so you can loop over tesseroids). It also has a ``__getitem__`` method to access individual elements in the mesh. In practice, it should be able to be passed to any function that asks for a list of tesseroids, like :func:`fatiando.gravmag.tesseroid.gz`. To make the mesh incorporate a topography, use :meth:`~fatiando.mesher.TesseroidMesh.carvetopo` Parameters: * bounds : list = [w, e, s, n, top, bottom] Boundaries of the mesh. ``w, e, s, n`` in degrees, ``top`` and ``bottom`` are heights (positive upward) and in meters. * shape : tuple = (nr, nlat, nlon) Number of tesseroids in the radial, latitude, and longitude directions. * props : dict Physical properties of each tesseroid in the mesh. Each key should be the name of a physical property. The corresponding value should be a list with the values of that particular property on each tesseroid of the mesh. Examples: >>> from fatiando.mesher import TesseroidMesh >>> mesh = TesseroidMesh((0, 1, 0, 2, 3, 0), (1, 2, 2)) >>> for p in mesh: ... print p w:0 | e:0.5 | s:0 | n:1 | top:3 | bottom:0 w:0.5 | e:1 | s:0 | n:1 | top:3 | bottom:0 w:0 | e:0.5 | s:1 | n:2 | top:3 | bottom:0 w:0.5 | e:1 | s:1 | n:2 | top:3 | bottom:0 >>> print mesh[0] w:0 | e:0.5 | s:0 | n:1 | top:3 | bottom:0 >>> print mesh[-1] w:0.5 | e:1 | s:1 | n:2 | top:3 | bottom:0 One with physical properties:: >>> props = {'density':[2670.0, 1000.0]} >>> mesh = TesseroidMesh((0, 2, 0, 4, 3, 0), (1, 1, 2), props=props) >>> for p in mesh: ... print p w:0 | e:1 | s:0 | n:4 | top:3 | bottom:0 | density:2670 w:1 | e:2 | s:0 | n:4 | top:3 | bottom:0 | density:1000 or equivalently:: >>> mesh = TesseroidMesh((0, 2, 0, 4, 3, 0), (1, 1, 2)) >>> mesh.addprop('density', [200, -1000.0]) >>> for p in mesh: ... print p w:0 | e:1 | s:0 | n:4 | top:3 | bottom:0 | density:200 w:1 | e:2 | s:0 | n:4 | top:3 | bottom:0 | density:-1000 You can use :meth:`~fatiando.mesher.PrismMesh.get_xs` (and similar methods for y and z) to get the x coordinates of the tesseroidss in the mesh:: >>> mesh = TesseroidMesh((0, 2, 0, 4, 3, 0), (1, 1, 2)) >>> print mesh.get_xs() [ 0. 1. 2.] >>> print mesh.get_ys() [ 0. 4.] >>> print mesh.get_zs() [ 3. 0.] You can iterate over the layers of the mesh:: >>> mesh = TesseroidMesh((0, 2, 0, 2, 2, 0), (2, 2, 2)) >>> for layer in mesh.layers(): ... for p in layer: ... print p w:0 | e:1 | s:0 | n:1 | top:2 | bottom:1 w:1 | e:2 | s:0 | n:1 | top:2 | bottom:1 w:0 | e:1 | s:1 | n:2 | top:2 | bottom:1 w:1 | e:2 | s:1 | n:2 | top:2 | bottom:1 w:0 | e:1 | s:0 | n:1 | top:1 | bottom:0 w:1 | e:2 | s:0 | n:1 | top:1 | bottom:0 w:0 | e:1 | s:1 | n:2 | top:1 | bottom:0 w:1 | e:2 | s:1 | n:2 | top:1 | bottom:0 The ``shape`` of the mesh must be integer! >>> mesh = TesseroidMesh((0, 2, 0, 4, 0, 3), (1, 1, 2.5)) Traceback (most recent call last): ... AttributeError: Invalid mesh shape (1, 1, 2.5). shape must be integers """ celltype = Tesseroid def __init__(self, bounds, shape, props=None): PrismMesh.__init__(self, bounds, shape, props) self.zdown = False self.dump = None
[docs]def extract(prop, prisms): """ Extract the values of a physical property from the cells in a list. If a cell is ``None`` or doesn't have the physical property, a value of ``None`` will be put in it's place. Parameters: * prop : str The name of the physical property to extract * cells : list A list of cells (e.g., :class:`~fatiando.mesher.Prism`, :class:`~fatiando.mesher.PolygonalPrism`, etc) Returns: * values : array The extracted values Examples: >>> cells = [Prism(1, 2, 3, 4, 5, 6, {'foo':1}), ... Prism(1, 2, 3, 4, 5, 6, {'foo':10}), ... None, ... Prism(1, 2, 3, 4, 5, 6, {'bar':2000})] >>> print extract('foo', cells) [1, 10, None, None] """ def getprop(p): if p is None or prop not in p.props: return None return p.props[prop] return [getprop(p) for p in prisms]
[docs]def vfilter(vmin, vmax, prop, cells): """ Remove cells whose physical property value falls outside a given range. If a cell is `None` or doesn't have the physical property, it will be not be included in the result. Parameters: * vmin : float Minimum value * vmax : float Maximum value * prop : str The name of the physical property used to filter * cells : list A list of cells (e.g., :class:`~fatiando.mesher.Prism`, :class:`~fatiando.mesher.PolygonalPrism`, etc) Returns: * filtered : list The cells that fall within the desired range Examples: >>> cells = [Prism(1, 2, 3, 4, 5, 6, {'foo':1}), ... Prism(1, 2, 3, 4, 5, 6, {'foo':20}), ... Prism(1, 2, 3, 4, 5, 6, {'foo':3}), ... None, ... Prism(1, 2, 3, 4, 5, 6, {'foo':4}), ... Prism(1, 2, 3, 4, 5, 6, {'foo':200}), ... Prism(1, 2, 3, 4, 5, 6, {'bar':1000})] >>> for cell in vfilter(0, 10, 'foo', cells): ... print cell x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6 | foo:1 x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6 | foo:3 x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6 | foo:4 """ def isin(cell): if cell is None or prop not in cell.props: return False value = cell.props[prop] if not isinstance(value, float) and not isinstance(value, int): value = numpy.linalg.norm(cell.props[prop]) if value < vmin or value > vmax: return False return True return [c for c in cells if isin(c)]
[docs]def vremove(value, prop, cells): """ Remove the cells with a given physical property value. If a cell is `None` it will be not be included in the result. If a cell doesn't have the physical property, it will be included in the result. Parameters: * value : float The value of the physical property to remove. If the physical property is a vector, will compare the norm of the vector to **value**. * prop : str The name of the physical property to remove * cells : list A list of cells (e.g., :class:`~fatiando.mesher.Prism`, :class:`~fatiando.mesher.PolygonalPrism`, etc) Returns: * removed : list A list of cells that have *prop* != *value* Examples: >>> cells = [Prism(1, 2, 3, 4, 5, 6, {'foo':1}), ... Prism(1, 2, 3, 4, 5, 6, {'foo':20}), ... Prism(1, 2, 3, 4, 5, 6, {'foo':3}), ... None, ... Prism(1, 2, 3, 4, 5, 6, {'foo':1}), ... Prism(1, 2, 3, 4, 5, 6, {'foo':200}), ... Prism(1, 2, 3, 4, 5, 6, {'bar':1000})] >>> for cell in vremove(1, 'foo', cells): ... print cell x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6 | foo:20 x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6 | foo:3 x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6 | foo:200 x1:1 | x2:2 | y1:3 | y2:4 | z1:5 | z2:6 | bar:1000 """ def keep(cell): if cell is None: return False if prop not in cell.props: return True p = cell.props[prop] if not isinstance(p, float) and not isinstance(p, int): p = numpy.linalg.norm(cell.props[prop]) if p != value: return True return False return [c for c in cells if keep(c)]