Add variable density tesseroids

doc/conf.py

@@ -42,6 +42,7 @@
 intersphinx_mapping = {
  "python": ("https://docs.python.org/3/", None),
  "numpy": ("https://docs.scipy.org/doc/numpy/", None),
+ "numba": ("https://numba.pydata.org/numba-doc/latest/", None),
  "scipy": ("https://docs.scipy.org/doc/scipy/reference", None),
  "pandas": ("http://pandas.pydata.org/pandas-docs/stable/", None),
  "xarray": ("http://xarray.pydata.org/en/stable/", None),

doc/install.rst

@@ -22,6 +22,7 @@ Dependencies
 * `numpy <http://www.numpy.org/>`__
 * `pandas <http://pandas.pydata.org/>`__
 * `numba <https://numba.pydata.org/>`__
+* `scipy <https://www.scipy.org/>`__
 * `xarray <https://xarray.pydata.org/>`__
 * `scikit-learn <https://scikit-learn.org>`__
 * `pooch <http://www.fatiando.org/pooch/>`__

doc/references.rst

@@ -12,6 +12,7 @@ References
 .. [Hofmann-WellenhofMoritz2006] Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd, corr. ed. 2006 edition ed.). Wien ; New York: Springer.
 .. [Nagy2000] Nagy, D., Papp, G. & Benedek, J.(2000). The gravitational potential and its derivatives for the prism. Journal of Geodesy 74: 552. doi:`10.1007/s001900000116 <https://doi.org/10.1007/s001900000116>`__
 .. [Nagy2002] Nagy, D., Papp, G. & Benedek, J.(2002). Corrections to "The gravitational potential and its derivatives for the prism". Journal of Geodesy. doi:`10.1007/s00190-002-0264-7 <https://doi.org/10.1007/s00190-002-0264-7>`__
+.. [Soler2019] Soler, S. R., Pesce, A., Gimenez, M. E., & Uieda, L. (2019). Gravitational field calculation in spherical coordinates using variable densities in depth, Geophysical Journal International. doi: `10.1093/gji/ggz277 <https://doi.org/10.1093/gji/ggz277>`__
 .. [Soler2021] Soler, S. R. and Uieda, L. (2021). Gradient-boosted equivalent sources, Geophysical Journal International. doi:`10.1093/gji/ggab297 <https://doi.org/10.1093/gji/ggab297>`__
 .. [Vajda2004] Vajda, P., Vaníček, P., Novák, P. and Meurers, B. (2004). On evaluation of Newton integrals in geodetic coordinates: Exact formulation and spherical approximation. Contributions to Geophysics and Geodesy, 34(4), 289-314.
 .. [TurcotteSchubert2014] Turcotte, D. L., & Schubert, G. (2014). Geodynamics (3 edition). Cambridge, United Kingdom: Cambridge University Press.

environment.yml

@@ -9,6 +9,7 @@ dependencies:
  - numpy
  - pandas
  - numba
+ - scipy
  - scikit-learn
  - pooch>=0.7.0
  - verde>=1.5.0

examples/forward/tesseroid_variable_density.py

@@ -0,0 +1,82 @@
+# 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)
+#
+"""
+Tesseroids with variable density
+================================
+
+The :func:`harmonica.tesseroid_gravity` is capable of computing the
+gravitational effects of tesseroids whose density is defined through
+a continuous function of the radial coordinate. This is achieved by the
+application of the method introduced in [Soler2021]_.
+
+To do so we need to define a regular Python function for the density, which
+should have a single argument (the ``radius`` coordinate) and return the
+density of the tesseroids at that radial coordinate.
+In addition, we need to decorate the density function with
+:func:`numba.jit(nopython=True)` or ``numba.njit`` for short.
+
+On this example we will show how we can compute the gravitational effect of
+four tesseroids whose densities are given by a custom linear ``density``
+function.
+"""
+import matplotlib.pyplot as plt
+import cartopy.crs as ccrs
+import verde as vd
+import boule as bl
+import harmonica as hm
+from numba import njit
+
+
+# Use the WGS84 ellipsoid to obtain the mean Earth radius which we'll use to
+# reference the tesseroid
+ellipsoid = bl.WGS84
+mean_radius = ellipsoid.mean_radius
+
+# Define tesseroid with top surface at the mean Earth radius, a thickness of
+# 10km and a linear density function
+tesseroids = (
+ [-70, -60, -40, -30, mean_radius - 3e3, mean_radius],
+ [-70, -60, -30, -20, mean_radius - 5e3, mean_radius],
+ [-60, -50, -40, -30, mean_radius - 7e3, mean_radius],
+ [-60, -50, -30, -20, mean_radius - 10e3, mean_radius],
+)
+
+# Define a linear density function. We should use the jit decorator so Numba
+# can run the forward model efficiently.
+
+
+@njit
+def density(radius):
+ """Linear density function"""
+ top = mean_radius
+ bottom = mean_radius - 10e3
+ density_top = 2670
+ density_bottom = 3000
+ slope = (density_top - density_bottom) / (top - bottom)
+ return slope * (radius - bottom) + density_bottom
+
+
+# Define computation points on a regular grid at 100km above the mean Earth
+# radius
+coordinates = vd.grid_coordinates(
+ region=[-80, -40, -50, -10],
+ shape=(80, 80),
+ extra_coords=100e3 + ellipsoid.mean_radius,
+)
+
+# Compute the radial component of the acceleration
+gravity = hm.tesseroid_gravity(coordinates, tesseroids, density, field="g_z")
+print(gravity)
+
+# Plot the gravitational field
+fig = plt.figure(figsize=(8, 9))
+ax = plt.axes(projection=ccrs.Orthographic(central_longitude=-60))
+img = ax.pcolormesh(*coordinates[:2], gravity, transform=ccrs.PlateCarree())
+plt.colorbar(img, ax=ax, pad=0, aspect=50, orientation="horizontal", label="mGal")
+ax.coastlines()
+ax.set_title("Downward component of gravitational acceleration")
+plt.show()