added total gradient amplitude function fatiando#470

- implemented the tga function in `transformations.py`
- created the tga.py example
- added the section Total gradient amplitude (aka Analytic Signal) to transformations.rst
- added one test to test_transformations.py

doc/user_guide/transformations.rst

@@ -1,4 +1,4 @@
-.. _transformations:
+	.. _transformations:
 
 Grid transformations
 ====================
@@ -396,6 +396,46 @@ Let's plot the results side by side:
     )
     plt.show()
 
+
+Total gradient amplitude (aka Analytic Signal)
+----------------------------------------------
+
+We can also calculate the total gradient amplitude of any magnetic anomaly grid.
+This transformation consists in obtaining the amplitude of the gradient of the
+magnetic field in all the three spatial directions by applying  
+
+.. math::
+    |A(x, y)| = \sqrt{\left( \frac{\partial M}{\partial x} \right)^2 + \left( \frac{\partial M}{\partial y} \right)^2 + \left( \frac{\partial M}{\partial z} \right)^2}
+
+
+We can apply it through the :func:`harmonica.total_gradient_amplitude` function.
+
+.. jupyter-execute::
+
+    tga_grid = hm.total_gradient_amplitude(
+        magnetic_grid_padded
+    )
+
+    # Unpad the total gradient amplitude grid
+    tga_grid = xrft.unpad(tga_grid, pad_width)
+    tga_grid
+
+And plot it:
+
+.. jupyter-execute::
+
+    import verde as vd
+
+    maxabs = vd.maxabs(tga_grid)
+    kwargs = dict(cmap="seismic", vmin=0, vmax=maxabs, add_colorbar=False)
+
+    tmp = tga_grid.plot(cmap="seismic", center=0, add_colorbar=False)
+    plt.gca().set_aspect("equal")
+    plt.title("Magnetic anomaly total gradient amplitude")
+    plt.gca().ticklabel_format(style="sci", scilimits=(0, 0))
+    plt.colorbar(tmp, label="nT")
+    plt.show()
+
 ----
 
 .. grid:: 2

examples/transformations/tga.py

@@ -5,7 +5,7 @@
 # This code is part of the Fatiando a Terra project (https://www.fatiando.org)
 #
 """
-Upward derivative of a regular grid
+Total gradient amplitude of a regular grid
 ===================================
 """
 import ensaio
@@ -30,16 +30,16 @@
 magnetic_grid_no_height = magnetic_grid.drop_vars("height")
 magnetic_grid_padded = xrft.pad(magnetic_grid_no_height, pad_width)
 
-# Compute the upward derivative of the grid
+# Compute the total gradient amplitude of the grid
 tga = hm.total_gradient_amplitude(magnetic_grid_padded)
 
-# Unpad the derivative grid
+# Unpad the total gradient amplitude grid
 tga = xrft.unpad(tga, pad_width)
 
-# Show the upward derivative
+# Show the total gradient amplitude
 print("\nTotal Gradient Amplitude:\n", tga)
 
-# Plot original magnetic anomaly and the upward derivative
+# Plot original magnetic anomaly and the total gradient amplitude
 fig = pygmt.Figure()
 with fig.subplot(nrows=1, ncols=2, figsize=("28c", "15c"), sharey="l"):
     with fig.set_panel(panel=0):
@@ -58,16 +58,14 @@
             position="JBC+h+o0/1c+e",
         )
     with fig.set_panel(panel=1):
-        # Make colormap for upward derivative (saturate it a little bit)
+        # Make colormap for total gradient amplitude (saturate it a little bit)
         scale = 0.6 * vd.maxabs(tga)
-        pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
-        # Plot upward derivative
+        pygmt.makecpt(cmap="polar+h", series=[0, scale], background=True)
+        # Plot total gradient amplitude
         fig.grdimage(grid=tga, projection="X?", cmap=True)
         # Add colorbar
         fig.colorbar(
             frame='af+l"Total Gradient Amplitude [nT/m]"',
             position="JBC+h+o0/1c+e",
         )
