From 38cd7465a5c7dc0a8cc9ec2303de50821a555c1e Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 14:22:00 -0700 Subject: [PATCH 01/13] Reorganize prism forward modelling into a submodule Group together prism-related pieces of code inside a new `harmonica._forward/prisms` module, organized in a similar way as the one for `ellipsoids`. --- src/harmonica/__init__.py | 9 ++-- src/harmonica/_forward/prisms/__init__.py | 13 ++++++ .../{prism_gravity.py => prisms/gravity.py} | 41 ++---------------- .../{prism_layer.py => prisms/layer.py} | 4 +- .../{prism_magnetic.py => prisms/magnetic.py} | 8 ++-- src/harmonica/_forward/prisms/utils.py | 43 +++++++++++++++++++ 6 files changed, 72 insertions(+), 46 deletions(-) create mode 100644 src/harmonica/_forward/prisms/__init__.py rename src/harmonica/_forward/{prism_gravity.py => prisms/gravity.py} (93%) rename src/harmonica/_forward/{prism_layer.py => prisms/layer.py} (99%) rename src/harmonica/_forward/{prism_magnetic.py => prisms/magnetic.py} (99%) create mode 100644 src/harmonica/_forward/prisms/utils.py diff --git a/src/harmonica/__init__.py b/src/harmonica/__init__.py index 6881b9cc0..191a5eca1 100644 --- a/src/harmonica/__init__.py +++ b/src/harmonica/__init__.py @@ -18,9 +18,12 @@ from ._forward.dipole import dipole_magnetic from ._forward.ellipsoids import Ellipsoid, ellipsoid_gravity, ellipsoid_magnetic from ._forward.point import point_gravity -from ._forward.prism_gravity import prism_gravity -from ._forward.prism_layer import DatasetAccessorPrismLayer, prism_layer -from ._forward.prism_magnetic import prism_magnetic +from ._forward.prisms import ( + DatasetAccessorPrismLayer, + prism_gravity, + prism_layer, + prism_magnetic, +) from ._forward.tesseroid_gravity import tesseroid_gravity from ._forward.tesseroid_layer import DatasetAccessorTesseroidLayer, tesseroid_layer from ._gravity_corrections import bouguer_correction diff --git a/src/harmonica/_forward/prisms/__init__.py b/src/harmonica/_forward/prisms/__init__.py new file mode 100644 index 000000000..23640b8dd --- /dev/null +++ b/src/harmonica/_forward/prisms/__init__.py @@ -0,0 +1,13 @@ +# 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) +# +""" +Forward modelling of rectangular prisms. +""" + +from .gravity import prism_gravity +from .magnetic import prism_magnetic +from .layer import DatasetAccessorPrismLayer, prism_layer diff --git a/src/harmonica/_forward/prism_gravity.py b/src/harmonica/_forward/prisms/gravity.py similarity index 93% rename from src/harmonica/_forward/prism_gravity.py rename to src/harmonica/_forward/prisms/gravity.py index fe98794e9..ee64393d8 100644 --- a/src/harmonica/_forward/prism_gravity.py +++ b/src/harmonica/_forward/prisms/gravity.py @@ -30,7 +30,8 @@ ) from numba import jit, prange -from .utils import initialize_progressbar +from ..utils import initialize_progressbar +from .utils import check_prisms # Define dictionary with available gravity fields for prisms FIELDS = { @@ -209,7 +210,7 @@ def prism_gravity( f"Number of elements in density ({density.size}) " + f"mismatch the number of prisms ({prisms.shape[0]})" ) - _check_prisms(prisms) + check_prisms(prisms) _check_singular_points(coordinates, prisms, field) # Discard null prisms (zero volume or zero density) prisms, density = _discard_null_prisms(prisms, density) @@ -448,40 +449,6 @@ def _any_singular_point_g_nz(coordinates, prisms): return False -def _check_prisms(prisms): - """ - Check if prisms boundaries are well defined. - - Parameters - ---------- - prisms : 2d-array - Array containing the boundaries of the prisms in the following order: - ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``. - The array must have the following shape: (``n_prisms``, 6), where - ``n_prisms`` is the total number of prisms. - """ - west, east, south, north, bottom, top = tuple(prisms[:, i] for i in range(6)) - err_msg = "Invalid prism or prisms. " - bad_we = west > east - bad_sn = south > north - bad_bt = bottom > top - if bad_we.any(): - err_msg += "The west boundary can't be greater than the east one.\n" - for prism in prisms[bad_we]: - err_msg += f"\tInvalid prism: {prism}\n" - raise ValueError(err_msg) - if bad_sn.any(): - err_msg += "The south boundary can't be greater than the north one.\n" - for prism in prisms[bad_sn]: - err_msg += f"\tInvalid prism: {prism}\n" - raise ValueError(err_msg) - if bad_bt.any(): - err_msg += "The bottom boundary can't be greater than the top one.\n" - for prism in prisms[bad_bt]: - err_msg += f"\tInvalid tesseroid: {prism}\n" - raise ValueError(err_msg) - - def _discard_null_prisms(prisms, density): """ Discard prisms with zero volume or zero density. @@ -494,7 +461,7 @@ def _discard_null_prisms(prisms, density): The array must have the following shape: (``n_prisms``, 6), where ``n_prisms`` is the total number of prisms. This array of prisms must have valid boundaries. - Run ``_check_prisms`` before. + Run ``check_prisms`` before. density : 1d-array Array containing the density of each prism in kg/m^3. Must have the same size as the number of prisms. diff --git a/src/harmonica/_forward/prism_layer.py b/src/harmonica/_forward/prisms/layer.py similarity index 99% rename from src/harmonica/_forward/prism_layer.py rename to src/harmonica/_forward/prisms/layer.py index 975960b50..f75e0ce37 100644 --- a/src/harmonica/_forward/prism_layer.py +++ b/src/harmonica/_forward/prisms/layer.py @@ -14,8 +14,8 @@ import verde as vd import xarray as xr -from ..visualization import prism_to_pyvista -from .prism_gravity import prism_gravity +from ...visualization import prism_to_pyvista +from .gravity import prism_gravity def prism_layer( diff --git a/src/harmonica/_forward/prism_magnetic.py b/src/harmonica/_forward/prisms/magnetic.py similarity index 99% rename from src/harmonica/_forward/prism_magnetic.py rename to src/harmonica/_forward/prisms/magnetic.py index 65b139389..43fc5e344 100644 --- a/src/harmonica/_forward/prism_magnetic.py +++ b/src/harmonica/_forward/prisms/magnetic.py @@ -14,8 +14,8 @@ from choclo.prism import magnetic_e, magnetic_field, magnetic_n, magnetic_u from numba import jit, prange -from .prism_gravity import _check_prisms -from .utils import initialize_progressbar +from ..utils import initialize_progressbar +from .utils import check_prisms VALID_FIELDS = ("b", "b_e", "b_n", "b_u") FORWARD_FUNCTIONS = { @@ -412,7 +412,7 @@ def _discard_null_prisms(prisms, magnetization): The array must have the following shape: (``n_prisms``, 6), where ``n_prisms`` is the total number of prisms. This array of prisms must have valid boundaries. - Run ``_check_prisms`` before. + Run ``check_prisms`` before. magnetization : 2d-array Array containing the magnetization vector of each prism in :math:`Am^{-1}`. Each vector will be a row in the 2d-array. @@ -454,7 +454,7 @@ def _run_sanity_checks(prisms, magnetization): f"Number of magnetization vectors ({magnetization[0].size}) " + f"mismatch the number of prisms ({prisms.shape[0]})" ) - _check_prisms(prisms) + check_prisms(prisms) # Define jitted versions of the forward modelling function diff --git a/src/harmonica/_forward/prisms/utils.py b/src/harmonica/_forward/prisms/utils.py new file mode 100644 index 000000000..9a2468295 --- /dev/null +++ b/src/harmonica/_forward/prisms/utils.py @@ -0,0 +1,43 @@ +# 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) +# +""" +Utilities for forward modelling of rectangular prisms. +""" + + +def check_prisms(prisms): + """ + Check if prisms boundaries are well defined. + + Parameters + ---------- + prisms : 2d-array + Array containing the boundaries of the prisms in the following order: + ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``. + The array must have the following shape: (``n_prisms``, 6), where + ``n_prisms`` is the total number of prisms. + """ + west, east, south, north, bottom, top = tuple(prisms[:, i] for i in range(6)) + err_msg = "Invalid prism or prisms. " + bad_we = west > east + bad_sn = south > north + bad_bt = bottom > top + if bad_we.any(): + err_msg += "The west boundary can't be greater than the east one.\n" + for prism in prisms[bad_we]: + err_msg += f"\tInvalid prism: {prism}\n" + raise ValueError(err_msg) + if bad_sn.any(): + err_msg += "The south boundary can't be greater than the north one.\n" + for prism in prisms[bad_sn]: + err_msg += f"\tInvalid prism: {prism}\n" + raise ValueError(err_msg) + if bad_bt.any(): + err_msg += "The bottom boundary can't be greater than the top one.\n" + for prism in prisms[bad_bt]: + err_msg += f"\tInvalid tesseroid: {prism}\n" + raise ValueError(err_msg) From ef33716554dbb011c8888fd1dc9005261b721d2e Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 14:22:59 -0700 Subject: [PATCH 02/13] Update tests for the new structure --- test/test_prism.py | 25 +++++++++++-------------- test/test_prism_layer.py | 2 +- test/test_prism_magnetic.py | 2 +- 3 files changed, 13 insertions(+), 16 deletions(-) diff --git a/test/test_prism.py b/test/test_prism.py index 077e1470f..daeb950ac 100644 --- a/test/test_prism.py +++ b/test/test_prism.py @@ -33,12 +33,9 @@ except ImportError: # pragma: no cover ProgressBar = None -from harmonica import bouguer_correction -from harmonica._forward.prism_gravity import ( - _check_prisms, - _discard_null_prisms, - prism_gravity, -) +from harmonica import bouguer_correction, prism_gravity +from harmonica._forward.prisms.gravity import _discard_null_prisms +from harmonica._forward.prisms.utils import check_prisms from .utils import run_only_with_numba @@ -75,24 +72,24 @@ def test_invalid_density_array(): def test_invalid_prisms(): - """Check if invalid prism boundaries are caught by _check_prisms.""" + """Check if invalid prism boundaries are caught by check_prisms.""" w, e, s, n, bottom, top = -100, 100, -100, 100, -200, -100 # Check if it works properly on valid prisms - _check_prisms(np.atleast_2d([w, e, s, n, bottom, top])) + check_prisms(np.atleast_2d([w, e, s, n, bottom, top])) # Check if it works properly on valid prisms with zero volume - _check_prisms(np.atleast_2d([w, w, s, n, bottom, top])) - _check_prisms(np.atleast_2d([w, e, s, s, bottom, top])) - _check_prisms(np.atleast_2d([w, e, s, n, bottom, bottom])) + check_prisms(np.atleast_2d([w, w, s, n, bottom, top])) + check_prisms(np.atleast_2d([w, e, s, s, bottom, top])) + check_prisms(np.atleast_2d([w, e, s, n, bottom, bottom])) # Test invalid boundaries msg = "The west boundary can't be greater than the east one" with pytest.raises(ValueError, match=msg): - _check_prisms(np.atleast_2d([e, w, s, n, bottom, top])) + check_prisms(np.atleast_2d([e, w, s, n, bottom, top])) msg = "The south boundary can't be greater than the north one" with pytest.raises(ValueError, match=msg): - _check_prisms(np.atleast_2d([w, e, n, s, bottom, top])) + check_prisms(np.atleast_2d([w, e, n, s, bottom, top])) msg = "The bottom boundary can't be greater than the top one" with pytest.raises(ValueError, match=msg): - _check_prisms(np.atleast_2d([w, e, s, n, top, bottom])) + check_prisms(np.atleast_2d([w, e, s, n, top, bottom])) def test_discard_null_prisms(): diff --git a/test/test_prism_layer.py b/test/test_prism_layer.py index acd560dd0..0192fe655 100644 --- a/test/test_prism_layer.py +++ b/test/test_prism_layer.py @@ -18,7 +18,7 @@ import xarray as xr from harmonica import prism_gravity, prism_layer -from harmonica._forward.prism_layer import _discard_thin_prisms +from harmonica._forward.prisms.layer import _discard_thin_prisms try: import pyvista diff --git a/test/test_prism_magnetic.py b/test/test_prism_magnetic.py index 89f847ead..505ccf9e0 100644 --- a/test/test_prism_magnetic.py +++ b/test/test_prism_magnetic.py @@ -21,7 +21,7 @@ ProgressBar = None from harmonica import prism_magnetic -from harmonica._forward.prism_magnetic import VALID_FIELDS +from harmonica._forward.prisms.magnetic import VALID_FIELDS from .utils import run_only_with_numba From a125cd1dd78ed41680187ab6d885632084c46613 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 16:00:22 -0700 Subject: [PATCH 03/13] Fix typo --- src/harmonica/_forward/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/harmonica/_forward/utils.py b/src/harmonica/_forward/utils.py index 1831a5f63..aa1b46264 100644 --- a/src/harmonica/_forward/utils.py +++ b/src/harmonica/_forward/utils.py @@ -349,7 +349,7 @@ def initialize_progressbar(total, use_progressbar): Number of iterations that will be assumed as the total for the :class:`numba_progress.ProgressBar`. use_progressbar : bool - Weather to initialize a progressbar or not. If True, then the function + Whether to initialize a progressbar or not. If True, then the function will return None. Returns From 15796e858f347de7910942278a3b763604cbd7b2 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 16:00:51 -0700 Subject: [PATCH 04/13] Refactor `gravity` method of the prism layer Reimplement the forward modeling computation of the gravity fields of a prism layer through a custom Numba function that iterates directly over the 1D horizontal coordinates of the layer. Replace the old implementation that used to build the whole set of prism boundaries and call the `prism_gravity` function. This new implementation doesn't require to store all those prisms in memory, significantly lowering the memory requirements, especially for large problems. --- src/harmonica/_forward/prisms/layer.py | 191 ++++++++++++++++++++----- 1 file changed, 158 insertions(+), 33 deletions(-) diff --git a/src/harmonica/_forward/prisms/layer.py b/src/harmonica/_forward/prisms/layer.py index f75e0ce37..959fdcc53 100644 --- a/src/harmonica/_forward/prisms/layer.py +++ b/src/harmonica/_forward/prisms/layer.py @@ -10,12 +10,14 @@ import warnings +import numba import numpy as np import verde as vd import xarray as xr from ...visualization import prism_to_pyvista -from .gravity import prism_gravity +from ..utils import initialize_progressbar +from .gravity import FIELDS, prism_gravity def prism_layer( @@ -313,23 +315,23 @@ def gravity( self, coordinates, field, - progressbar=False, + *, density_name="density", thickness_threshold=None, - **kwargs, + parallel=True, + progressbar=False, ): r""" Compute the gravity generated by the layer of prisms. - Uses :func:`harmonica.prism_gravity` for computing the gravity field - generated by the prisms of the layer. + Compute the gravity field generated by the layer of prisms on a set of + observation points. The density of the prisms will be assigned from the ``data_var`` chosen through the ``density_name`` argument. Ignores the prisms which ``top`` or ``bottom`` boundaries are ``np.nan``'s. Prisms thinner than a given threshold can be optionally ignored through the ``thickness_threshold`` argument. - All ``kwargs`` will be passed to :func:`harmonica.prism_gravity`. Parameters ---------- @@ -346,18 +348,20 @@ def gravity( - Downward acceleration: ``g_z`` - Diagonal tensor components: ``g_ee``, ``g_nn``, ``g_zz`` - Non-diagonal tensor components: ``g_en``, ``g_ez``, ``g_nz`` - progressbar : bool (optional) - If True, a progress bar of the computation will be printed to - standard error (stderr). Requires :mod:`numba_progress` to be - installed. Default to ``False``. - density_name : str (optional) + density_name : str, optional Name of the property layer (or ``data_var`` of the :class:`xarray.Dataset`) that will be used for the density of each prism in the layer. Default to ``"density"`` - thickness_threshold : float or None + thickness_threshold : float or None, optional Prisms thinner than this threshold will be ignored in the forward gravity calculation. If None, every prism with non-zero volume will be considered. Default to None. + parallel : bool, optional + Whether to run the computation of gravity fields in parallel or in serial. + progressbar : bool, optional + If True, a progress bar of the computation will be printed to + standard error (stderr). Requires :mod:`numba_progress` to be + installed. Default to ``False``. Returns ------- @@ -369,31 +373,53 @@ def gravity( -------- harmonica.prism_gravity """ - # Get boundaries and density of the prisms - boundaries = self._to_prisms() + # Sanity check the field + if field not in FIELDS: + msg = f"Gravitational field {field} not recognized" + raise ValueError(msg) + + # Check if prism layer defines a regular grid in horizontal coords + _check_regular_grid(self._obj.easting.values, self._obj.northing.values) + + # Ravel coordinates to 1D + cast = np.broadcast(*coordinates[:3]) + coordinates = tuple(c.ravel() for c in coordinates) + + # Determine parallel or serial forward modelling function + numba_function = ( + _forward_gravity_prism_layer_parallel + if parallel + else _forward_gravity_prism_layer_serial + ) + + # Forward model the prism layer + thickness_threshold = ( + 0.0 if thickness_threshold is None else thickness_threshold + ) density = self._obj[density_name].values - # Get the mask for selecting only the prisms whose top boundary, bottom - # boundary and density have no nans - mask = self._get_nonans_mask(property_name=density_name) - # Select only the boundaries and density elements for masked prisms - boundaries = boundaries[mask.ravel()] - density = density[mask] - # Discard thin prisms and their densities - if thickness_threshold is not None: - boundaries, density = _discard_thin_prisms( - boundaries, + with initialize_progressbar(coordinates[0].size, progressbar) as progress_proxy: + result = numba_function( + coordinates, + self._obj.easting.values, + self._obj.northing.values, + self._obj.bottom.values, + self._obj.top.values, density, + FIELDS[field], thickness_threshold, + progress_proxy, ) - # Return gravity field of prisms - return prism_gravity( - coordinates, - prisms=boundaries, - density=density, - field=field, - progressbar=progressbar, - **kwargs, - ) + + # Invert sign of gravity_u, gravity_eu, gravity_nu + if field in ("g_z", "g_ez", "g_nz"): + result *= -1 + # Convert to more convenient units + if field in ("g_e", "g_n", "g_z"): + result *= 1e5 # SI to mGal + # Convert to more convenient units + if field in ("g_ee", "g_nn", "g_zz", "g_en", "g_ez", "g_nz"): + result *= 1e9 # SI to Eotvos + return result.reshape(cast.shape) def _get_nonans_mask(self, property_name=None): """ @@ -556,3 +582,102 @@ def _discard_thin_prisms( prisms = prisms[np.logical_not(null_prisms), :] density = density[np.logical_not(null_prisms)] return prisms, density + + +def _forward_gravity_prism_layer( + coordinates, + prisms_easting, + prisms_northing, + prisms_bottom, + prisms_top, + densities, + forward_func, + thickness_threshold=0.0, + progress_proxy=None, +): + """ + Forward model the gravity fields of prisms in a prism layer. + + Parameters + ---------- + coordinates : tuple of (n,) array + Tuple with the easting, northing, and upward coordinates of the observation + points. + prisms_easting : (n_e,) array + 1D array with the easting coordinates of the center of the prisms. + Use the easting dimensional coordinate of the prism layer. + prisms_northing : (n_n,) array + 1D array with the northing coordinates of the center of the prisms. + Use the northing dimensional coordinate of the prism layer. + prisms_bottom, prisms_top : (n_n, n_e) array + 2D arrays with the bottom and top boundaries of the prisms. + densities : (n_n, n_e) array + 2D array with the densities of the prisms in kg/m3. + forward_func : callable + Choclo function to forward model a gravity field of prisms. + thickness_threshold : float + Prisms thinner than this threshold will be ignored in the + forward gravity calculation. If None, every prism with non-zero + volume will be considered. + progress_proxy : :class:`numba_progress.ProgressBar` or None + Instance of :class:`numba_progress.ProgressBar` that gets updated after + each iteration on the observation points. Use None if no progress bar + is should be used. + + Returns + ------- + result : array + Gravity field in SI units. + """ + # Remember to: + # - ignore prisms with zero density + # - ignore prisms based on thickness threshold + + easting, northing, upward = coordinates + n_coords = easting.size + half_spacing_east = (prisms_easting[1] - prisms_easting[0]) / 2 + half_spacing_north = (prisms_northing[1] - prisms_northing[0]) / 2 + + result = np.zeros(n_coords, dtype=np.float64) + for i in numba.prange(n_coords): + for j, easting_center in enumerate(prisms_easting): + west = easting_center - half_spacing_east + east = easting_center + half_spacing_east + for k, northing_center in enumerate(prisms_northing): + # Ignore prisms with zero density + density = densities[k, j] + if density == 0.0: + continue + + # Ignore thin prisms + bottom, top = prisms_bottom[k, j], prisms_top[k, j] + if top - bottom < thickness_threshold: + continue + + # Ignore prisms with nan top or bottom + if np.isnan(top) or np.isnan(bottom): + continue + + south = northing_center - half_spacing_north + north = northing_center + half_spacing_north + result[i] += forward_func( + easting[i], + northing[i], + upward[i], + west, + east, + south, + north, + bottom, + top, + density, + ) + return result + + +_forward_gravity_prism_layer_parallel = numba.jit(nopython=True, parallel=True)( + _forward_gravity_prism_layer +) +_forward_gravity_prism_layer_serial = numba.jit(nopython=True, parallel=False)( + _forward_gravity_prism_layer +) From 62540e8fe39deef6c9e79e40d3ae2f89bedd88c5 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 16:12:28 -0700 Subject: [PATCH 05/13] Autoformat --- src/harmonica/_forward/prisms/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/harmonica/_forward/prisms/__init__.py b/src/harmonica/_forward/prisms/__init__.py index 23640b8dd..af872175b 100644 --- a/src/harmonica/_forward/prisms/__init__.py +++ b/src/harmonica/_forward/prisms/__init__.py @@ -9,5 +9,5 @@ """ from .gravity import prism_gravity -from .magnetic import prism_magnetic from .layer import DatasetAccessorPrismLayer, prism_layer +from .magnetic import prism_magnetic From d7bb19a1697aa331d11105b36d27f1e89bcb297e Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 16:15:09 -0700 Subject: [PATCH 06/13] Improve docs of numba function --- src/harmonica/_forward/prisms/layer.py | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/src/harmonica/_forward/prisms/layer.py b/src/harmonica/_forward/prisms/layer.py index 959fdcc53..24254f2cc 100644 --- a/src/harmonica/_forward/prisms/layer.py +++ b/src/harmonica/_forward/prisms/layer.py @@ -598,6 +598,17 @@ def _forward_gravity_prism_layer( """ Forward model the gravity fields of prisms in a prism layer. + Builds the boundaries of each one the prisms in the layer on the fly, iterating over + the easting and northing dimensional coordinates. This function is intended to avoid + allocating all prism boundaries in memory before computation. + + The function ignores prisms with zero density, prisms thinner than the + ``thickness_threshold``, and prisms with NaN as top or bottom boundaries. + + .. important:: + + This function is intended to be decorated with Numba before use. + Parameters ---------- coordinates : tuple of (n,) array @@ -629,10 +640,6 @@ def _forward_gravity_prism_layer( result : array Gravity field in SI units. """ - # Remember to: - # - ignore prisms with zero density - # - ignore prisms based on thickness threshold - easting, northing, upward = coordinates n_coords = easting.size half_spacing_east = (prisms_easting[1] - prisms_easting[0]) / 2 From d8d9d4660d8cdcf7ad7147a8ead5f9db7ea0c27a Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 16:33:30 -0700 Subject: [PATCH 07/13] Ignore prisms with density as NaN and raise warning --- src/harmonica/_forward/prisms/layer.py | 10 ++++++++-- test/test_prism_layer.py | 3 +-- 2 files changed, 9 insertions(+), 4 deletions(-) diff --git a/src/harmonica/_forward/prisms/layer.py b/src/harmonica/_forward/prisms/layer.py index 24254f2cc..a3001ce5d 100644 --- a/src/harmonica/_forward/prisms/layer.py +++ b/src/harmonica/_forward/prisms/layer.py @@ -397,6 +397,12 @@ def gravity( 0.0 if thickness_threshold is None else thickness_threshold ) density = self._obj[density_name].values + if np.isnan(density).any(): + msg = ( + "Found NaN values in 'density' property of the prisms layer. " + "Their respective prisms will be ignored." + ) + warnings.warn(msg, stacklevel=2) with initialize_progressbar(coordinates[0].size, progressbar) as progress_proxy: result = numba_function( coordinates, @@ -651,9 +657,9 @@ def _forward_gravity_prism_layer( west = easting_center - half_spacing_east east = easting_center + half_spacing_east for k, northing_center in enumerate(prisms_northing): - # Ignore prisms with zero density + # Ignore prisms with zero or NaN density density = densities[k, j] - if density == 0.0: + if density == 0.0 or np.isnan(density): continue # Ignore thin prisms diff --git a/test/test_prism_layer.py b/test/test_prism_layer.py index 0192fe655..d8d46d1f0 100644 --- a/test/test_prism_layer.py +++ b/test/test_prism_layer.py @@ -401,9 +401,8 @@ def test_prism_layer_gravity_density_nans(field, dummy_layer, prism_layer_with_h prisms_coords, surface, reference, properties={"density": density} ) # Check if warning is raised after passing density with nans - with warnings.catch_warnings(record=True) as warn: + with pytest.warns(): result = layer.prism_layer.gravity(coordinates, field=field) - assert len(warn) == 1 # Check if it generates the expected gravity field prisms, rho = prism_layer_with_holes npt.assert_allclose( From fb0f259600793ae5b406733735a43461e20f1ef4 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 16:34:10 -0700 Subject: [PATCH 08/13] Remove unused import --- src/harmonica/_forward/prisms/layer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/harmonica/_forward/prisms/layer.py b/src/harmonica/_forward/prisms/layer.py index a3001ce5d..fd53df113 100644 --- a/src/harmonica/_forward/prisms/layer.py +++ b/src/harmonica/_forward/prisms/layer.py @@ -17,7 +17,7 @@ from ...visualization import prism_to_pyvista from ..utils import initialize_progressbar -from .gravity import FIELDS, prism_gravity +from .gravity import FIELDS def prism_layer( From db7b626cf3f676a0557e7c0eef3b1167d09d007d Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 16:36:09 -0700 Subject: [PATCH 09/13] Update the progressbar if any --- src/harmonica/_forward/prisms/layer.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/harmonica/_forward/prisms/layer.py b/src/harmonica/_forward/prisms/layer.py index fd53df113..adf91c84b 100644 --- a/src/harmonica/_forward/prisms/layer.py +++ b/src/harmonica/_forward/prisms/layer.py @@ -651,6 +651,9 @@ def _forward_gravity_prism_layer( half_spacing_east = (prisms_easting[1] - prisms_easting[0]) / 2 half_spacing_north = (prisms_northing[1] - prisms_northing[0]) / 2 + # Check if we need to update the progressbar on each iteration + update_progressbar = progress_proxy is not None + result = np.zeros(n_coords, dtype=np.float64) for i in numba.prange(n_coords): for j, easting_center in enumerate(prisms_easting): @@ -685,6 +688,9 @@ def _forward_gravity_prism_layer( top, density, ) + # Update progress bar if called + if update_progressbar: + progress_proxy.update(1) return result From 7c6e315f4aa778627b97d03aeea85c07ce091d26 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 16:40:32 -0700 Subject: [PATCH 10/13] Remove obsolete private function and method --- src/harmonica/_forward/prisms/layer.py | 76 -------------------------- 1 file changed, 76 deletions(-) diff --git a/src/harmonica/_forward/prisms/layer.py b/src/harmonica/_forward/prisms/layer.py index adf91c84b..96d1c6e83 100644 --- a/src/harmonica/_forward/prisms/layer.py +++ b/src/harmonica/_forward/prisms/layer.py @@ -427,43 +427,6 @@ def gravity( result *= 1e9 # SI to Eotvos return result.reshape(cast.shape) - def _get_nonans_mask(self, property_name=None): - """ - Build a mask for prisms with no nans on top, bottom or a property. - - Parameters - ---------- - property_name : str (optional) - Name of the property layer (or ``data_var`` of the - :class:`xarray.Dataset`) that will be used for masking the prisms - in the layer. - - Returns - ------- - mask : 2d-array - Array of bools that can be used as a mask for selecting prisms with - no nans on top boundaries, bottom boundaries and the passed - property. - """ - # Mask the prisms that contains no nans on top and bottom boundaries - mask = np.logical_and( - np.logical_not(np.isnan(self._obj.top.values)), - np.logical_not(np.isnan(self._obj.bottom.values)), - ) - # Mask the prisms that contains nans on the selected property - if property_name is not None: - mask_property = np.logical_not(np.isnan(self._obj[property_name].values)) - # Warn if a nan is found within the masked property - if not mask_property[mask].all(): - warnings.warn( - f"Found missing values in '{property_name}' property " - + "of the prisms layer. The prisms with a nan as " - + f"'{property_name}' will be ignored.", - stacklevel=1, - ) - mask = np.logical_and(mask, mask_property) - return mask - def _to_prisms(self): """ Return the boundaries of each prism of the layer. @@ -551,45 +514,6 @@ def to_pyvista(self, drop_null_prisms=True): return prism_to_pyvista(prisms, properties=properties) -def _discard_thin_prisms( - prisms, - density, - thickness_threshold, -): - """ - Discard prisms with a thickness below a threshold. - - Parameters - ---------- - prisms : 2d-array - Array containing the boundaries of the prisms in the following order: - ``w``, ``e``, ``s``, ``n``, ``bottom``, ``top``. - The array must have the following shape: (``n_prisms``, 6), where - ``n_prisms`` is the total number of prisms. - density : 1d-array - Array containing the density of each prism in kg/m^3. Must have the - same size as the number of prisms. - thickness_threshold : float - Prisms thinner than this threshold will be discarded. - - Returns - ------- - prisms : 2d-array - A copy of the ``prisms`` array that doesn't include the thin prisms. - density : 1d-array - A copy of the ``density`` array that doesn't include the density values - for thin prisms. - """ - bottom, top = prisms[:, -2], prisms[:, -1] - # Mark prisms with thickness < threshold as null prisms - thickness = top - bottom - null_prisms = thickness < thickness_threshold - # Keep only thick prisms and their densities - prisms = prisms[np.logical_not(null_prisms), :] - density = density[np.logical_not(null_prisms)] - return prisms, density - - def _forward_gravity_prism_layer( coordinates, prisms_easting, From 71c77dc3ebfd7f7c11e504725a5b95f5728b067f Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 16:44:57 -0700 Subject: [PATCH 11/13] Remove tests for removed functions and methods --- test/test_prism_layer.py | 146 +-------------------------------------- 1 file changed, 3 insertions(+), 143 deletions(-) diff --git a/test/test_prism_layer.py b/test/test_prism_layer.py index d8d46d1f0..488c86a36 100644 --- a/test/test_prism_layer.py +++ b/test/test_prism_layer.py @@ -8,7 +8,7 @@ Test prisms layer. """ -import warnings +import re from unittest.mock import patch import numpy as np @@ -18,7 +18,6 @@ import xarray as xr from harmonica import prism_gravity, prism_layer -from harmonica._forward.prisms.layer import _discard_thin_prisms try: import pyvista @@ -238,107 +237,6 @@ def test_prism_layer_get_prism_by_index(): ) -def test_nonans_prisms_mask(dummy_layer): - """ - Check if the mask for nonans prism is correctly created. - """ - (easting, northing), surface, reference, _ = dummy_layer - shape = (northing.size, easting.size) - - # No nan in top nor bottom - # ------------------------ - layer = prism_layer((easting, northing), surface, reference) - expected_mask = np.ones(shape, dtype=bool) - mask = layer.prism_layer._get_nonans_mask() - npt.assert_allclose(mask, expected_mask) - - # Nans in top only - # ---------------- - layer = prism_layer((easting, northing), surface, reference) - expected_mask = np.ones(shape, dtype=bool) - for index in ((1, 2), (2, 3)): - layer.top[index] = np.nan - expected_mask[index] = False - mask = layer.prism_layer._get_nonans_mask() - npt.assert_allclose(mask, expected_mask) - - # Nans in bottom only - # ------------------- - layer = prism_layer((easting, northing), surface, reference) - expected_mask = np.ones(shape, dtype=bool) - for index in ((2, 1), (3, 2)): - layer.bottom[index] = np.nan - expected_mask[index] = False - mask = layer.prism_layer._get_nonans_mask() - npt.assert_allclose(mask, expected_mask) - - # Nans in top and bottom - # ---------------------- - layer = prism_layer((easting, northing), surface, reference) - expected_mask = np.ones(shape, dtype=bool) - for index in ((1, 2), (2, 3)): - layer.top[index] = np.nan - expected_mask[index] = False - for index in ((1, 2), (2, 1), (3, 2)): - layer.bottom[index] = np.nan - expected_mask[index] = False - mask = layer.prism_layer._get_nonans_mask() - npt.assert_allclose(mask, expected_mask) - - -def test_nonans_prisms_mask_property( - dummy_layer, -): - """ - Check if the method masks the property and raises a warning. - """ - (easting, northing), surface, reference, density = dummy_layer - shape = (northing.size, easting.size) - - # Nans in top and property (on the same prisms) - # --------------------------------------------- - expected_mask = np.ones_like(surface, dtype=bool) - indices = ((1, 2), (2, 3)) - # Set some elements of surface and density as nans - for index in indices: - surface[index] = np.nan - density[index] = np.nan - expected_mask[index] = False - layer = prism_layer( - (easting, northing), surface, reference, properties={"density": density} - ) - # Check if no warning is raised - with warnings.catch_warnings(record=True) as warn: - mask = layer.prism_layer._get_nonans_mask(property_name="density") - assert len(warn) == 0 - npt.assert_allclose(mask, expected_mask) - - # Nans in top and property (not precisely on the same prisms) - # ----------------------------------------------------------- - surface = np.arange(20, dtype=float).reshape(shape) - density = 2670 * np.ones_like(surface) - expected_mask = np.ones_like(surface, dtype=bool) - # Set some elements of surface as nans - indices = ((1, 2), (2, 3)) - for index in indices: - surface[index] = np.nan - expected_mask[index] = False - # Set a different set of elements of density as nans - indices = ((2, 2), (0, 1)) - for index in indices: - density[index] = np.nan - expected_mask[index] = False - layer = prism_layer( - (easting, northing), surface, reference, properties={"density": density} - ) - # Check if warning is raised - with warnings.catch_warnings(record=True) as warn: - mask = layer.prism_layer._get_nonans_mask(property_name="density") - assert len(warn) == 1 - assert issubclass(warn[-1].category, UserWarning) - npt.assert_allclose(mask, expected_mask) - - @pytest.mark.use_numba @pytest.mark.parametrize("field", ["potential", "g_z"]) def test_prism_layer_gravity(field, dummy_layer): @@ -401,7 +299,8 @@ def test_prism_layer_gravity_density_nans(field, dummy_layer, prism_layer_with_h prisms_coords, surface, reference, properties={"density": density} ) # Check if warning is raised after passing density with nans - with pytest.warns(): + msg = re.escape("Found NaN values in 'density' property of the prisms layer.") + with pytest.warns(match=msg): result = layer.prism_layer.gravity(coordinates, field=field) # Check if it generates the expected gravity field prisms, rho = prism_layer_with_holes @@ -543,42 +442,3 @@ def test_gravity_discarded_thin_prisms(dummy_layer): coordinates, field="g_z", thickness_threshold=5 ) npt.assert_allclose(gravity_manually_removed, gravity_threshold_removed) - - -def test_discard_thin_prisms(): - """ - Check if thin prisms are properly discarded. - """ - # create set of 4 prisms (west, east, south, north, bottom, top) - prism_boundaries = np.array( - [ - [-5000.0, 5000.0, -5000.0, 5000.0, 0.0, 55.1], - [5000.0, 15000.0, -5000.0, 5000.0, 0.0, 55.01], - [-5000.0, 5000.0, 5000.0, 15000.0, 0.0, 35.0], - [5000.0, 15000.0, 5000.0, 15000.0, 0.0, 84.0], - ] - ) - - # assign densities to each prism - densities = np.array([2306, 2122, 2190, 2069]) - - # drop prisms and respective densities thinner than 55.05 - # (2nd and 3rd prisms) - thick_prisms, thick_densities = _discard_thin_prisms( - prism_boundaries, - densities, - thickness_threshold=55.05, - ) - - # manually remove prisms and densities of thin prisms - expected_prisms = np.array( - [ - [-5000.0, 5000.0, -5000.0, 5000.0, 0.0, 55.1], - [5000.0, 15000.0, 5000.0, 15000.0, 0.0, 84.0], - ] - ) - expected_densities = np.array([2306, 2069]) - - # check the correct prisms and densities were discarded - npt.assert_allclose(expected_prisms, thick_prisms) - npt.assert_allclose(expected_densities, thick_densities) From e676ec817aded54beefcc06ef242dbed6f1fe212 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 16:59:59 -0700 Subject: [PATCH 12/13] Fix raveling of coordinates --- src/harmonica/_forward/prisms/layer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/harmonica/_forward/prisms/layer.py b/src/harmonica/_forward/prisms/layer.py index 96d1c6e83..64b431289 100644 --- a/src/harmonica/_forward/prisms/layer.py +++ b/src/harmonica/_forward/prisms/layer.py @@ -383,7 +383,7 @@ def gravity( # Ravel coordinates to 1D cast = np.broadcast(*coordinates[:3]) - coordinates = tuple(c.ravel() for c in coordinates) + coordinates = tuple(np.atleast_1d(c).ravel() for c in coordinates[:3]) # Determine parallel or serial forward modelling function numba_function = ( From 6a80fada8e20b92658f72f19a829eda858efa614 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 21 Apr 2026 17:09:27 -0700 Subject: [PATCH 13/13] Add tests to increase coverage --- src/harmonica/_forward/prisms/layer.py | 2 +- test/test_prism_layer.py | 22 +++++++++++++++++++++- 2 files changed, 22 insertions(+), 2 deletions(-) diff --git a/src/harmonica/_forward/prisms/layer.py b/src/harmonica/_forward/prisms/layer.py index 64b431289..fa72ab0c2 100644 --- a/src/harmonica/_forward/prisms/layer.py +++ b/src/harmonica/_forward/prisms/layer.py @@ -375,7 +375,7 @@ def gravity( """ # Sanity check the field if field not in FIELDS: - msg = f"Gravitational field {field} not recognized" + msg = f"Gravitational field '{field}' not recognized." raise ValueError(msg) # Check if prism layer defines a regular grid in horizontal coords diff --git a/test/test_prism_layer.py b/test/test_prism_layer.py index 488c86a36..e6f2ab939 100644 --- a/test/test_prism_layer.py +++ b/test/test_prism_layer.py @@ -11,6 +11,7 @@ import re from unittest.mock import patch +import bordado as bd import numpy as np import numpy.testing as npt import pytest @@ -18,6 +19,7 @@ import xarray as xr from harmonica import prism_gravity, prism_layer +from harmonica._forward.prisms.gravity import FIELDS try: import pyvista @@ -238,7 +240,7 @@ def test_prism_layer_get_prism_by_index(): @pytest.mark.use_numba -@pytest.mark.parametrize("field", ["potential", "g_z"]) +@pytest.mark.parametrize("field", FIELDS) def test_prism_layer_gravity(field, dummy_layer): """ Check if gravity method works as expected. @@ -259,6 +261,24 @@ def test_prism_layer_gravity(field, dummy_layer): ) +def test_prism_layer_invalid_field(dummy_layer): + """ + Test error after passing invalid field. + """ + coordinates = bd.grid_coordinates( + (1, 3, 7, 10), spacing=1, non_dimensional_coords=30.0 + ) + (easting, northing), surface, reference, density = dummy_layer + layer = prism_layer( + (easting, northing), surface, reference, properties={"density": density} + ) + + invalid_field = "invalid field" + msg = re.escape(f"Gravitational field '{invalid_field}' not recognized.") + with pytest.raises(ValueError, match=msg): + layer.prism_layer.gravity(coordinates, field=invalid_field) + + @pytest.mark.use_numba @pytest.mark.parametrize("field", ["potential", "g_z"]) def test_prism_layer_gravity_surface_nans(field, dummy_layer, prism_layer_with_holes):