Update and refactor code and docs for tesseroids

doc/gallery_src/forward/tesseroid_layer.py

@@ -27,10 +27,20 @@
 ellipsoid = bl.WGS84
 
 longitude, latitude = np.meshgrid(topo.longitude, topo.latitude)
+
+# Compute radius of WGS84 ellipsoid at each lat/lon coordinates
 reference = ellipsoid.geocentric_radius(latitude)
+
+# Compute surface topography with respect to the center of the Earth
 surface = topo + reference
+
+# tesseroids above sea level have density of 2670 kg/m³ and tesseroids located below sea
+# level have density of -1630 kg/m³
 density = xr.where(topo > 0, 2670.0, 1040.0 - 2670.0)
 
+# Create a layer of tesseroids representing the topography
+# These have tops and bottoms defined by Earth's topography with respect to the center
+# of the Earth.
 tesseroids = hm.tesseroid_layer(
     coordinates=(topo.longitude, topo.latitude),
     surface=surface,
@@ -53,16 +63,32 @@
     extra_coords_names="radius",
 )
 
-# Plot gravity field
 fig = pygmt.Figure()
+
+# Plot tesseroid thickness
+fig.grdimage(
+    tesseroids.top - tesseroids.bottom,
+    projection="M15c",
+    nan_transparent=True,
+    cmap="viridis",
+    frame="+tTesseroid thickness",
+)
+fig.basemap(frame=True)
+fig.colorbar(frame="af+lthickness (m)")
+fig.coast(shorelines="0.5p,black", borders=["1/0.5p,black"])
+
+fig.shift_origin(xshift="17c")
+
+# Plot gravity field
 maxabs = vd.maxabs(gravity.g_z)
-pygmt.makecpt(cmap="polar", series=(-maxabs, maxabs))
+pygmt.makecpt(cmap="balance+h0", series=(-maxabs, maxabs))
 fig.grdimage(
     gravity.g_z,
     projection="M15c",
     nan_transparent=True,
+    frame="+tForward gravity",
 )
 fig.basemap(frame=True)
-fig.colorbar(frame='af+l"Gravity [mGal]"', position="JCR")
+fig.colorbar(frame="af+lGravity (mGal)")
 fig.coast(shorelines="0.5p,black", borders=["1/0.5p,black"])
 fig.show()

doc/user_guide/forward_modelling/tesseroid.rst

@@ -36,6 +36,11 @@ boundaries in the following order: *west*, *east*, *south*, *north*, *bottom*,
 *top*, where the former four are its longitudinal and latitudinal boundaries in
 decimal degrees and the latter two are the two radii given in meters.
 
+These two radii represent the top and bottom surfaces of the tesseroid, and should be
+given as distances from the center of the Earth. Note this is different from the
+vertical boundaries used for **prisms** in Cartesian coordinates, which are given as
+heights above or below some reference level (e.g., mean sea level or a reference ellipsoid).
+
 .. note::
 
    The :func:`harmonica.tesseroid_gravity` numerically computed the
@@ -45,7 +50,7 @@ decimal degrees and the latter two are the two radii given in meters.
 
 
 Lets define a single tesseroid and compute the gravitational potential
-it generates on a regular grid of computation points located at 10 km  above
+it generates on a regular grid of computation points located at 10 km above
 its *top* boundary.
 
 Get the WGS84 reference ellipsoid from :mod:`boule` so we can obtain its mean
@@ -127,6 +132,15 @@ And finally plot the computed gravitational field
 
    fig.colorbar(cmap=True, frame=["a200f50", "x+lmGal"])
    fig.coast(shorelines="1p,black")
+
+   # Plot edges of tesseroid
+   fig.plot(
+      x=[tesseroid[0], tesseroid[1], tesseroid[1], tesseroid[0], tesseroid[0]],
+      y=[tesseroid[2], tesseroid[2], tesseroid[3], tesseroid[3], tesseroid[2]],
+      pen="1p,red",
+      label="Tesseroid boundary",
+   )
+   fig.legend()
 
    fig.show()
 
