diff --git a/harmonica/_forward/_tesseroid_utils.py b/harmonica/_forward/_tesseroid_utils.py index eaa0b5210..2d73bdb07 100644 --- a/harmonica/_forward/_tesseroid_utils.py +++ b/harmonica/_forward/_tesseroid_utils.py @@ -402,7 +402,7 @@ def _check_tesseroids(tesseroids): return tesseroids -def _check_points_outside_tesseroids(coordinates, tesseroids): +def check_points_outside_tesseroids(coordinates, tesseroids): """ Check if computation points are not inside the tesseroids @@ -423,44 +423,53 @@ def _check_points_outside_tesseroids(coordinates, tesseroids): This array of tesseroids must have longitude continuity and valid boundaries. Run ``_check_tesseroids`` before. + + Raises + ------ + ValueError + If any computation point falls inside any tesseroid. """ - longitude, latitude, radius = coordinates[:] - west, east, south, north, bottom, top = tuple(tesseroids[:, i] for i in range(6)) - # Longitudinal boundaries of the tesseroid must be compared with - # longitudinal coordinates of computation points when moved to - # [0, 360) and [-180, 180). - longitude_360 = longitude % 360 - longitude_180 = ((longitude + 180) % 360) - 180 - inside_longitude = np.logical_or( - np.logical_and( - west < longitude_360[:, np.newaxis], longitude_360[:, np.newaxis] < east - ), - np.logical_and( - west < longitude_180[:, np.newaxis], longitude_180[:, np.newaxis] < east - ), - ) - inside_latitude = np.logical_and( - south < latitude[:, np.newaxis], latitude[:, np.newaxis] < north - ) - inside_radius = np.logical_and( - bottom < radius[:, np.newaxis], radius[:, np.newaxis] < top - ) - # Build array of booleans. - # The (i, j) element is True if the computation point i is inside the - # tesseroid j. - inside = inside_longitude * inside_latitude * inside_radius - if inside.any(): + longitude, latitude, radius = coordinates + conflicting = _check_points_outside_tesseroids(coordinates, tesseroids) + if conflicting: err_msg = ( - "Found computation point inside tesseroid. " - + "Computation points must be outside of tesseroids.\n" + "Found computation point(s) inside tesseroid(s). " + "Computation points must be outside of tesseroids.\n" ) - for point_i, tess_i in np.argwhere(inside): - err_msg += "\tComputation point '{}' found inside tesseroid '{}'\n".format( - coordinates[:, point_i], tesseroids[tess_i, :] + for i, j in conflicting: + west, east, south, north, bottom, top = tesseroids[j, :] + err_msg += ( + f" - Computation point '({longitude[i]}, {latitude[i]}, {radius[i]})' " + "inside tesseroid " + f"'({west}, {east}, {south}, {north}, {bottom}, {top})'.\n" ) raise ValueError(err_msg) +@jit(nopython=True) +def _check_points_outside_tesseroids(coordinates, tesseroids): + """ + Check if observation points fall inside tesseroids. + """ + longitude, latitude, radius = coordinates[:] + conflicting = [] + for i in range(longitude.size): + for j in range(tesseroids.shape[0]): + # Longitudinal boundaries of the tesseroid must be compared with + # longitudinal coordinates of computation points when moved to + # [0, 360) and [-180, 180). + longitude_360 = longitude[i] % 360 + longitude_180 = ((longitude[i] + 180) % 360) - 180 + west, east, south, north, bottom, top = tesseroids[j, :] + if ( + (west < longitude_180 < east or west < longitude_360 < east) + and south < latitude[i] < north + and bottom < radius[i] < top + ): + conflicting.append((i, j)) + return conflicting + + def _longitude_continuity(tesseroids): """ Modify longitudinal boundaries of tesseroids to ensure longitude continuity diff --git a/harmonica/_forward/tesseroid.py b/harmonica/_forward/tesseroid.py index 4154f2e29..3472695c9 100644 --- a/harmonica/_forward/tesseroid.py +++ b/harmonica/_forward/tesseroid.py @@ -13,9 +13,9 @@ from ..constants import GRAVITATIONAL_CONST from ._tesseroid_utils import ( _adaptive_discretization, - _check_points_outside_tesseroids, _check_tesseroids, _discard_null_tesseroids, + check_points_outside_tesseroids, gauss_legendre_quadrature, glq_nodes_weights, ) @@ -163,7 +163,7 @@ def tesseroid_gravity( # Sanity checks for tesseroids and computation points if not disable_checks: tesseroids = _check_tesseroids(tesseroids) - _check_points_outside_tesseroids(coordinates, tesseroids) + check_points_outside_tesseroids(coordinates, tesseroids) # Check if density is a function or constant values if callable(density): # Run density-based discretization on each tesseroid diff --git a/harmonica/tests/test_tesseroid.py b/harmonica/tests/test_tesseroid.py index 91b064f16..e4ab0711e 100644 --- a/harmonica/tests/test_tesseroid.py +++ b/harmonica/tests/test_tesseroid.py @@ -24,8 +24,8 @@ MAX_DISCRETIZATIONS, STACK_SIZE, _adaptive_discretization, - _check_points_outside_tesseroids, _check_tesseroids, + check_points_outside_tesseroids, tesseroid_gravity, ) from ..constants import GRAVITATIONAL_CONST @@ -129,7 +129,7 @@ def test_valid_tesseroid(): _check_tesseroids(np.atleast_2d([w, e, s, n, bottom, bottom])) # Check if valid tesseroid with west > east is not caught _check_tesseroids(np.atleast_2d([350, 10, s, n, bottom, top])) - # Check if tesseroids on [-180, 180) are not caught + # Check if tesseroids on -180, 180 are not caught _check_tesseroids(np.atleast_2d([-70, -60, s, n, bottom, top])) _check_tesseroids(np.atleast_2d([-10, 10, s, n, bottom, top])) _check_tesseroids(np.atleast_2d([-150, 150, s, n, bottom, top])) @@ -207,35 +207,71 @@ def test_disable_checks(): npt.assert_allclose(invalid_result, -valid_result) -def test_point_inside_tesseroid(): - "Check if a computation point inside the tesseroid is caught" - tesseroids = np.atleast_2d([-10, 10, -10, 10, 100, 200]) - # Test if outside point is not caught - points = [ - np.atleast_2d([0, 0, 250]).T, # outside point on radius - np.atleast_2d([20, 0, 150]).T, # outside point on longitude - np.atleast_2d([0, 20, 150]).T, # outside point on latitude - np.atleast_2d([0, 0, 200]).T, # point on top surface - np.atleast_2d([0, 0, 100]).T, # point on bottom surface - np.atleast_2d([-10, 0, 150]).T, # point on western surface - np.atleast_2d([10, 0, 150]).T, # point on eastern surface - np.atleast_2d([0, -10, 150]).T, # point on southern surface - np.atleast_2d([0, 10, 150]).T, # point on northern surface - ] - for coordinates in points: - _check_points_outside_tesseroids(coordinates, tesseroids) - # Test if computation point is inside the tesseroid - coordinates = np.atleast_2d([0, 0, 150]).T - with pytest.raises(ValueError): - _check_points_outside_tesseroids(coordinates, tesseroids) - # Test if computation point with phased longitude is inside the tesseroid - coordinates = np.atleast_2d([360, 0, 150]).T - with pytest.raises(ValueError): - _check_points_outside_tesseroids(coordinates, tesseroids) - tesseroids = np.atleast_2d([260, 280, -10, 10, 100, 200]) - coordinates = np.atleast_2d([-90, 0, 150]).T - with pytest.raises(ValueError): - _check_points_outside_tesseroids(coordinates, tesseroids) +class TestPointInsideTesseroid: + """Test checks for points inside tesseroids.""" + + @pytest.mark.use_numba + def test_points_outside_tesseroid(self): + """Check if error is not raised when no point is inside a tesseroid.""" + points = np.array( + [ + [0, 0, 250], # outside point on radius + [20, 0, 150], # outside point on longitude + [0, 20, 150], # outside point on latitude + [0, 0, 200], # point on top surface + [0, 0, 100], # point on bottom surface + [-10, 0, 150], # point on western surface + [10, 0, 150], # point on eastern surface + [0, -10, 150], # point on southern surface + [0, 10, 150], # point on northern surface + ] + ).T + tesseroid = np.atleast_2d([-10, 10, -10, 10, 100, 200]) + check_points_outside_tesseroids(points, tesseroid) + + @pytest.mark.use_numba + @pytest.mark.parametrize("points", [[0, 0, 150], [360, 0, 150]]) + def test_point_inside_tesseroid(self, points): + """Check if error is raised when point fall inside the tesseroid.""" + points = np.atleast_2d(points).T + tesseroid = np.atleast_2d([-10, 10, -10, 10, 100, 200]) + with pytest.raises(ValueError): + check_points_outside_tesseroids(points, tesseroid) + + @pytest.mark.use_numba + @pytest.mark.parametrize( + "point, tesseroid", + [ + ([360, 0, 150], [-10, 10, -10, 10, 100, 200]), + ([-90, 0, 150], [260, 280, -10, 10, 100, 200]), + ], + ) + def test_point_phased_longitude(self, point, tesseroid): + """Test if error is raised when the longitude coords are phased.""" + point = np.atleast_2d(point).T + tesseroid = np.atleast_2d(tesseroid) + with pytest.raises(ValueError): + check_points_outside_tesseroids(point, tesseroid) + + @pytest.mark.use_numba + def test_multiple_points_and_tesseroids(self): + """Check if error is raised with multiple points and tesseroids.""" + tesseroids = np.atleast_2d( + [ + [-10, 10, -10, 10, 100, 200], + [20, 30, 20, 30, 400, 500], + [-50, -40, -30, -20, 100, 500], + ] + ) + points = np.array( + [ + [0, 0, 150], + [80, 82, 4000], + [10, 10, 450], + ] + ).T + with pytest.raises(ValueError): + check_points_outside_tesseroids(points, tesseroids) @pytest.mark.use_numba