Add function to create a tesseroid layer

examples/forward/tesseroid_layer.py

@@ -0,0 +1,67 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Layer of tesseroids
+===================
+"""
+import boule as bl
+import ensaio
+import numpy as np
+import pygmt
+import verde as vd
+import xarray as xr
+
+import harmonica as hm
+
+fname = ensaio.fetch_earth_topography(version=1)
+topo = xr.load_dataarray(fname)
+
+region = (-78, -53, -57, -20)
+topo = topo.sel(latitude=slice(*region[2:]), longitude=slice(*region[:2]))
+
+ellipsoid = bl.WGS84
+
+longitude, latitude = np.meshgrid(topo.longitude, topo.latitude)
+reference = ellipsoid.geocentric_radius(latitude)
+surface = topo + reference
+density = xr.where(topo > 0, 2670.0, 2670.0 - 1040.0)
+
+tesseroids = hm.tesseroid_layer(
+    coordinates=(topo.longitude, topo.latitude),
+    surface=surface,
+    reference=reference,
+    properties={"density": density},
+)
+
+# Create a regular grid of computation points located at 10km above reference
+grid_longitude, grid_latitude = vd.grid_coordinates(region=region, spacing=0.5)
+grid_radius = ellipsoid.geocentric_radius(grid_latitude) + 10e3
+grid_coords = (grid_longitude, grid_latitude, grid_radius)
+
+# Compute gravity field of tesseroids on a regular grid of observation points
+gravity = tesseroids.tesseroid_layer.gravity(grid_coords, field="g_z")
+
+grid = pygmt.xyz2grd(
+    x=grid_longitude.flatten(),
+    y=grid_latitude.flatten(),
+    z=gravity.flatten(),
+    region=region,
+    spacing=(0.5, 0.5),
+)
+
+# Plot gravity field
+fig = pygmt.Figure()
+fig.grdimage(
+    grid,
+    projection="M15c",
+    cmap="viridis",
+    nan_transparent=True,
+)
+fig.basemap(frame=True)
+fig.colorbar(frame='af+l"Gravity [mGal]"', position="JCR")
+fig.coast(shorelines="0.5p,black", borders=["1/0.5p,black"])
+fig.show()

harmonica/__init__.py

@@ -15,6 +15,7 @@
 from .forward.prism import prism_gravity
 from .forward.prism_layer import DatasetAccessorPrismLayer, prism_layer
 from .forward.tesseroid import tesseroid_gravity
+from .forward.tesseroid_layer import DatasetAccessorTesseroidLayer, tesseroid_layer
 from .gravity_corrections import bouguer_correction
 from .io import load_icgem_gdf
 from .isostasy import isostasy_airy, isostatic_moho_airy