From c001529d0134367172f7650a850adf4c5af5b21b Mon Sep 17 00:00:00 2001 From: Leonardo Miquelutti Date: Mon, 12 Aug 2024 17:09:18 -0300 Subject: [PATCH] Add new `tilt_angle` transformation function (#486) Add new `tilt_angle` transformation function that computes the tilt from a given a grid of potential field data. Add the function to the API index. Include tests and an example for the new function. --------- Co-authored-by: Leonardo Uieda Co-authored-by: Santiago Soler --- doc/api/index.rst | 1 + doc/references.rst | 1 + examples/transformations/tilt.py | 136 ++++++++++++++++++++++++ harmonica/__init__.py | 1 + harmonica/_transformations.py | 64 ++++++++++- harmonica/tests/test_transformations.py | 72 +++++++++++++ 6 files changed, 274 insertions(+), 1 deletion(-) create mode 100644 examples/transformations/tilt.py diff --git a/doc/api/index.rst b/doc/api/index.rst index c3a1c20fe..4764bbd58 100644 --- a/doc/api/index.rst +++ b/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 ------------------------ diff --git a/doc/references.rst b/doc/references.rst index e8f496db9..4ecfd447a 100644 --- a/doc/references.rst +++ b/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 `__ .. [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 `__ +.. [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 `__ .. [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 `__ .. [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 `__ .. [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 `__ diff --git a/examples/transformations/tilt.py b/examples/transformations/tilt.py new file mode 100644 index 000000000..93e6afd29 --- /dev/null +++ b/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() diff --git a/harmonica/__init__.py b/harmonica/__init__.py index d5782867a..367a96af1 100644 --- a/harmonica/__init__.py +++ b/harmonica/__init__.py @@ -28,6 +28,7 @@ gaussian_highpass, gaussian_lowpass, reduction_to_pole, + tilt_angle, total_gradient_amplitude, upward_continuation, ) diff --git a/harmonica/_transformations.py b/harmonica/_transformations.py index 3957e73ad..ea7c253b3 100644 --- a/harmonica/_transformations.py +++ b/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 diff --git a/harmonica/tests/test_transformations.py b/harmonica/tests/test_transformations.py index 05f04d1fe..e1cc7c613 100644 --- a/harmonica/tests/test_transformations.py +++ b/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