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,8 +63,23 @@
     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+lm")
+fig.coast(shorelines="0.5p,black", borders=["1/0.5p,black"])
+
+fig.shift_origin(xshift="16c")
+
+# Plot gravity field
 maxabs = vd.maxabs(gravity.g_z)
 pygmt.makecpt(cmap="polar", series=(-maxabs, maxabs))
 fig.grdimage(
@@ -63,6 +88,6 @@
     nan_transparent=True,
 )
 fig.basemap(frame=True)
-fig.colorbar(frame='af+l"Gravity [mGal]"', position="JCR")
+fig.colorbar(frame='af+l"Gravity [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 a sphere. 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 an 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
@@ -125,6 +130,15 @@ And finally plot the computed gravitational field
          cmap="viridis",
       )
 
+   # 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 boundaries",
+   )
+   fig.legend()
+
    fig.colorbar(cmap=True, frame=["a200f50", "x+lmGal"])
    fig.coast(shorelines="1p,black")
 
@@ -179,6 +193,20 @@ And plot the results:
          cmap="viridis",
       )
 
+   # 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.colorbar(cmap=True, frame=["a1000f500", "x+lmGal"])
    fig.coast(shorelines="1p,black")
 
@@ -265,6 +293,20 @@ Finally, lets plot it:
          cmap="viridis",
       )
 
+   # 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.colorbar(cmap=True, frame=["a200f100", "x+lmGal"])
    fig.coast(shorelines="1p,black")
 

harmonica/__init__.py

@@ -27,7 +27,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 = []