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_angle
 
 Frequency domain filters
 ------------------------

doc/references.rst

@@ -11,6 +11,7 @@ References
 .. [Grombein2013] Grombein, T., Seitz, K., Heck, B. (2013), Optimized formulas for the gravitational field of a tesseroid, Journal of Geodesy. doi:`10.1007/s00190-013-0636-1 <https://doi.org/10.1007/s00190-013-0636-1>`__
 .. [Hofmann-WellenhofMoritz2006] Hofmann-Wellenhof, B., & Moritz, H. (2006). Physical Geodesy (2nd, corr. ed. 2006 edition ed.). Wien ; New York: Springer.
 .. [LiGotze2001] Li, X. and H. J. Gotze, 2001, Tutorial: Ellipsoid, geoid, gravity, geodesy, and geophysics, Geophysics, 66(6), p. 1660-1668, doi:`10.1190/1.1487109 <https://doi.org/10.1190/1.1487109>`__
+.. [MillerSingh1994] Miller, H. G., & Singh, V. (1994). Potential field tilt — A new concept for location of potential field sources. Journal of Applied Geophysics, 32(3), 213–217. doi:`10.1016/0926-9851(94)90002-7 <https://doi.org/10.1016/0926-9851(94)90002-7>`__
 .. [Nagy2000] Nagy, D., Papp, G. & Benedek, J.(2000). The gravitational potential and its derivatives for the prism. Journal of Geodesy 74: 552. doi:`10.1007/s001900000116 <https://doi.org/10.1007/s001900000116>`__
 .. [Nagy2002] Nagy, D., Papp, G. & Benedek, J.(2002). Corrections to "The gravitational potential and its derivatives for the prism". Journal of Geodesy. doi:`10.1007/s00190-002-0264-7 <https://doi.org/10.1007/s00190-002-0264-7>`__
 .. [Oliveira2021] Oliveira Jr, Vanderlei C. and Uieda, Leonardo and Barbosa, Valeria C. F.. Sketch of three coordinate systems: Geocentric Cartesian, Geocentric Geodetic, and Topocentric Cartesian. figshare. doi: `10.6084/m9.figshare.15044241.v1 <https://doi.org/10.6084/m9.figshare.15044241.v1>`__

examples/transformations/tilt.py

@@ -0,0 +1,136 @@
+# 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_angle(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_angle(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
+region = (
+ magnetic_grid.easting.values.min(),
+ magnetic_grid.easting.values.max(),
+ magnetic_grid.northing.values.min(),
+ magnetic_grid.northing.values.max(),
+)
+fig = pygmt.Figure()
+with fig.subplot(
+ nrows=2,
+ ncols=2,
+ subsize=("20c", "20c"),
+ sharex="b",
+ sharey="l",
+ margins=["1c", "1c"],
+):
+ 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", "+tTotal field anomaly grid"],
+ )
+ 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", "+tReduced to the pole (RTP)"],
+ )
+ # Add colorbar
+ label = "nT"
+ fig.colorbar(
+ frame=f"af+l{label}",
+ position="JMR+o1/-0.25c+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 of total field anomaly grid"],
+ )
+ with fig.set_panel(panel=3):
+ # Make colormap for tilt rtp (saturate it a little bit)
+ 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 of RTP grid"],
+ )
+ # Add colorbar
+ label = "rad"
+ fig.colorbar(
+ frame=f"af+l{label}",
+ position="JMR+o1/-0.25c+e",
+ )
+fig.show()

harmonica/__init__.py

@@ -28,6 +28,7 @@
  gaussian_highpass,
  gaussian_lowpass,
  reduction_to_pole,
+ tilt_angle,
  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,68 @@ def total_gradient_amplitude(grid):
  return np.sqrt(gradient[0] ** 2 + gradient[1] ** 2 + gradient[2] ** 2)
 
 
+def tilt_angle(grid):
+ r"""
+ Calculates the tilt angle of a potential field grid
+
+ Compute the tilt of a regular gridded potential field
+ :math:`M`. The horizontal derivatives are calculated
+ through 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` with the calculated tilt
+ in radians.
+
+ Notes
+ -----
+ The tilt is calculated as:
+
+ .. math::
+
+ \text{tilt}(f) = \tan^{-1} \left(
+ \frac{
+ \frac{\partial M}{\partial z}
+ }{
+ \sqrt{
+ \left( \frac{\partial M}{\partial x} \right)^2
+ +
+ \left( \frac{\partial M}{\partial y} \right)^2
+ }
+ }
+ \right)
+
+ where :math:`M` is the regularly gridded potential field.
+
+ References
+ ----------
+ [Blakely1995]_
+ [MillerSingh1994]_
+ """
+ # 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_angle,
  total_gradient_amplitude,
  upward_continuation,
 )
@@ -601,6 +602,77 @@ 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_angle(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]
+ g_n = sample_g_n[trim:-trim, trim:-trim]
+ g_z = sample_g_z[trim:-trim, trim:-trim]
+ g_horiz_deriv = np.sqrt(g_e**2 + g_n**2)
+ g_tilt = np.arctan2(
+ -g_z, g_horiz_deriv
+ ) # use -g_z to use the _upward_ derivative
+ 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_angle(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_angle(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_angle(sample_potential)
+
+
 class Testfilter:
  """
  Test filter result against the output from oasis montaj