-fig.show()
-
-a = 2
+fig.show()

harmonica/_transformations.py

@@ -15,10 +15,10 @@
     gaussian_lowpass_kernel,
     reduction_to_pole_kernel,
     upward_continuation_kernel,
-    total_gradient_amplitude_kernel,
 )
 from .filters._utils import apply_filter
-# from harmonica._transformations import total_gradient_amplitude
+
+import numpy as np
 
 
 def derivative_upward(grid, order=1):
@@ -340,8 +340,59 @@ def reduction_to_pole(
         magnetization_declination=magnetization_declination,
     )
 
+
 def total_gradient_amplitude(grid):
-    return apply_filter(grid, total_gradient_amplitude_kernel)
+    """
+    Calculate the total gradient amplitude a magnetic field grid
+
+    Compute the total gradient amplitude of a regular gridded potential
+    field `M` through `A(x, y) = sqrt((dM/dx)^2 + (dM/dy)^2 + (dM/dz)^2)`.
+    So far the horizontal derivatives are calculated though finite-differences
+    while the upward derivative is calculated using fft.
+
+    Parameters
+    ----------
+    grid : :class:`xarray.DataArray`
+        A two dimensional :class:`xarray.DataArray` whose coordinates are
+        evenly spaced (regular grid). Its dimensions should be in the following
+        order: *northing*, *easting*. Its coordinates should be defined in the
+        same units.
+
+    Returns
+    -------
+    total_gradient_amplitude_grid : :class:`xarray.DataArray`
+        A :class:`xarray.DataArray` after calculating the
+        total gradient amplitude of the passed ``grid``.
+
+    References
+    ----------
+    [Blakely1995]_
+    """
+
+   # Catch the dims of the grid
+    dims = grid.dims
+    # Check if the array has two dimensions
+    if len(dims) != 2:
+        raise ValueError(
+            f"Invalid grid with {len(dims)} dimensions. "
+            + "The passed grid must be a 2 dimensional array."
+        )
+    # Check if the grid has nans
+    if np.isnan(grid).any():
+        raise ValueError(
+            "Found nan(s) on the passed grid. "
+            + "The grid must not have missing values before computing the "
+            + "Fast Fourier Transform."
+        )
+    # Calculate the gradients of the grid
+    gradient = (
+        derivative_easting(grid, order=1),
+        derivative_northing(grid, order=1),
+        derivative_upward(grid, order=1),
+    )
+    # return the total gradient amplitude
+    return np.sqrt(gradient[0]**2 + gradient[1]**2 + gradient[2]**2)
+
 
 def _get_dataarray_coordinate(grid, dimension_index):
     """
@@ -369,4 +420,4 @@ def _get_dataarray_coordinate(grid, dimension_index):
             f"Grid contains more than one coordinate along the '{direction}' "
             f"direction: '{coords}'."
         )
-    return coords[0]
+    return coords[0]

harmonica/filters/_filters.py

@@ -119,60 +119,6 @@ def derivative_easting_kernel(fft_grid, order=1):
 
 
 def derivative_northing_kernel(fft_grid, order=1):
-    r"""
-    Filter for northing derivative in frequency domain
-
-    Return a :class:`xarray.DataArray` with the values of the frequency domain
-    filter for computing the northing derivative. The filter is built upon the
-    frequency coordinates of the passed ``fft_grid`` and is defined as follows:
-
-    .. math::
-
-        g(\mathbf{k}) = (i k_n)^n
-
-    where :math:`\mathbf{k}` is the wavenumber vector
-    (:math:`\mathbf{k} = 2\pi \mathbf{f}` where :math:`\mathbf{f}` is the
-    frequency vector), :math:`k_n` is the northing wavenumber component of
-    :math:`\mathbf{k}`, :math:`i` is the imaginary unit and :math:`n` is the
-    order of the derivative.
-
-    Parameters
-    ----------
-    fft_grid : :class:`xarray.DataArray`
-        Array with the Fourier transform of the original grid.
-        Its dimensions should be in the following order:
-        *freq_northing*, *freq_easting*.
-        Use :func:`xrft.xrft.fft` and :func:`xrft.xrft.ifft` functions to
-        compute the Fourier Transform and its inverse, respectively.
-    order : int
-        The order of the derivative. Default to 1.
-
-    Returns
-    -------
-    da_filter : :class:`xarray.DataArray`
-        Array with the kernel for the northing derivative filter in frequency
-        domain.
-
-    References
-    ----------
-    [Blakely1995]_
-
-    See also
-    --------
-    harmonica.derivative_northing
-    """
-    # Catch the dims of the Fourier transformed grid
-    dims = fft_grid.dims
-    # Grab the coordinates of the Fourier transformed grid
-    freq_northing = fft_grid.coords[dims[0]]
-    # Convert frequencies to wavenumbers
-    k_northing = 2 * np.pi * freq_northing
-    # Compute the filter for northing derivative in frequency domain
-    da_filter = (k_northing * 1j) ** order
-    return da_filter
-
-
-def total_gradient_amplitude_kernel(fft_grid):
 
 
     # Catch the dims of the Fourier transformed grid
@@ -183,7 +129,7 @@ def total_gradient_amplitude_kernel(fft_grid):
     # Convert frequencies to wavenumbers
     k_easting = 2 * np.pi * freq_easting
     k_northing = 2 * np.pi * freq_northing    
-    # Compute the filter for upward derivative in frequency domain
+    # Compute the filter for derivatives in frequency domain
     da_filter_up = np.sqrt(k_easting**2 + k_northing**2)
     da_filter_easting = (k_easting * 1j)
     da_filter_northing = (k_northing * 1j)

harmonica/tests/test_transformations.py

@@ -26,7 +26,7 @@
     gaussian_lowpass,
     reduction_to_pole,
     upward_continuation,
-    # total_gradient_amplitude,
+    total_gradient_amplitude,
 )
 from .utils import root_mean_square_error
 
@@ -520,6 +520,34 @@ def test_upward_continuation(sample_g_z, sample_g_z_upward):
     xrt.assert_allclose(continuation, g_z_upward, atol=1e-8)
 
 
+def test_total_gradient_amplitude(sample_potential, sample_g_n, sample_g_e, sample_g_z):
+    """
+    Test total_gradient_amplitude function against the synthetic model
+    """
+    # Pad the potential field grid to improve accuracy
+    pad_width = {
+        "easting": sample_potential.easting.size // 3,
+        "northing": sample_potential.northing.size // 3,
+    }
+    # need to drop upward coordinate (bug in xrft)
+    potential_padded = xrft.pad(
+        sample_potential.drop_vars("upward"),
+        pad_width=pad_width,
+    )
+    # Calculate total gradient amplitude and unpad it
+    tga = total_gradient_amplitude(potential_padded)
+    tga = xrft.unpad(tga, pad_width)
+    # Compare against g_tga (trim the borders to ignore boundary effects)
+    trim = 6
+    tga = tga[trim:-trim, trim:-trim]
+    g_n = sample_g_n[trim:-trim, trim:-trim] * 1e-5  # convert to SI units
+    g_e = sample_g_e[trim:-trim, trim:-trim] * 1e-5  # convert to SI units
+    g_z = sample_g_z[trim:-trim, trim:-trim] * 1e-5  # convert to SI units
+    g_tga = np.sqrt(g_n**2 + g_e**2 + g_z**2)
+    rms = root_mean_square_error(tga, g_tga)
+    assert rms / np.abs(g_tga).max() < 0.1
+
+
 class Testfilter:
     """
     Test filter result against the output from oasis montaj