harmonica.tesseroid_layer

harmonica.tesseroid_layer(coordinates, surface, reference, properties=None)

Create a layer of tesseroids of equal size.

Build a regular grid of tesseroids (spherical prisms) with variable top and bottom boundaries and properties like density, magnetization, etc. The function returns a xarray.Dataset containing longitude, latitude, top and bottom coordinates, and all physical properties as data_var s. The longitude and latitude coordinates correspond to the location of the center of each tesseroid.

Unlike the harmonica.prism_layer, the top and bottom boundaries of tesseroids, provided via surface and reference, are not heights relative to a reference level (e.g., sea level). Instead, they are radii from the center of the body (ellipsoid or sphere) to the top and bottom boundaries of the tesseroids. To convert top and bottom boundaries from heights to radii, you can use the boule.geocentric_radius method from boule to calculate the radius based on the latitude, and add these values to your height values.

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

Parameters:

coordinatestuple

List containing the coordinates of the centers of the tesseroids in spherical coordinates in the following order longitude and latitude. All coordinates should be in degrees.

surface2d-array

Array used to create the uppermost boundary of the tesseroids layer. Values should be radii (not heights) in meters. On every point where surface is below reference, the surface value will be used to set the bottom boundary of that tesseroid, while the reference value will be used to set the top boundary of the tesseroid.

reference2d-array or float

Reference surface used to create the lowermost boundary of the tesseroids layer. It can be either a plane or an irregular surface passed as 2d array. Values must be radii (not heights) in meters.

propertiesdict or None

Dictionary containing the physical properties of the tesseroids. 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 tesseroid, the height of its top and bottom boundaries ans its corresponding physical properties.

See also

harmonica.DatasetAccessorTesseroidLayer

Examples

>>> import numpy as np
>>> import boule as bl
>>> # Create a synthetic relief
>>> longitude = np.linspace(0, 10, 5)
>>> latitude = np.linspace(2, 8, 4)
>>> mesh_longitude, mesh_latitude = np.meshgrid(longitude, latitude)
>>> surface_heights = np.arange(20, dtype=float).reshape((4, 5))
>>> # Convert heights to radii
>>> reference_radii = bl.WGS84.geocentric_radius(mesh_latitude)
>>> surface_radii = surface_heights + reference_radii
>>> # Define constant densities
>>> density = 2670.0 * np.ones_like(surface_radii)
>>> # Define a layer of tesseroids
>>> tesseroids = tesseroid_layer(
...     (longitude, latitude),
...     surface_radii,
...     reference=reference_radii,
...     properties={"density": density},
... )
>>> print(tesseroids) 
<xarray.Dataset>
Dimensions:   (latitude: 4, longitude: 5)
Coordinates:
  * latitude   (latitude) float64 2.0 4.0 6.0 8.0
  * longitude  (longitude) float64 0.0 2.5 5.0 7.5 10.0
    top       (latitude, longitude) float64 6.378e+06 6.378e+06 ... 6.378e+06
    bottom    (latitude, longitude) float64 6.378e+06 6.378e+06 ... 6.378e+06
Data variables:
    density   (latitude, longitude) float64 2.67e+03 2.67e+03 ... 2.67e+03
Attributes:
    longitude_units:  degrees
    latitude_units:   degrees
    radius_units:     meters
    properties_units:  SI
>>> # Get the boundaries of the layer (will exceed the region)
>>> boundaries = tesseroids.tesseroid_layer.boundaries
>>> list(float(b) for b in boundaries)
[-1.25, 11.25, 1.0, 9.0]
>>> # Get the boundaries of one of the tesseroids
>>> tesseroid = tesseroids.tesseroid_layer.get_tesseroid((0, 2))
>>> list(float(b) for b in tesseroid)
[3.75, 6.25, 1.0, 3.0, 6378111.2, 6378113.2]

Examples using harmonica.tesseroid_layer

Layer of tesseroids

Layer of tesseroids