@@ -182,6 +196,20 @@ And plot the results:
    fig.colorbar(cmap=True, frame=["a1000f500", "x+lmGal"])
    fig.coast(shorelines="1p,black")
 
+   # Plot edges of tesseroids
+   for i, tesseroid in enumerate(tesseroids):
+      if i == 0:
+         label="Tesseroid boundaries"
+      else:
+         label=None
+      fig.plot(
+         x=[tesseroid[0], tesseroid[1], tesseroid[1], tesseroid[0], tesseroid[0]],
+         y=[tesseroid[2], tesseroid[2], tesseroid[3], tesseroid[3], tesseroid[2]],
+         pen="1p,red",
+         label=label,
+      )
+   fig.legend()
+
    fig.show()
 
 
@@ -264,10 +292,23 @@ Finally, lets plot it:
          frame=["a", f"+t{title}"],
          cmap="viridis",
       )
-
    fig.colorbar(cmap=True, frame=["a200f100", "x+lmGal"])
    fig.coast(shorelines="1p,black")
 
+   # Plot edges of tesseroids
+   for i, tesseroid in enumerate(tesseroids):
+      if i == 0:
+         label="Tesseroid boundaries"
+      else:
+         label=None
+      fig.plot(
+         x=[tesseroid[0], tesseroid[1], tesseroid[1], tesseroid[0], tesseroid[0]],
+         y=[tesseroid[2], tesseroid[2], tesseroid[3], tesseroid[3], tesseroid[2]],
+         pen="1p,red",
+         label=label,
+      )
+   fig.legend()
+
    fig.show()
 
 ----

harmonica/__init__.py

@@ -21,7 +21,7 @@
 from ._forward.prism_gravity import prism_gravity
 from ._forward.prism_layer import DatasetAccessorPrismLayer, prism_layer
 from ._forward.prism_magnetic import prism_magnetic
-from ._forward.tesseroid import tesseroid_gravity
+from ._forward.tesseroid_gravity import tesseroid_gravity
 from ._forward.tesseroid_layer import DatasetAccessorTesseroidLayer, tesseroid_layer
 from ._gravity_corrections import bouguer_correction
 from ._io.icgem_gdf import load_icgem_gdf

harmonica/_forward/_tesseroid_utils.py

@@ -61,7 +61,7 @@ def gauss_legendre_quadrature(
     glq_nodes : list
         Unscaled location of GLQ nodes for each direction.
     glq_weights : list
-        GLQ weigths for each node for each direction.
+        GLQ weights for each node for each direction.
     kernel : func
         Kernel function for the gravitational field of point masses.
 

harmonica/_forward/_tesseroid_variable_density.py

@@ -61,7 +61,7 @@ def gauss_legendre_quadrature_variable_density(
     glq_nodes : list
         Unscaled location of GLQ nodes for each direction.
     glq_weights : list
-        GLQ weigths for each node for each direction.
+        GLQ weights for each node for each direction.
     kernel : func
         Kernel function for the gravitational field of point masses.
 
@@ -136,7 +136,7 @@ def density_based_discretization(tesseroids, density):
     Returns
     -------
     discretized_tesseroids : 2d-array
-        Array containing the coordinates of radially discretized tesseriods.
+        Array containing the coordinates of radially discretized tesseroids.
         Each row of the array will have the boundaries for each new tesseroid.
     """
     discretized_tesseroids = []

harmonica/_forward/prism_layer.py

@@ -326,7 +326,7 @@ def gravity(
         The density of the prisms will be assigned from the ``data_var`` chosen
         through the ``density_name`` argument.
         Ignores the prisms which ``top`` or ``bottom`` boundaries are
-        ``np.nan``s.
+        ``np.nan``'s.
         Prisms thinner than a given threshold can be optionally ignored through
         the ``thickness_threshold`` argument.
         All ``kwargs`` will be passed to :func:`harmonica.prism_gravity`.