Add new tilt_angle transformation function

doc/api/index.rst

@@ -33,6 +33,7 @@ Apply well known transformations regular gridded potential fields data.
     gaussian_highpass
     reduction_to_pole
     total_gradient_amplitude
+    tilt
 
 Frequency domain filters
 ------------------------

examples/transformations/tilt.py

@@ -0,0 +1,122 @@
+# 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)
+#
+"""
+Tilt 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 tilt of the grid
+tilt_grid = hm.tilt(magnetic_grid_padded)
+
+# Unpad the tilt grid
+tilt_grid = xrft.unpad(tilt_grid, pad_width)
+
+# # Show the tilt
+# print("\nTilt:\n", tilt_grid)
+
+# Define the inclination and declination of the region by the time of the data
+# acquisition (1990).
+inclination, declination = -52.98, 6.51
+
+# Apply a reduction to the pole over the magnetic anomaly grid. We will assume
+# that the sources share the same inclination and declination as the
+# geomagnetic field.
+rtp_grid_padded = hm.reduction_to_pole(
+    magnetic_grid_padded, inclination=inclination, declination=declination
+)
+
+# Unpad the reduced to the pole grid
+rtp_grid = xrft.unpad(rtp_grid_padded, pad_width)
+
+# Compute the tilt of the padded rtp grid
+tilt_rtp_grid = hm.tilt(rtp_grid_padded)
+
+# Unpad the tilt grid
+tilt_rtp_grid = xrft.unpad(tilt_rtp_grid, pad_width)
+
+# # Show the tilt from RTP
+# print("\nTilt from RTP:\n", tilt_rtp_grid)
+
+# Plot original magnetic anomaly, its RTP, and the tilt of both
+fig = pygmt.Figure()
+with fig.subplot(nrows=2, ncols=2, figsize=("28c", "30c"), sharey="l"):
+    scale = 0.5 * vd.maxabs(magnetic_grid, rtp_grid)
+    with fig.set_panel(panel=0):
+        # Make colormap of data
+        pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+        # Plot magnetic anomaly grid
+        fig.grdimage(
+            grid=magnetic_grid,
+            projection="X?",
+            cmap=True,
+            # frame=["a", "+tMagnetic anomaly intensity (TMI)"],
+        )
+    with fig.set_panel(panel=1):
+        # Make colormap of data
+        pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+        # Plot reduced to the pole magnetic anomaly grid
+        fig.grdimage(
+            grid=rtp_grid,
+            projection="X?",
+            cmap=True,
+            # frame=["a", "+tTMI reduced to the pole (RTP)"],
+        )
+        # Add colorbar
+        fig.colorbar(
+            frame='af+l"Magnetic anomaly intensity [nT]"',
+            position="JMR+o1/1c+e",
+        )
+    scale = 0.6 * vd.maxabs(tilt_grid, tilt_rtp_grid)
+    with fig.set_panel(panel=2):
+        # Make colormap for tilt (saturate it a little bit)
+        pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+        # Plot tilt
+        fig.grdimage(
+            grid=tilt_grid,
+            projection="X?",
+            cmap=True,
+            # frame=["a", "+tTilt from TMI"],
+        )
+        # Add colorbar
+    with fig.set_panel(panel=3):
+        # Make colormap for tilt rtp (saturate it a little bit)
+        scale = 0.6 * vd.maxabs(tilt_rtp_grid)
+        pygmt.makecpt(cmap="polar+h", series=[-scale, scale], background=True)
+        # Plot tilt
+        fig.grdimage(
+            grid=tilt_rtp_grid,
+            projection="X?",
+            cmap=True,
+            # frame=["a", "+tTilt from RTP"],
+        )
+        # Add colorbar
+        fig.colorbar(
+            frame='af+l"Tilt [rad]"',
+            position="JMR+o1/1c+e",
+        )
+fig.show()

harmonica/__init__.py

@@ -27,6 +27,7 @@
     gaussian_highpass,
     gaussian_lowpass,
     reduction_to_pole,
+    tilt,
     total_gradient_amplitude,
     upward_continuation,
 )

harmonica/_transformations.py

@@ -343,7 +343,7 @@ def reduction_to_pole(
 
 def total_gradient_amplitude(grid):
     r"""
-    Calculate the total gradient amplitude a magnetic field grid
+    Calculates the total gradient amplitude of a potential field grid
 
     Compute the total gradient amplitude of a regular gridded potential field
     `M`. The horizontal derivatives are calculated though finite-differences
@@ -393,6 +393,69 @@ def total_gradient_amplitude(grid):
     return np.sqrt(gradient[0] ** 2 + gradient[1] ** 2 + gradient[2] ** 2)
 
 
