Add a finite-difference upward derivative

harmonica/_transformations.py

@@ -7,6 +7,8 @@
 """
 Apply transformations to regular grids of potential fields
 """
+import numpy as np
+
 from .filters._filters import (
     derivative_easting_kernel,
     derivative_northing_kernel,
@@ -19,7 +21,7 @@
 from .filters._utils import apply_filter
 
 
-def derivative_upward(grid, order=1):
+def derivative_upward(grid, order=1, method="finite-diff"):
     """
     Calculate the derivative of a potential field grid in the upward direction
 
@@ -51,7 +53,22 @@ def derivative_upward(grid, order=1):
     --------
     harmonica.filters.derivative_upward_kernel
     """
-    return apply_filter(grid, derivative_upward_kernel, order=order)
+    if method == "finite-diff":
+        spacing = np.mean([
+            np.abs(grid[dim][1] - grid[dim][0]) for dim in grid.dims
+        ])
+        for _ in range(order):
+            up = upward_continuation(grid, spacing / 2)
+            down = upward_continuation(grid, -spacing / 2)
+            grid = (up - down) / spacing
+    elif method == "fft":
+        grid = apply_filter(grid, derivative_upward_kernel, order=order)
+    else:
+        raise ValueError(
+            f"Invalid method '{method}'. "
+            "Please select one from 'finite-diff' or 'fft'."
+        )
+    return grid
 
 
 def derivative_easting(grid, order=1, method="finite-diff"):

harmonica/tests/test_transformations.py

@@ -258,7 +258,7 @@ def test_derivative_upward(sample_potential, sample_g_z):
         pad_width=pad_width,
     )
     # Calculate upward derivative and unpad it
-    derivative = derivative_upward(potential_padded)
+    derivative = derivative_upward(potential_padded, method="fft")
     derivative = xrft.unpad(derivative, pad_width)
     # Compare against g_z (trim the borders to ignore boundary effects)
     trim = 6
@@ -283,7 +283,7 @@ def test_derivative_upward_order2(sample_potential, sample_g_zz):
         pad_width=pad_width,
     )
     # Calculate second upward derivative and unpad it
-    second_deriv = derivative_upward(potential_padded, order=2)
+    second_deriv = derivative_upward(potential_padded, order=2, method="fft")
     second_deriv = xrft.unpad(second_deriv, pad_width)
     # Compare against g_zz (trim the borders to ignore boundary effects)
     trim = 6
@@ -293,6 +293,34 @@ def test_derivative_upward_order2(sample_potential, sample_g_zz):
     assert rms / np.abs(g_zz).max() < 0.015
 
 
+def test_derivative_upward_finite_diff(sample_potential, sample_g_z):
+    """
+    Test derivative_upward function with finite differences against the
+    synthetic model
+    """
+    derivative = derivative_upward(sample_potential, method="finite-diff")
+    # Compare against g_z (trim the borders to ignore boundary effects)
+    trim = 6
+    derivative = derivative[trim:-trim, trim:-trim]
+    g_z = sample_g_z[trim:-trim, trim:-trim] * 1e-5  # convert to SI units
+    rms = root_mean_square_error(derivative, g_z)
+    assert rms / np.abs(g_z).max() < 0.015
+
+
+def test_derivative_upward_finite_diff_order2(sample_potential, sample_g_zz):
+    """
+    Test higher order of derivative_upward function with finite differences
+    against the sample grid
+    """
+    second_deriv = derivative_upward(sample_potential, order=2, method="finite-diff")
+    # Compare against g_zz (trim the borders to ignore boundary effects)
+    trim = 6
+    second_deriv = second_deriv[trim:-trim, trim:-trim]
+    g_zz = sample_g_zz[trim:-trim, trim:-trim] * 1e-9  # convert to SI units
+    rms = root_mean_square_error(second_deriv, g_zz)
+    assert rms / np.abs(g_zz).max() < 0.015
+
+
 @pytest.mark.parametrize(
     "derivative_func",
     [derivative_easting, derivative_northing],