Repository: fatiando/harmonica
Pull request #378 (2023-02-12 00:40:55)

Compute horizontal derivatives of grids using finite differences
Make functions `derivative_easting` and `derivative_northing` to compute the horizontal derivatives using central finite-differences through the `xarray.DataArray.differentiate` method by default. Add a `method` argument to let users choose if they want to use the FFT implementation or the finite differences one. This implementation mimics the one we had in the old `fatiando`.


santisoler: @LL-Geo, I'm planning to merge this PR during this week. If you have the chance to take a look at it, that would be nice. Otherwise, don't worry, we can always change this code afterwards.

Since horizontal derivatives using finite differences in the grid are more accurate than their FFT counterpart, I think it might be better to offer the finite differences option as the default one. This was also the default behavior in the old [`fatiando` horizontal derivatives functions](https://legacy.fatiando.org/api/gravmag.transform.html#fatiando.gravmag.transform.derivx).
santisoler: @LL-Geo, I'm merging this so I can start the release of Harmonica. Don't worry about the review. If you find something we should change, feel free to open an issue or a pr! Thanks again for all the work with the fft filters!

Patch:
diff --git a/harmonica/_transformations.py b/harmonica/_transformations.py
index 0ff87c8..4396a42 100644
--- a/harmonica/_transformations.py
+++ b/harmonica/_transformations.py
@@ -54,12 +54,19 @@ def derivative_upward(grid, order=1):
     return apply_filter(grid, derivative_upward_kernel, order=order)
 
 
-def derivative_easting(grid, order=1):
+def derivative_easting(grid, order=1, method="finite-diff"):
     """
-    Calculate the derivative of a potential field grid in the easting direction
+    Calculate the derivative of a regular grid in the easting direction
 
     Compute the spatial derivative in the easting direction of regular gridded
-    data using frequency domain calculations through Fast Fourier Transform.
+    data. It can compute using accurate central differences using
+    :func:`xarray.differentiate` or through frequency domain calculations
+    through Fast Fourier Transform.
+
+    .. important::
+
+        Choosing the finite differences option produces more accurate results
+        without border effects.
 
     Parameters
     ----------
@@ -70,13 +77,19 @@ def derivative_easting(grid, order=1):
         same units.
     order : int
         The order of the derivative. Default to 1.
+    method : str (optional)
+        Method that will be used for computing the easting derivative. It can
+        be either ``"finite-diff"``, for computing using
+        :func:`xarray.differentiate`, or ``"fft"``, for using FFT-based
+        filters.
+        Default ``"finite-diff"``.
 
     Returns
     -------
     derivative : :class:`xarray.DataArray`
         A :class:`xarray.DataArray` with the easting derivatives of the passed
         ``grid``. Its units are the same units of the ``grid`` per units of its
-        coordinates.
+        coordinates to the power of the passed ``order``.
 
     References
     ----------
@@ -86,16 +99,33 @@ def derivative_easting(grid, order=1):
     --------
     harmonica.filters.derivative_easting_kernel
     """
-    return apply_filter(grid, derivative_easting_kernel, order=order)
-
-
-def derivative_northing(grid, order=1):
+    if method == "finite-diff":
+        coordinate = _get_dataarray_coordinate(grid, dimension_index=1)
+        for _ in range(order):
+            grid = grid.differentiate(coord=coordinate)
+    elif method == "fft":
+        grid = apply_filter(grid, derivative_easting_kernel, order=order)
+    else:
+        raise ValueError(
+            f"Invalid method '{method}'. "
+            "Please select one from 'finite-diff' or 'fft'."
+        )
+    return grid
+
+
+def derivative_northing(grid, order=1, method="finite-diff"):
     """
-    Calculate the derivative of a potential field grid in the northing
-    direction
+    Calculate the derivative of a regular grid in the northing direction
 
     Compute the spatial derivative in the northing direction of regular gridded
-    data using frequency domain calculations through Fast Fourier Transform.
+    data. It can compute using accurate central differences using
+    :func:`xarray.differentiate` or through frequency domain calculations
+    through Fast Fourier Transform.
+
+    .. important::
+
+        Choosing the finite differences option produces more accurate results
+        without border effects.
 
     Parameters
     ----------
@@ -106,13 +136,19 @@ def derivative_northing(grid, order=1):
         same units.
     order : int
         The order of the derivative. Default to 1.
+    method : str (optional)
+        Method that will be used for computing the easting derivative. It can
+        be either ``"finite-diff"``, for computing using
+        :func:`xarray.differentiate`, or ``"fft"``, for using FFT-based
+        filters.
+        Default ``"finite-diff"``.
 
     Returns
     -------
     derivative : :class:`xarray.DataArray`
         A :class:`xarray.DataArray` with the northing derivatives of the passed
         ``grid``. Its units are the same units of the ``grid`` per units of its
-        coordinates.
+        coordinates to the power of the passed ``order``.
 
     References
     ----------
@@ -122,7 +158,18 @@ def derivative_northing(grid, order=1):
     --------
     harmonica.filters.derivative_northing_kernel
     """
-    return apply_filter(grid, derivative_northing_kernel, order=order)
+    if method == "finite-diff":
+        coordinate = _get_dataarray_coordinate(grid, dimension_index=0)
+        for _ in range(order):
+            grid = grid.differentiate(coord=coordinate)
+    elif method == "fft":
+        return apply_filter(grid, derivative_northing_kernel, order=order)
+    else:
+        raise ValueError(
+            f"Invalid method '{method}'. "
+            "Please select one from 'finite-diff' or 'fft'."
+        )
+    return grid
 
 
 def upward_continuation(grid, height_displacement):
@@ -290,3 +337,32 @@ def reduction_to_pole(
         magnetization_inclination=magnetization_inclination,
         magnetization_declination=magnetization_declination,
     )
+
+
+def _get_dataarray_coordinate(grid, dimension_index):
+    """
+    Return the name of the easting or northing coordinate in the grid
+
+    Parameters
+    ----------
+    grid : :class:`xarray.DataArray`
+        Regular grid
+    dimension_index : int
+        Index of the dimension for the desired coordinate of the regular grid.
+        Since the dimensions of the grid should be in the order of
+        *northing*, *easting*, then 0 corresponds to *northing* and 1 to
+        *easting*.
+    """
+    dim_name = grid.dims[dimension_index]
+    coords = [c for c in grid.coords if grid[c].dims == (dim_name,)]
+    if len(coords) > 1:
+        if dimension_index == 0:
+            direction = "northing"
+        else:
+            direction = "easting"
+        coords = "', '".join(coords)
+        raise ValueError(
+            f"Grid contains more than one coordinate along the '{direction}' "
+            f"direction: '{coords}'."
+        )
+    return coords[0]


Test patch:
diff --git a/harmonica/tests/test_transformations.py b/harmonica/tests/test_transformations.py
index 0f33d90..2906e05 100644
--- a/harmonica/tests/test_transformations.py
+++ b/harmonica/tests/test_transformations.py
@@ -18,6 +18,7 @@ import xrft
 
 from .. import point_gravity
 from .._transformations import (
+    _get_dataarray_coordinate,
     derivative_easting,
     derivative_northing,
     derivative_upward,
@@ -197,6 +198,51 @@ def fixture_sample_g_ee(sample_grid_coords, sample_sources):
     return g_ee.g_ee
 
 
+@pytest.mark.parametrize("index, expected_dimension", ([1, "easting"], [0, "northing"]))
+def test_get_dataarray_coordinate(index, expected_dimension, sample_potential):
+    """
+    Test the _get_dataarray_coordinate private function
+    """
+    dimension = _get_dataarray_coordinate(sample_potential, index)
+    assert dimension == expected_dimension
+
+
+@pytest.mark.parametrize("index, dimension", ([1, "easting"], [0, "northing"]))
+def test_get_dataarray_coordinate_invalid_grid(index, dimension, sample_potential):
+    """
+    Test if _get_dataarray_coordinate raises error when grid have an additional
+    coordinate that share one of the horizontal dimensions
+    """
+    # Add another horizontal coordinate that shares the same dimension
+    extra_coord = np.ones_like(sample_potential[dimension])
+    grid = sample_potential.assign_coords({"extra_coord": (dimension, extra_coord)})
+    # Check if function raises an error
+    err_msg = "Grid contains more than one coordinate along the"
+    with pytest.raises(ValueError, match=err_msg):
+        _get_dataarray_coordinate(grid, index)
+
+
+@pytest.mark.parametrize(
+    "dimension, derivative_func",
+    (["easting", derivative_easting], ["northing", derivative_northing]),
+)
+def test_horizontal_derivative_with_invalid_grid(
+    dimension, derivative_func, sample_potential
+):
+    """
+    Test if the horizontal derivative functions raise an error when passing
+    a grid that has an additional coordinate that share the horizontal
+    dimension and the "finite-diff" method is selected.
+    """
+    # Add another horizontal coordinate that shares the same dimension
+    extra_coord = np.ones_like(sample_potential[dimension])
+    grid = sample_potential.assign_coords({"extra_coord": (dimension, extra_coord)})
+    # Check if function raises an error
+    err_msg = "Grid contains more than one coordinate along the"
+    with pytest.raises(ValueError, match=err_msg):
+        derivative_func(grid, method="finite-diff")
+
+
 def test_derivative_upward(sample_potential, sample_g_z):
     """
     Test derivative_upward function against the synthetic model
@@ -247,9 +293,50 @@ def test_derivative_upward_order2(sample_potential, sample_g_zz):
     assert rms / np.abs(g_zz).max() < 0.015
 
 
-def test_derivative_easting(sample_potential, sample_g_e):
+@pytest.mark.parametrize(
+    "derivative_func",
+    [derivative_easting, derivative_northing],
+    ids=["derivative_easting", "derivative_northing"],
+)
+def test_invalid_method_horizontal_derivatives(sample_potential, derivative_func):
+    """
+    Test if passing and invalid method to horizontal derivatives raise an error
+    """
+    method = "bla"
+    err_msg = f"Invalid method '{method}'."
+    with pytest.raises(ValueError, match=err_msg):
+        derivative_func(sample_potential, method=method)
+
+
+def test_derivative_easting_finite_diff(sample_potential, sample_g_e):
     """
-    Test derivative_easting function against the synthetic model
+    Test derivative_easting function against the synthetic model using finite
+    differences
+    """
+    # Calculate easting derivative
+    derivative = derivative_easting(sample_potential)
+    # Compare against g_e
+    g_e = sample_g_e * 1e-5  # convert to SI units
+    rms = root_mean_square_error(derivative, g_e)
+    assert rms / np.abs(g_e).max() < 0.01
+
+
+def test_derivative_easting_finite_diff_order_2(sample_potential, sample_g_ee):
+    """
+    Test higher order of derivative_easting function against the sample grid
+    using finite differences
+    """
+    # Calculate second easting derivative
+    second_deriv = derivative_easting(sample_potential, order=2)
+    # Compare against g_e
+    g_ee = sample_g_ee * 1e-9  # convert to SI units
+    rms = root_mean_square_error(second_deriv, g_ee)
+    assert rms / np.abs(g_ee).max() < 0.1
+
+
+def test_derivative_easting_fft(sample_potential, sample_g_e):
+    """
+    Test derivative_easting function against the synthetic model using FFTs
     """
     # Pad the potential field grid to improve accuracy
     pad_width = {
@@ -297,6 +384,32 @@ def test_derivative_easting_order2(sample_potential, sample_g_ee):
     assert rms / np.abs(g_ee).max() < 0.1
 
 
+def test_derivative_northing_finite_diff(sample_potential, sample_g_n):
+    """
+    Test derivative_northing function against the synthetic model using finite
+    differences
+    """
+    # Calculate northing derivative
+    derivative = derivative_northing(sample_potential)
+    # Compare against g_e
+    g_n = sample_g_n * 1e-5  # convert to SI units
+    rms = root_mean_square_error(derivative, g_n)
+    assert rms / np.abs(g_n).max() < 0.01
+
+
+def test_derivative_northing_finite_diff_order_2(sample_potential, sample_g_nn):
+    """
+    Test higher order of derivative_northing function against the sample grid
+    using finite differences
+    """
+    # Calculate second northing derivative
+    second_deriv = derivative_northing(sample_potential, order=2)
+    # Compare against g_e
+    g_nn = sample_g_nn * 1e-9  # convert to SI units
+    rms = root_mean_square_error(second_deriv, g_nn)
+    assert rms / np.abs(g_nn).max() < 0.1
+
+
 def test_derivative_northing(sample_potential, sample_g_n):
     """
     Test derivative_northing function against the synthetic model
@@ -347,9 +460,11 @@ def test_derivative_northing_order2(sample_potential, sample_g_nn):
     assert rms / np.abs(g_nn).max() < 0.1
 
 
-def test_laplace(sample_potential):
+def test_laplace_fft(sample_potential):
     """
-    Test second order of derivative fullfill laplace equation
+    Test if second order of derivatives fulfill Laplace equation
+
+    We will use FFT computations only.
     """
     # Pad the potential field grid to improve accuracy
     pad_width = {
@@ -362,12 +477,13 @@ def test_laplace(sample_potential):
         pad_width=pad_width,
     )
     # Calculate second northing derivative and unpad it
+    method = "fft"
+    second_deriv_ee = derivative_easting(potential_padded, order=2, method=method)
+    second_deriv_nn = derivative_northing(potential_padded, order=2, method=method)
     second_deriv_zz = derivative_upward(potential_padded, order=2)
-    second_deriv_ee = derivative_easting(potential_padded, order=2)
-    second_deriv_nn = derivative_northing(potential_padded, order=2)
-    second_deriv_zz = xrft.unpad(second_deriv_zz, pad_width)
     second_deriv_ee = xrft.unpad(second_deriv_ee, pad_width)
     second_deriv_nn = xrft.unpad(second_deriv_nn, pad_width)
+    second_deriv_zz = xrft.unpad(second_deriv_zz, pad_width)
     # Compare g_nn + g_ee against -g_zz (trim the borders to ignore boundary
     # effects)
     trim = 6