Add a function to ignore the tesseroid with zero density or volume (#339

)

Discard null tesseroids before passing them to the forward modelling
functions. A null tesseroid is the one that has either zero volume or zero
density. Since they don't contribute to the generated field, we can save
some computation time by discarding them.

harmonica/forward/_tesseroid_utils.py

@@ -491,3 +491,37 @@ def _longitude_continuity(tesseroids):
     east[tess_to_be_changed] = ((east[tess_to_be_changed] + 180) % 360) - 180
     west[tess_to_be_changed] = ((west[tess_to_be_changed] + 180) % 360) - 180
     return tesseroids
+
+
+def _discard_null_tesseroids(tesseroids, density):
+    """
+    Discard tesseroid with zero volume or zero density
+
+    Parameters
+    ----------
+    tesseroids : 2d-array
+        Array containing the boundaries of the tesserois in the following
+        order: ``west``, ``east``, ``south``, ``north``, ``bottom``, ``top``
+        defined in a geocentric spherical coordinate system.
+    density : 1d-array
+        Array containing the density of each tesseroid in kg/m^3. Must have the
+        same size as the number of tesseroids.
+
+    Returns
+    -------
+    tesseroids : 2d-array
+        A copy of the ``tesseroids`` array that doesn't include the null
+        tesseroids (tesseroids with zero density or zero volume).
+    density : 1d-array
+        A copy of the ``density`` array that doesn't include the density values
+        for the null tesseroids (tesseroid with zero density or zero volume).
+    """
+    west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6))
+    # Mark prisms with zero volume as null prisms
+    null_tesseroids = (west == east) | (south == north) | (bottom == top)
+    # Mark prisms with zero density as null prisms
+    null_tesseroids[density == 0] = True
+    # Keep only non null prisms
+    tesseroids = tesseroids[np.logical_not(null_tesseroids), :]
+    density = density[np.logical_not(null_tesseroids)]
+    return tesseroids, density

harmonica/forward/tesseroid.py

@@ -15,6 +15,7 @@
     _adaptive_discretization,
     _check_points_outside_tesseroids,
     _check_tesseroids,
+    _discard_null_tesseroids,
     gauss_legendre_quadrature,
     glq_nodes_weights,
 )
@@ -167,6 +168,8 @@ def tesseroid_gravity(
                 "Number of elements in density ({}) ".format(density.size)
                 + "mismatch the number of tesseroids ({})".format(tesseroids.shape[0])
             )
+        # Discard null tesseroids (zero density or zero volume)
+        tesseroids, density = _discard_null_tesseroids(tesseroids, density)
     # Get GLQ unscaled nodes, weights and number of nodes for each small
     # tesseroid
     glq_nodes, glq_weights = glq_nodes_weights(GLQ_DEGREES)

harmonica/tests/test_tesseroid.py

@@ -15,6 +15,7 @@
 
 from ..constants import GRAVITATIONAL_CONST
 from ..forward._tesseroid_utils import (
+    _discard_null_tesseroids,
     _distance_tesseroid_point,
     _longitude_continuity,
     _split_tesseroid,
@@ -259,6 +260,41 @@ def test_stack_overflow_on_adaptive_discretization():
         )
 
 
+def test_discard_null_tesseroids():
+    """
+    Test if discarding invalid tesseroid works as expected
+    """
+    # Define a set of sample tesseroids including invalid ones
+    ellipsoid = boule.WGS84
+    top = ellipsoid.mean_radius
+    bottom = top - 1e3
+    tesseroids = np.array(
+        [
+            [-10, -5, -10, -5, bottom, top],  # ok tesseroid
+            [-10, -5, -5, 0, bottom, top],  # ok tesseroid (will set zero density)
+            [-10, -10, 0, 5, bottom, top],  # no volume due to easting boundaries
+            [-10, -5, 5, 5, bottom, top],  # no volume due to noting boundaries
+            [-5, 0, -10, -5, top, top],  # no volume due to radial boundaries
+            [-5, 0, -5, 0, bottom, top],  # ok tesseroid
+            [-5, -5, 5, 5, top, top],  # no volume due to multiple boundaries
+            [-5, 0, 5, 10, bottom, top],  # ok tesseroid
+        ]
+    )
+    densities = np.array([2400, 0, 2500, 2600, 2700, 2800, 2900, 3000])
+    tesseroids, densities = _discard_null_tesseroids(tesseroids, densities)
+    npt.assert_allclose(
+        tesseroids,
+        np.array(
+            [
+                [-10, -5, -10, -5, bottom, top],
+                [-5, 0, -5, 0, bottom, top],
+                [-5, 0, 5, 10, bottom, top],
+            ]
+        ),
+    )
+    npt.assert_allclose(densities, np.array([2400, 2800, 3000]))
+
+
 # --------------------------------------
 # Test tesseroid distance and dimensions
 # --------------------------------------