Meshing (`fatiando.mesher`)¶

Generate and operate on various kinds of meshes and geometric elements

Geometric elements

Polygon
Square
Prism
PolygonalPrism
Sphere
Tesseroid

Meshes

SquareMesh
PrismMesh
PrismRelief
TesseroidMesh
PointGrid

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:

subgrids
: list
List of 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.])

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:

polygon
: 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.]

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.

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

Meshing (fatiando.mesher)¶

Meshing (`fatiando.mesher`)¶