Contents

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 XarrayAcessorPrismLayer for the definition of these methods and attributes.

Parameters

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

  • surface (2d-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.

  • reference (float 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.

  • properties (dict 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

dataset (xarray.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
:ivar ~prism_layer.coords_units: meters
:ivar ~prism_layer.properties_units: SI

>>> # Get the boundaries of the layer (will exceed the region)
>>> print(prisms.prism_layer.boundaries)
(-1.25, 11.25, 1.0, 9.0)
>>> # Get the boundaries of one of the prisms
>>> print(prisms.prism_layer.get_prism((0, 2)))
(3.75, 6.25, 1.0, 3.0, 0.0, 2.0)

Examples using harmonica.prism_layer

harmonica.tesseroid_gravity harmonica.DatasetAccessorPrismLayer