Repository: fatiando/harmonica (python)

Issue:
Ignore prisms with user-defined volume/thickness threshold
**Description of the desired feature:**
<!--
Add as much detail as you can in your description below. If possible, include an
example of how you would like to use this feature (even better if it's a code
example).
-->

As discussed in the meeting today ([notes](https://hackmd.io/@fatiando/development-calls-2022#2022-12-07)), it would be great to allow users to vary the threshold for prisms of low volume/thickness. Currently, as introduced in PR #334, prisms with no volume are ignored when `prism_gravity` is called to save some unnecessary computation. In my use case, there are many prisms that are thin, but not exactly equal to 0, and thus are contributing to the run time, but not affecting the results significantly. 

It would be helpful if users could change some threshold parameter to account for this. 

## Option 1: 
Currently, zero-volume prisms are defined by this [line](https://github.com/fatiando/harmonica/blob/9529b5ace1056128bce000d2f8b9e88dd2b9016c/harmonica/forward/prism.py#L242): 
```python
# Mark prisms with zero volume as null prisms
null_prisms = (west == east) | (south == north) | (bottom == top)
```
We could change this to something similar to: 
```python
null_prisms = (abs(west - east) <= ew_thres) | (abs(south - north) <= ns_thres) | (abs(bottom - top) <= thick_thres)
```
where `ew_thres`, `ns_thres`, and `thick_thres` all default to 0. This would allow users to vary if they wanted.

## Option 2: 
Suggested by @santisoler, we could also implement this within the `prism_layer.gravity()` function. I assume this could be implement the same as the Option 1 example above.

**Are you willing to help implement and maintain this feature?** 
<!--
Every feature we add is code that we will have to maintain and keep updated.
Please let us know if you're willing to help maintain this feature in the future.
-->
Yes I'm happy to attempt this!

Pull request:
Add option to discard thin prisms of a layer when forward modelling
Add a new `thickness_threshold` parameter to the `prism_layer.gravity()` method
that allows the user to ignore prisms below a certain thickness. This is
implemented within a new `_discard_thin_prisms()`, function which is similar to
the `_discard_null_prisms()` function. Add tests that check whether the
calculated gravity is the same from manually removing prisms and removing
prisms based on a threshold, and whether providing no threshold still produces
the correct gravity.

I've added a parameter `thickness_threshold` to the `prism_layer.gravity()` function which allows the user to ignore prisms below a certain thickness. This is implemented within a new function `_discard_thin_prisms()`, which is similar to the function `_discard_null_prisms()` [here](https://github.com/fatiando/harmonica/blob/9529b5ace1056128bce000d2f8b9e88dd2b9016c/harmonica/forward/prism.py#L214). 

I've also added a test that checks:
*  whether the calculated gravity is the same from manually removing prisms and removing prisms based on a threshold        

*  whether providing no threshold still produces the correct gravity. 

This shouldn't affect performance if no threshold is supplied since the discard is skipped if `thickness_threshold is None`.

I'm not sure what the performance hit is when a threshold is supplied.

**Relevant issues/PRs:**
<!--
Example: "Fixes #1234" / "See also #345" / "Relevant to #111"
Use keywords (e.g., Fixes, Closes) to create the links and automatically
close issues when this PR is merged. 
See https://github.com/blog/1506-closing-issues-via-pull-requests
-->
Fixes #371 

Patch:
diff --git a/harmonica/forward/prism_layer.py b/harmonica/forward/prism_layer.py
index f2f41cc..b9fb207 100644
--- a/harmonica/forward/prism_layer.py
+++ b/harmonica/forward/prism_layer.py
@@ -306,7 +306,13 @@ class DatasetAccessorPrismLayer:
         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 @@ class DatasetAccessorPrismLayer:
         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 @@ class DatasetAccessorPrismLayer:
             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 @@ class DatasetAccessorPrismLayer:
         # 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 @@ class DatasetAccessorPrismLayer:
                 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


Test patch:
diff --git a/harmonica/tests/test_prism_layer.py b/harmonica/tests/test_prism_layer.py
index 634e563..7de1b07 100644
--- a/harmonica/tests/test_prism_layer.py
+++ b/harmonica/tests/test_prism_layer.py
@@ -17,6 +17,7 @@ import verde as vd
 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)