From 8e3f55c9296200221ced3e97132c58afcdbf41e6 Mon Sep 17 00:00:00 2001 From: ErosKerouak Date: Sun, 11 May 2025 10:41:42 -0300 Subject: [PATCH 1/6] Add total_horizontal_gradient function to _transformations.py Implements the computation of the total horizontal derivative (THDR) of a potential field grid. This function calculates the square root of the sum of squares of the horizontal derivatives in the easting and northing directions. --- harmonica/_transformations.py | 52 +++++++++++++++++++++++++++++++++++ 1 file changed, 52 insertions(+) diff --git a/harmonica/_transformations.py b/harmonica/_transformations.py index ea7c253b3..7e0495b3a 100644 --- a/harmonica/_transformations.py +++ b/harmonica/_transformations.py @@ -341,6 +341,58 @@ def reduction_to_pole( ) +def total_horizontal_gradient(grid): + r""" + Calculate the total horizontal derivative of a potential field grid. + + Compute the amplitude of the horizontal gradient of a regular gridded + potential field `M`. This is a measure of the rate of change in the x and y + (horizontal) directions. . The horizontal derivatives are calculated though + finite-differences. + + 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*. The coordinates must be defined in the same units. + + Returns + ------- + horizontal_derivative_grid : :class:`xarray.DataArray` + A :class:`xarray.DataArray` containing the total horizontal derivative of + the input ``grid``. + + Notes + ----- + The total horizontal derivative is calculated as: + + .. math:: + + A(x, y) = \sqrt{ + \left( \frac{\partial M}{\partial x} \right)^2 + + \left( \frac{\partial M}{\partial y} \right)^2 + } + + where :math:`M` is the regularly gridded potential field. + + References + ---------- + [Blakely1995]_ + [CordellGrauch1985]_ + """ + + # Run sanity checks on the grid + grid_sanity_checks(grid) + # Calculate the horizontal gradients of the grid + horizontal_gradient = ( + derivative_easting(grid, order=1), + derivative_northing(grid, order=1) + ) + # return the total horizontal gradient + return np.sqrt(horizontal_gradient[0] ** 2 + horizontal_gradient[1] ** 2) + + def total_gradient_amplitude(grid): r""" Calculates the total gradient amplitude of a potential field grid From 395610d75f990c0c26623e82a061366308dc64db Mon Sep 17 00:00:00 2001 From: ErosKerouak Date: Sun, 11 May 2025 11:10:11 -0300 Subject: [PATCH 2/6] Add tests for total_horizontal_gradient function Includes tests for: - Agreement with a synthetic model - Rejection of 1D and 3D grids - Handling of grids containing NaN values Ensures correctness and robustness of the total_horizontal_gradient filter. --- harmonica/tests/test_transformations.py | 65 +++++++++++++++++++++++++ 1 file changed, 65 insertions(+) diff --git a/harmonica/tests/test_transformations.py b/harmonica/tests/test_transformations.py index e1cc7c613..5a8c5c2c0 100644 --- a/harmonica/tests/test_transformations.py +++ b/harmonica/tests/test_transformations.py @@ -27,6 +27,7 @@ reduction_to_pole, tilt_angle, total_gradient_amplitude, + total_horizontal_gradient, upward_continuation, ) from .utils import root_mean_square_error @@ -602,6 +603,70 @@ def test_invalid_grid_with_nans(self, sample_potential): total_gradient_amplitude(sample_potential) +class TestTotalHorizontalGradient: + """ + Test total_horizontal_gradient function + """ + + def test_against_synthetic( + self, sample_potential, sample_g_n, sample_g_e + ): + """ + Test total_horizontal_gradient function against the synthetic model + """ + pad_width = { + "easting": sample_potential.easting.size // 3, + "northing": sample_potential.northing.size // 3, + } + potential_padded = xrft.pad( + sample_potential.drop_vars("upward"), + pad_width=pad_width, + ) + thg = total_horizontal_gradient(potential_padded) + thg = xrft.unpad(thg, pad_width) + + trim = 6 + thg = thg[trim:-trim, trim:-trim] + g_e = sample_g_e[trim:-trim, trim:-trim] * 1e-5 # convert to SI + g_n = sample_g_n[trim:-trim, trim:-trim] * 1e-5 + g_thg = np.sqrt(g_e**2 + g_n**2) + rms = root_mean_square_error(thg, g_thg) + assert rms / np.abs(g_thg).max() < 0.1 + + def test_invalid_grid_single_dimension(self): + """ + Check if total_horizontal_gradient 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."): + total_horizontal_gradient(grid) + + def test_invalid_grid_three_dimensions(self): + """ + Check if total_horizontal_gradient 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."): + total_horizontal_gradient(grid) + + def test_invalid_grid_with_nans(self, sample_potential): + """ + Check if total_horizontal_gradient raises error if grid contains nans + """ + sample_potential.values[0, 0] = np.nan + with pytest.raises(ValueError, match="Found nan"): + total_horizontal_gradient(sample_potential) + + class TestTilt: """ Test tilt function From 0bb4a5b8909cc938ce37fb3574b8a91d510107b5 Mon Sep 17 00:00:00 2001 From: ErosKerouak Date: Sun, 11 May 2025 17:39:25 -0300 Subject: [PATCH 3/6] Add theta_map function Implements the theta map filter. The function computes the arctangent of the horizontal gradient magnitude over the total gradient amplitude. Includes a detailed docstring with formula and references. --- harmonica/_transformations.py | 66 +++++++++++++++++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/harmonica/_transformations.py b/harmonica/_transformations.py index 7e0495b3a..8ec1ab29b 100644 --- a/harmonica/_transformations.py +++ b/harmonica/_transformations.py @@ -507,6 +507,72 @@ def tilt_angle(grid): return tilt +def theta_map(grid): + r""" + Calculate the theta map of a potential field grid + + Compute the theta map of a regularly gridded potential field + :math:`M`, used to enhance edges in gravity and magnetic data. + The horizontal and vertical derivatives are calculated from the + input grid, and the theta angle is defined as the arctangent + of the ratio between the total horizontal derivative and the + total gradient amplitude. + + 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 + ------- + theta_grid : :class:`xarray.DataArray` + A :class:`xarray.DataArray` with the calculated theta angle + in radians. + + Notes + ----- + The theta angle is calculated as: + + .. math:: + + \theta(M) = \tan^{-1} \left( + \frac{ + \sqrt{ + \left( \frac{\partial M}{\partial x} \right)^2 + + \left( \frac{\partial M}{\partial y} \right)^2 + } + }{ + \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 + } + } + \right) + + where :math:`M` is the regularly gridded potential field. + + References + ---------- + [Wijns et al., 2005]_ + """ + # 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 theta map + total_gradient = np.sqrt(gradient[0] ** 2 + gradient[1] ** 2 + gradient[2] ** 2) + horiz_deriv = np.sqrt(gradient[0] ** 2 + gradient[1] ** 2) + return np.arctan2(horiz_deriv, total_gradient) + + def _get_dataarray_coordinate(grid, dimension_index): """ Return the name of the easting or northing coordinate in the grid From 53dc3bde671a4091909b738ea098bffc277f5243 Mon Sep 17 00:00:00 2001 From: ErosKerouak Date: Sun, 11 May 2025 17:49:06 -0300 Subject: [PATCH 4/6] Add tests for theta_map function in test_transformations.py Includes: - Validation against synthetic model using finite differences and FFT - Input error handling for 1D, 3D, and NaN-containing grids - Integration of theta_map import for test discovery --- harmonica/tests/test_transformations.py | 65 +++++++++++++++++++++++++ 1 file changed, 65 insertions(+) diff --git a/harmonica/tests/test_transformations.py b/harmonica/tests/test_transformations.py index 5a8c5c2c0..3050932c0 100644 --- a/harmonica/tests/test_transformations.py +++ b/harmonica/tests/test_transformations.py @@ -26,6 +26,7 @@ gaussian_lowpass, reduction_to_pole, tilt_angle, + theta_map, total_gradient_amplitude, total_horizontal_gradient, upward_continuation, @@ -737,6 +738,70 @@ def test_invalid_grid_with_nans(self, sample_potential): with pytest.raises(ValueError, match="Found nan"): tilt_angle(sample_potential) +class TestThetaMap: + """ + Test theta_map function + """ + + def test_against_synthetic( + self, sample_potential, sample_g_n, sample_g_e, sample_g_z + ): + """ + Test theta_map function against the synthetic model + """ + pad_width = { + "easting": sample_potential.easting.size // 3, + "northing": sample_potential.northing.size // 3, + } + potential_padded = xrft.pad( + sample_potential.drop_vars("upward"), + pad_width=pad_width, + ) + theta_grid = theta_map(potential_padded) + theta_grid = xrft.unpad(theta_grid, pad_width) + + trim = 6 + theta_grid = theta_grid[trim:-trim, trim:-trim] + g_e = sample_g_e[trim:-trim, trim:-trim] * 1e-5 # SI units + g_n = sample_g_n[trim:-trim, trim:-trim] * 1e-5 + g_z = sample_g_z[trim:-trim, trim:-trim] * 1e-5 + g_thdr = np.sqrt(g_e**2 + g_n**2) + g_total = np.sqrt(g_thdr**2 + g_z**2) + g_theta = np.arctan2(g_thdr, g_total) + rms = root_mean_square_error(theta_grid, g_theta) + assert rms / np.abs(g_theta).max() < 0.1 + + def test_invalid_grid_single_dimension(self): + """ + Check if theta_map 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."): + theta_map(grid) + + def test_invalid_grid_three_dimensions(self): + """ + Check if theta_map 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."): + theta_map(grid) + + def test_invalid_grid_with_nans(self, sample_potential): + """ + Check if theta_map raises error if grid contains nans + """ + sample_potential.values[0, 0] = np.nan + with pytest.raises(ValueError, match="Found nan"): + theta_map(sample_potential) + class Testfilter: """ From 27fa37f2c26397e4ba0bc025e7176dc5c46ba1e9 Mon Sep 17 00:00:00 2001 From: ErosKerouak Date: Sun, 11 May 2025 18:19:34 -0300 Subject: [PATCH 5/6] Add horizontal_tilt_angle function to _transformations.py Implements the horizontal tilt angle (HTA) filter, which computes the arctangent of the horizontal gradient magnitude over the absolute value of the upward derivative. Includes a full docstring with LaTeX formula and reference. --- harmonica/_transformations.py | 59 +++++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) diff --git a/harmonica/_transformations.py b/harmonica/_transformations.py index 8ec1ab29b..61bd9ac28 100644 --- a/harmonica/_transformations.py +++ b/harmonica/_transformations.py @@ -572,6 +572,65 @@ def theta_map(grid): horiz_deriv = np.sqrt(gradient[0] ** 2 + gradient[1] ** 2) return np.arctan2(horiz_deriv, total_gradient) +def horizontal_tilt_angle(grid): + r""" + Calculate the horizontal tilt angle of a potential field grid + + Compute the horizontal tilt angle of a regularly gridded potential field + :math:`M`. This filter computing the arctangent of the ratio between the + horizontal gradient magnitude and the absolute value of the upward derivative. + + 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 + ------- + horiz_tilt_grid : :class:`xarray.DataArray` + A :class:`xarray.DataArray` with the calculated horizontal tilt angle + in radians. + + Notes + ----- + The horizontal tilt angle is calculated as: + + .. math:: + + \text{HTA}(M) = \tan^{-1} \left( + \frac{ + \sqrt{ + \left( \frac{\partial M}{\partial x} \right)^2 + + \left( \frac{\partial M}{\partial y} \right)^2 + } + }{ + \left| \frac{\partial M}{\partial z} \right| + } + \right) + + where :math:`M` is the regularly gridded potential field. + + References + ---------- + [Cooper & Cowan, 2006]_ + """ + # 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 horizontal tilt angle + horiz_deriv = np.sqrt(gradient[0] ** 2 + gradient[1] ** 2) + return np.arctan2(horiz_deriv, np.abs(gradient[2])) + + + def _get_dataarray_coordinate(grid, dimension_index): """ From fdf1df9f0e43f17ea6bb8e14bdae6cf51c865382 Mon Sep 17 00:00:00 2001 From: ErosKerouak Date: Sun, 11 May 2025 18:29:08 -0300 Subject: [PATCH 6/6] Add unit tests for horizontal_tilt_angle function Includes: - Validation against synthetic model - Input error handling for 1D, 3D, and NaN-filled grids --- harmonica/tests/test_transformations.py | 66 +++++++++++++++++++++++++ 1 file changed, 66 insertions(+) diff --git a/harmonica/tests/test_transformations.py b/harmonica/tests/test_transformations.py index 3050932c0..4a53112a5 100644 --- a/harmonica/tests/test_transformations.py +++ b/harmonica/tests/test_transformations.py @@ -27,6 +27,7 @@ reduction_to_pole, tilt_angle, theta_map, + horizontal_tilt_angle, total_gradient_amplitude, total_horizontal_gradient, upward_continuation, @@ -738,6 +739,7 @@ def test_invalid_grid_with_nans(self, sample_potential): with pytest.raises(ValueError, match="Found nan"): tilt_angle(sample_potential) + class TestThetaMap: """ Test theta_map function @@ -802,6 +804,70 @@ def test_invalid_grid_with_nans(self, sample_potential): with pytest.raises(ValueError, match="Found nan"): theta_map(sample_potential) +class TestHorizontalTiltAngle: + """ + Test horizontal_tilt_angle function + """ + + def test_against_synthetic( + self, sample_potential, sample_g_n, sample_g_e, sample_g_z + ): + """ + Test horizontal_tilt_angle function against the synthetic model + """ + pad_width = { + "easting": sample_potential.easting.size // 3, + "northing": sample_potential.northing.size // 3, + } + potential_padded = xrft.pad( + sample_potential.drop_vars("upward"), + pad_width=pad_width, + ) + hta_grid = horizontal_tilt_angle(potential_padded) + hta_grid = xrft.unpad(hta_grid, pad_width) + + trim = 6 + hta_grid = hta_grid[trim:-trim, trim:-trim] + g_e = sample_g_e[trim:-trim, trim:-trim] * 1e-5 + g_n = sample_g_n[trim:-trim, trim:-trim] * 1e-5 + g_z = sample_g_z[trim:-trim, trim:-trim] * 1e-5 + g_thdr = np.sqrt(g_e**2 + g_n**2) + g_hta = np.arctan2(g_thdr, np.abs(g_z)) + rms = root_mean_square_error(hta_grid, g_hta) + assert rms / np.abs(g_hta).max() < 0.1 + + def test_invalid_grid_single_dimension(self): + """ + Check if horizontal_tilt_angle 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."): + horizontal_tilt_angle(grid) + + def test_invalid_grid_three_dimensions(self): + """ + Check if horizontal_tilt_angle 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."): + horizontal_tilt_angle(grid) + + def test_invalid_grid_with_nans(self, sample_potential): + """ + Check if horizontal_tilt_angle raises error if grid contains nans + """ + sample_potential.values[0, 0] = np.nan + with pytest.raises(ValueError, match="Found nan"): + horizontal_tilt_angle(sample_potential) + + class Testfilter: """