+def tilt(grid):
+    r"""
+    Calculates the tilt of a potential field grid as defined
+    by Miller and Singh (1994)
+
+    Compute the tilt of a regular gridded potential field
+    `M`. 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
+    -------
+    tilt_grid : :class:`xarray.DataArray`
+        A :class:`xarray.DataArray` after calculating the
+        tilt of the passed ``grid``.
+
+    Notes
+    -----
+    The tilt is calculated as:
+
+    .. math::
+
+        tilt(f) = tan^{-1}\left(
+            \frac{
+                \frac{\partial M}{\partial z}}{
+                \sqrt{\frac{\partial M}{\partial x}^2 +
+                      \frac{\partial M}{\partial y}^2}}
+            \right)
+
+    where :math:`M` is the regularly gridded potential field.
+
+    When used on magnetic total field anomaly data, works best if the data is
+    reduced to the pole.
+
+    It's useful to plot the zero contour line of the tilt to represent possible
+    outlines of the source bodies. Use :func:`matplotlib.pyplot.contour` or
+    :func:`matplotlib.pyplot.tricontour` for this.
+
+    References
+    ----------
+    [Blakely1995]_
+    """
+    # Run sanity checks on the grid
+    grid_sanity_checks(grid)
+    # Calculate the gradients of the grid
+    gradient = (
+        derivative_easting(grid, order=1),
+        derivative_northing(grid, order=1),
+        derivative_upward(grid, order=1),
+    )
+    # Calculate and return the tilt
+    horiz_deriv = np.sqrt(gradient[0] ** 2 + gradient[1] ** 2)
+    tilt = np.arctan2(gradient[2], horiz_deriv)
+    return tilt
+
+
 def _get_dataarray_coordinate(grid, dimension_index):
     """
     Return the name of the easting or northing coordinate in the grid

harmonica/tests/test_transformations.py

@@ -25,6 +25,7 @@
     gaussian_highpass,
     gaussian_lowpass,
     reduction_to_pole,
+    tilt,
     total_gradient_amplitude,
     upward_continuation,
 )
@@ -591,6 +592,75 @@ def test_invalid_grid_with_nans(self, sample_potential):
             total_gradient_amplitude(sample_potential)
 
 
+class TestTilt:
+    """
+    Test tilt function
+    """
+
+    def test_against_synthetic(
+        self, sample_potential, sample_g_n, sample_g_e, sample_g_z
+    ):
+        """
+        Test tilt 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 the tilt and unpad it
+        tilt_grid = tilt(potential_padded)
+        tilt_grid = xrft.unpad(tilt_grid, pad_width)
+        # Compare against g_tilt (trim the borders to ignore boundary effects)
+        trim = 6
+        tilt_grid = tilt_grid[trim:-trim, trim:-trim]
+        g_e = sample_g_e[trim:-trim, trim:-trim] * 1e-5  # convert to SI units
+        g_n = sample_g_n[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_horiz_deriv = np.sqrt(g_e**2 + g_n**2)
+        g_tilt = np.arctan2(-g_z, g_horiz_deriv)
+        rms = root_mean_square_error(tilt_grid, g_tilt)
+        assert rms / np.abs(tilt_grid).max() < 0.1
+
+    def test_invalid_grid_single_dimension(self):
+        """
+        Check if tilt raises error on grid with single
+        dimension
+        """
+        x = np.linspace(0, 10, 11)
+        y = x**2
+        grid = xr.DataArray(y, coords={"x": x}, dims=("x",))
+        with pytest.raises(ValueError, match="Invalid grid with 1 dimensions."):
+            tilt(grid)
+
+    def test_invalid_grid_three_dimensions(self):
+        """
+        Check if tilt raises error on grid with three
+        dimensions
+        """
+        x = np.linspace(0, 10, 11)
+        y = np.linspace(-4, 4, 9)
+        z = np.linspace(20, 30, 5)
+        xx, yy, zz = np.meshgrid(x, y, z)
+        data = xx + yy + zz
+        grid = xr.DataArray(data, coords={"x": x, "y": y, "z": z}, dims=("y", "x", "z"))
+        with pytest.raises(ValueError, match="Invalid grid with 3 dimensions."):
+            tilt(grid)
+
+    def test_invalid_grid_with_nans(self, sample_potential):
+        """
+        Check if tilt raises error if grid contains nans
+        """
+        sample_potential.values[0, 0] = np.nan
+        with pytest.raises(ValueError, match="Found nan"):
+            tilt(sample_potential)
+
+
 class Testfilter:
     """
     Test filter result against the output from oasis montaj