harmonica.prism_layer

harmonica.prism_layer#

harmonica.prism_layer(coordinates, surface, reference, properties=None)[source]#

Create a layer of prisms of equal size

Build a regular grid of prisms of equal size on the horizontal directions with variable top and bottom boundaries and properties like density, magnetization, etc. The function returns a xarray.Dataset containing easting, northing, top and bottom coordinates, and all physical properties as data_var s. The easting and northing coordinates correspond to the location of the center of each prism.

The prism_layer dataset accessor can be used to access special methods and attributes for the layer of prisms, like the horizontal dimensions of the prisms, getting the boundaries of each prisms, etc. See DatasetAccessorPrismLayer for the definition of these methods and attributes.

Parameters:
coordinatestuple

List containing the coordinates of the centers of the prisms in the following order: easting, northing. The arrays must be 1d arrays containing the coordinates of the centers per axis, or could be 2d arrays as the ones returned by numpy.meshgrid. All coordinates should be in meters and should define a regular grid.

surface2d-array

Array used to create the uppermost boundary of the prisms layer. All heights should be in meters. On every point where surface is below reference, the surface value will be used to set the bottom boundary of that prism, while the reference value will be used to set the top boundary of the prism.

referencefloat or 2d-array

Reference surface used to create the lowermost boundary of the prisms layer. It can be either a plane or an irregular surface passed as 2d array. Height(s) must be in meters.

propertiesdict or None

Dictionary containing the physical properties of the prisms. The keys must be strings that will be used to name the corresponding data_var inside the xarray.Dataset, while the values must be 2d-arrays. All physical properties must be passed in SI units. If None, no data_var will be added to the xarray.Dataset. Default is None.

Returns:
datasetxarray.Dataset

Dataset containing the coordinates of the center of each prism, the height of its top and bottom boundaries and its corresponding physical properties.

See also

harmonica.DatasetAccessorPrismLayer

Examples

>>> # Create a synthetic relief
>>> import numpy as np
>>> easting = np.linspace(0, 10, 5)
>>> northing = np.linspace(2, 8, 4)
>>> surface = np.arange(20, dtype=float).reshape((4, 5))
>>> density = 2670.0 * np.ones_like(surface)
>>> # Define a layer of prisms
>>> prisms = prism_layer(
...     (easting, northing),
...     surface,
...     reference=0,
...     properties={"density": density},
... )
>>> print(prisms) 
<xarray.Dataset>
Dimensions:   (northing: 4, easting: 5)
Coordinates:
  * easting   (easting) float64 0.0 2.5 5.0 7.5 10.0
  * northing  (northing) float64 2.0 4.0 6.0 8.0
    top       (northing, easting) float64 0.0 1.0 2.0 3.0 ... 17.0 18.0 19.0
    bottom    (northing, easting) float64 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
Data variables:
    density   (northing, easting) float64 2.67e+03 2.67e+03 ... 2.67e+03
Attributes:
    coords_units:      meters
    properties_units:  SI
>>> # Get the boundaries of the layer (will exceed the region)
>>> boundaries = prisms.prism_layer.boundaries
>>> list(float(b) for b in boundaries)
[-1.25, 11.25, 1.0, 9.0]
>>> # Get the boundaries of one of the prisms
>>> prism = prisms.prism_layer.get_prism((0, 2))
>>> list(float(b) for b in prism)
[3.75, 6.25, 1.0, 3.0, 0.0, 2.0]

Examples using harmonica.prism_layer#

Layer of prisms

Layer of prisms

Gravitational effect of topography

Gravitational effect of topography

Plot a prism layer in 3D

Plot a prism layer in 3D
Contents