Add option to discard thin prisms of a layer when forward modelling

harmonica/forward/prism_layer.py

@@ -306,7 +306,13 @@ def update_top_bottom(self, surface, reference):
  self._obj.coords["bottom"] = (self.dims, bottom)
 
  def gravity(
- self, coordinates, field, progressbar=False, density_name="density", **kwargs
+ self,
+ coordinates,
+ field,
+ progressbar=False,
+ density_name="density",
+ thickness_threshold=None,
+ **kwargs,
  ):
  """
  Computes the gravity generated by the layer of prisms
@@ -317,6 +323,8 @@ def gravity(
  through the ``density_name`` argument.
  Ignores the prisms which ``top`` or ``bottom`` boundaries are
  ``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`.
 
  Parameters
@@ -334,6 +342,10 @@ def gravity(
  Name of the property layer (or ``data_var`` of the
  :class:`xarray.Dataset`) that will be used for the density of each
  prism in the layer. Default to ``"density"``
+ thickness_threshold : float or None
+ Prisms thinner than this threshold will be ignored in the
+ forward gravity calculation. If None, every prism with non-zero
+ volume will be considered. Default to None.
 
  Returns
  -------
@@ -354,6 +366,13 @@ def gravity(
  # Select only the boundaries and density elements for masked prisms
  boundaries = boundaries[mask.ravel()]
  density = density[mask]
+ # Discard thin prisms and their densities
+ if thickness_threshold is not None:
+ boundaries, density = _discard_thin_prisms(
+ boundaries,
+ density,
+ thickness_threshold,
+ )
  # Return gravity field of prisms
  return prism_gravity(
  coordinates,
@@ -469,3 +488,42 @@ def to_pyvista(self):
  for data_var in self._obj.data_vars
  }
  return prism_to_pyvista(prisms, properties=properties)
+
+
+def _discard_thin_prisms(
+ prisms,
+ density,
+ thickness_threshold,
+):
+ """
+ Discard prisms with a thickness below a threshold
+
+ Parameters
+ ----------
+ prisms : 2d-array
+ Array containing the boundaries of the prisms in the following order:
+ ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``.
+ The array must have the following shape: (``n_prisms``, 6), where
+ ``n_prisms`` is the total number of prisms.
+ density : 1d-array
+ Array containing the density of each prism in kg/m^3. Must have the
+ same size as the number of prisms.
+ thickness_threshold : float
+ Prisms thinner than this threshold will be discarded.
+
+ Returns
+ -------
+ prisms : 2d-array
+ A copy of the ``prisms`` array that doesn't include the thin prisms.
+ density : 1d-array
+ A copy of the ``density`` array that doesn't include the density values
+ for thin prisms.
+ """
+ west, east, south, north, bottom, top = tuple(prisms[:, i] for i in range(6))
+ # Mark prisms with thickness < threshold as null prisms
+ thickness = top - bottom
+ null_prisms = thickness < thickness_threshold
+ # Keep only thick prisms and their densities
+ prisms = prisms[np.logical_not(null_prisms), :]
+ density = density[np.logical_not(null_prisms)]
+ return prisms, density

harmonica/tests/test_prism_layer.py

@@ -17,6 +17,7 @@
 import xarray as xr
 
 from .. import prism_gravity, prism_layer
+from ..forward.prism_layer import _discard_thin_prisms
 
 try:
  import pyvista
@@ -465,3 +466,77 @@ def test_numba_progress_missing_error(dummy_layer):
  # Check if error is raised
  with pytest.raises(ImportError):
  layer.prism_layer.gravity(coordinates, field="g_z", progressbar=True)
+
+
+def test_gravity_discarded_thin_prisms(dummy_layer):
+ """
+ Check if gravity of prism layer after discarding thin prisms is correct.
+ """
+ coordinates = vd.grid_coordinates((1, 3, 7, 10), spacing=1, extra_coords=30.0)
+ prism_coords, surface, reference, density = dummy_layer
+
+ layer = prism_layer(
+ prism_coords, surface, reference, properties={"density": density}
+ )
+ # Check that result with no threshold is the same as with a threshold of 0
+ gravity_prisms_nothres = layer.prism_layer.gravity(coordinates, field="g_z")
+ gravity_prisms_0thres = layer.prism_layer.gravity(
+ coordinates, field="g_z", thickness_threshold=0
+ )
+ npt.assert_allclose(gravity_prisms_nothres, gravity_prisms_0thres)
+
+ # Check that gravity from manually removed prisms is the same as using a
+ # threshold
+ manually_removed_prisms = []
+ for _, j in enumerate(layer.prism_layer._to_prisms()):
+ if abs(j[5] - j[4]) >= 5.0:
+ manually_removed_prisms.append(j)
+ gravity_manually_removed = prism_gravity(
+ coordinates,
+ prisms=manually_removed_prisms,
+ density=[2670] * len(manually_removed_prisms),
+ field="g_z",
+ )
+ gravity_threshold_removed = layer.prism_layer.gravity(
+ coordinates, field="g_z", thickness_threshold=5
+ )
+ npt.assert_allclose(gravity_manually_removed, gravity_threshold_removed)
+
+
+def test_discard_thin_prisms():
+ """
+ Check if thin prisms are properly discarded.
+ """
+ # create set of 4 prisms (west, east, south, north, bottom, top)
+ prism_boundaries = np.array(
+ [
+ [-5000.0, 5000.0, -5000.0, 5000.0, 0.0, 55.1],
+ [5000.0, 15000.0, -5000.0, 5000.0, 0.0, 55.01],
+ [-5000.0, 5000.0, 5000.0, 15000.0, 0.0, 35.0],
+ [5000.0, 15000.0, 5000.0, 15000.0, 0.0, 84.0],
+ ]
+ )
+
+ # assign densities to each prism
+ densities = np.array([2306, 2122, 2190, 2069])
+
+ # drop prisms and respective densities thinner than 55.05
+ # (2nd and 3rd prisms)
+ thick_prisms, thick_densities = _discard_thin_prisms(
+ prism_boundaries,
+ densities,
+ thickness_threshold=55.05,
+ )
+
+ # manually remove prisms and densities of thin prisms
+ expected_prisms = np.array(
+ [
+ [-5000.0, 5000.0, -5000.0, 5000.0, 0.0, 55.1],
+ [5000.0, 15000.0, 5000.0, 15000.0, 0.0, 84.0],
+ ]
+ )
+ expected_densities = np.array([2306, 2069])
+
+ # check the correct prisms and densities were discarded
+ npt.assert_allclose(expected_prisms, thick_prisms)
+ npt.assert_allclose(expected_densities, thick_densities)