Meshing (fatiando.mesher)

Generate and operate on various kinds of meshes and geometric elements

Geometric elements

Meshes

Utility functions

  • extract: Extract the values of a physical property from the cells in a list
  • vfilter: Remove cells whose physical property value falls outside a given range
  • vremove: Remove the cells with a given physical property value

class fatiando.mesher.GeometricElement(props)[source]

Bases: object

Base class for all geometric elements.

addprop(prop, value)[source]

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.

copy()[source]

Return a deep copy of the current instance.

class fatiando.mesher.PointGrid(area, z, shape, props=None)[source]

Bases: object

Create a regular grid of 3D point sources (spheres of unit volume).

Use this as a 1D list of 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)
addprop(prop, values)[source]

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

copy()[source]

Return a deep copy of the current instance.

split(shape)[source]

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:

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.])
class fatiando.mesher.Polygon(vertices, props=None)[source]

Bases: fatiando.mesher.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])
addprop(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.

copy()

Return a deep copy of the current instance.

class fatiando.mesher.PolygonalPrism(vertices, z1, z2, props=None)[source]

Bases: fatiando.mesher.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
addprop(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.

copy()

Return a deep copy of the current instance.

topolygon()[source]

Get the polygon describing the prism viewed from above.

Returns:

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.]
class fatiando.mesher.Prism(x1, x2, y1, y2, z1, z2, props=None)[source]

Bases: fatiando.mesher.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
addprop(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.

center()[source]

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. ]
copy()

Return a deep copy of the current instance.

get_bounds()[source]

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]
class fatiando.mesher.PrismMesh(bounds, shape, props=None)[source]

Bases: 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).

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, PrismMesh should be able to be passed to any function that asks for a list of prisms, like fatiando.gravmag.prism.gz.

To make the mesh incorporate a topography, use 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 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
addprop(prop, values)[source]

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 PrismMesh

carvetopo(x, y, height)[source]

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

celltype

alias of Prism

copy()[source]

Return a deep copy of the current instance.

dump(meshfile, propfile, prop)[source]

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
get_layer(i)[source]

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 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
get_xs()[source]

Return an array with the x coordinates of the prisms in mesh.

get_ys()[source]

Return an array with the y coordinates of the prisms in mesh.

get_zs()[source]

Return an array with the z coordinates of the prisms in mesh.

layers()[source]

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
class fatiando.mesher.PrismRelief(ref, dims, nodes)[source]

Bases: 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 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.

addprop(prop, values)[source]

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.

copy()[source]

Return a deep copy of the current instance.

class fatiando.mesher.Sphere(x, y, z, radius, props=None)[source]

Bases: fatiando.mesher.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
addprop(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.

copy()

Return a deep copy of the current instance.

class fatiando.mesher.Square(bounds, props=None)[source]

Bases: fatiando.mesher.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 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
addprop(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.

bounds

The x, y boundaries of the square as [xmin, xmax, ymin, ymax]

copy()

Return a deep copy of the current instance.

vertices

The vertices of the square.

class fatiando.mesher.SquareMesh(bounds, shape, props=None)[source]

Bases: object

Generate a 2D regular mesh of squares.

For all purposes, SquareMesh can be used as a list of 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
addprop(prop, values)[source]

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 SquareMesh

copy()[source]

Return a deep copy of the current instance.

get_xs()[source]

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.

get_ys()[source]

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.

img2prop(fname, vmin, vmax, prop)[source]

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

class fatiando.mesher.Tesseroid(w, e, s, n, top, bottom, props=None)[source]

Bases: fatiando.mesher.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
addprop(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.

copy()

Return a deep copy of the current instance.

get_bounds()[source]

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]
half(lon=True, lat=True, r=True)[source]

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
split(nlon, nlat, nh)[source]

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
class fatiando.mesher.TesseroidMesh(bounds, shape, props=None)[source]

Bases: fatiando.mesher.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 fatiando.gravmag.tesseroid.gz.

To make the mesh incorporate a topography, use 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 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
addprop(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 PrismMesh

carvetopo(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

celltype

alias of Tesseroid

copy()

Return a deep copy of the current instance.

dump(meshfile, propfile, prop)

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
get_layer(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 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
get_xs()

Return an array with the x coordinates of the prisms in mesh.

get_ys()

Return an array with the y coordinates of the prisms in mesh.

get_zs()

Return an array with the z coordinates of the prisms in mesh.

layers()

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
fatiando.mesher.extract(prop, prisms)[source]

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., Prism, 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]
fatiando.mesher.vfilter(vmin, vmax, prop, cells)[source]

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., Prism, 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
fatiando.mesher.vremove(value, prop, cells)[source]

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., Prism, 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