Add percentile parameter to maxabs() function (#524)

Allow returning a percentile of the absolute value of the provided
values, instead of just the maximum absolute value. For example, using
the 95% percentile instead of the maximum will create a colorbar that is
less washed out and reduces the effects of outliers in the data.
If multiple arrays are provided, they need to be concatenated before
calculating the quantiles, so it can be significantly slower than just
the maximum absolute value determination.

---------

Co-authored-by: Santiago Soler <santisoler@fastmail.com>

Author: mdtanker

pyproject.toml

@@ -22,3 +22,6 @@ notice = '''
 #
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #'''
+
+[tool.pytest.ini_options]
+doctest_optionflags = "NUMBER"

verde/tests/test_utils.py

@@ -23,6 +23,7 @@
     grid_to_table,
     kdtree,
     make_xarray_grid,
+    maxabs,
     meshgrid_from_1d,
     meshgrid_to_1d,
     parse_engine,
@@ -334,3 +335,60 @@ def test_check_ndim_easting_northing():
     northing = np.linspace(-5, 5, 16).reshape(4, 4)
     with pytest.raises(ValueError):
         get_ndim_horizontal_coords(easting, northing)
+
+
+def test_maxabs_nans():
+    """
+    Test maxabs handles nans correctly
+    """
+    assert float(maxabs((0, 100, 1, 2, np.nan))) == 100
+    assert float(maxabs((np.nan, -3.2, -1, -2, 3.1))) == 3.2
+    assert np.isnan(maxabs((np.nan, -3, -1, 3), nan=False))
+
+
+def test_maxabs_percentile():
+    """
+    Test maxabs with percentile option
+    """
+    # test generic functionality
+    data = np.arange(1, 101)
+    assert float(maxabs(data, percentile=100)) == 100
+    assert pytest.approx(float(maxabs(data, percentile=90)), 0.1) == 90
+    assert pytest.approx(float(maxabs(data, percentile=50)), 0.1) == 50
+
+    # test with nans
+    data_with_nans = np.append(data, np.nan)
+    assert float(maxabs(data_with_nans, percentile=100)) == 100
+    assert pytest.approx(float(maxabs(data_with_nans, percentile=90)), 0.1) == 90
+    assert pytest.approx(float(maxabs(data_with_nans, percentile=50)), 0.1) == 50
+    assert (
+        pytest.approx(float(maxabs(data_with_nans, percentile=90, nan=True)), 0.1) == 90
+    )
+    assert np.isnan(float(maxabs(data_with_nans, percentile=90, nan=False)))
+
+    # test with varying array sizes
+    assert (
+        pytest.approx(
+            float(maxabs([0, 1, 2, 3, 4], [[-2, 2], [0, 5]], percentile=80)), 0.1
+        )
+        == 3.4
+    )
+
+    # test invalid percentile types
+    msg = "Invalid 'percentile' of type"
+    with pytest.raises(TypeError, match=msg):
+        maxabs(data, percentile="90")
+    msg = "Invalid 'percentile' of type"
+    with pytest.raises(TypeError, match=msg):
+        maxabs(data, percentile=[90])
+    msg = "Invalid 'percentile' of type"
+    with pytest.raises(TypeError, match=msg):
+        maxabs(data, percentile=None)
+
+    # test invalid percentile values
+    msg = "Invalid 'percentile' value of"
+    with pytest.raises(ValueError, match=msg):
+        maxabs(data, percentile=-10)
+    msg = "Invalid 'percentile' value of"
+    with pytest.raises(ValueError, match=msg):
+        maxabs(data, percentile=110)

verde/utils.py

@@ -176,7 +176,7 @@ def variance_to_weights(variance, tol=1e-15, dtype="float64"):
     return tuple(weights)
 
 
-def maxabs(*args, nan=True):
+def maxabs(*args, nan=True, percentile=100):
     """
     Calculate the maximum absolute value of the given array(s).
 
@@ -187,11 +187,19 @@ def maxabs(*args, nan=True):
     args
         One or more arrays. If more than one are given, a single maximum will
         be calculated across all arrays.
+    nan : bool, optional
+        If True, will use the ``nan`` version of numpy functions to ignore
+        NaNs.
+    percentile : float, optional
+        Instead of return the maximum absolute value, return a given
+        percentile of the absolute values. Must be between 0 and 100. A
+        value of 100 (default) will give the maximum absolute value, while
+        a value of 50 will give the median of the absolute values.
 
     Returns
     -------
     maxabs : float
-        The maximum absolute value across all arrays.
+        The maximum (or percentile) absolute value across all arrays.
 
     Examples
     --------
@@ -216,14 +224,41 @@ def maxabs(*args, nan=True):
     >>> float(result)
     nan
 
+    If a more robust statistic is desired, you can use ``percentile`` to get
+    the value at a given percentile instead of the maximum.
+
+    >>> result = maxabs((1, -10, 25, 2, 3), percentile=95)
+    >>> float(result)
+    21.99
+    >>> result = maxabs((1, -10, 25, 2, 3), percentile=100)
+    >>> float(result)
+    25.0
+
     """
     arrays = [np.atleast_1d(i) for i in args]
-    if nan:
-        npmin, npmax = np.nanmin, np.nanmax
+
+    if percentile == 100:
+        if nan:
+            npmin, npmax = np.nanmin, np.nanmax
+        else:
+            npmin, npmax = np.min, np.max
+        absolute = [npmax(np.abs([npmin(i), npmax(i)])) for i in arrays]
+        return npmax(absolute)
     else:
-        npmin, npmax = np.min, np.max
-    absolute = [npmax(np.abs([npmin(i), npmax(i)])) for i in arrays]
-    return npmax(absolute)
+        if not isinstance(percentile, (int, float)):
+            raise TypeError(
+                f"Invalid 'percentile' of type '{type(percentile).__name__}'. Percentile must be a float or an integer."
+            )
+        if percentile < 0 or percentile > 100:
+            raise ValueError(
+                f"Invalid 'percentile' value of '{percentile}'. It must be between 0 and 100."
+            )
+        if nan:
+            nppercentile = np.nanpercentile
+        else:
+            nppercentile = np.percentile
+        combined_array = np.concatenate([np.abs(a.ravel()) for a in arrays])
+        return nppercentile(combined_array, percentile)
 
 
 def make_xarray_grid(