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

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
 # --------------------------------------