Add total gradient amplitude transformation

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

@@ -0,0 +1,71 @@
+# Copyright (c) 2018 The Harmonica Developers.
+# Distributed under the terms of the BSD 3-Clause License.
+# SPDX-License-Identifier: BSD-3-Clause
+#
+# This code is part of the Fatiando a Terra project (https://www.fatiando.org)
+#
+"""
+Total gradient amplitude of a regular grid
+===================================
+"""
+import ensaio
+import pygmt
+import verde as vd
+import xarray as xr
+import xrft
+
+import harmonica as hm
+
+# Fetch magnetic grid over the Lightning Creek Sill Complex, Australia using
+# Ensaio and load it with Xarray
+fname = ensaio.fetch_lightning_creek_magnetic(version=1)
+magnetic_grid = xr.load_dataarray(fname)
+
+# Pad the grid to increase accuracy of the FFT filter
+pad_width = {
+ "easting": magnetic_grid.easting.size // 3,
+ "northing": magnetic_grid.northing.size // 3,
+}
+# drop the extra height coordinate
+magnetic_grid_no_height = magnetic_grid.drop_vars("height")
+magnetic_grid_padded = xrft.pad(magnetic_grid_no_height, pad_width)
+
+# Compute the total gradient amplitude of the grid
+tga = hm.total_gradient_amplitude(magnetic_grid_padded)
+
+# Unpad the total gradient amplitude grid
+tga = xrft.unpad(tga, pad_width)
+
+# Show the total gradient amplitude
+print("\nTotal Gradient Amplitude:\n", tga)
+
+# 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):
+ # Make colormap of data
+ scale = 2500
+ pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+ # Plot magnetic anomaly grid
+ fig.grdimage(
+ grid=magnetic_grid,
+ projection="X?",
+ cmap=True,
+ )
+ # Add colorbar
+ fig.colorbar(
+ frame='af+l"Magnetic anomaly [nT]"',
+ position="JBC+h+o0/1c+e",
+ )
+ with fig.set_panel(panel=1):
+ # Make colormap for total gradient amplitude (saturate it a little bit)
+ scale = 0.6 * vd.maxabs(tga)
+ 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()

harmonica/__init__.py

@@ -27,6 +27,7 @@
  gaussian_highpass,
  gaussian_lowpass,
  reduction_to_pole,
+ total_gradient_amplitude,
  upward_continuation,
 )
 from ._utils import magnetic_angles_to_vec, magnetic_vec_to_angles

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,
@@ -339,6 +341,59 @@ def reduction_to_pole(
  )
 
 
+def total_gradient_amplitude(grid):
+ """
+ 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):
  """
  Return the name of the easting or northing coordinate in the grid

harmonica/filters/_filters.py

@@ -164,12 +164,17 @@ def derivative_northing_kernel(fft_grid, order=1):
  # Catch the dims of the Fourier transformed grid
  dims = fft_grid.dims
  # Grab the coordinates of the Fourier transformed grid
+ freq_easting = fft_grid.coords[dims[1]]
  freq_northing = fft_grid.coords[dims[0]]
  # Convert frequencies to wavenumbers
+ k_easting = 2 * np.pi * freq_easting
  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
+ # 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
+
+ return np.sqrt(da_filter_up**2 + da_filter_easting**2 + da_filter_northing**2)
 
 
 def upward_continuation_kernel(fft_grid, height_displacement):

harmonica/tests/test_transformations.py

@@ -25,6 +25,7 @@
  gaussian_highpass,
  gaussian_lowpass,
  reduction_to_pole,
+ total_gradient_amplitude,
  upward_continuation,
 )
 from .utils import root_mean_square_error
@@ -519,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