From c957d4746e2ec929e0073dbaff8fcaab60559873 Mon Sep 17 00:00:00 2001 From: "Jean Pierre Mandujano G." Date: Fri, 24 Jul 2026 11:37:53 -0500 Subject: [PATCH 1/2] Add terrain correction in spherical coordinates to User Guide Add a section to the topographic correction guide showing how to forward model the topographic masses with a tesseroid layer in geocentric spherical coordinates: build the layer from the topography grid using radii, express the observation points consistently with the layer, clamp their radii so every computation point falls outside the tesseroids, and compute the topography-free gravity disturbance. --- doc/user_guide/topographic_correction.rst | 128 ++++++++++++++++++++++ 1 file changed, 128 insertions(+) diff --git a/doc/user_guide/topographic_correction.rst b/doc/user_guide/topographic_correction.rst index c0cb0bce0..4bbc65c63 100644 --- a/doc/user_guide/topographic_correction.rst +++ b/doc/user_guide/topographic_correction.rst @@ -319,6 +319,134 @@ topography. fig.colorbar(cmap=True, frame=["af", "x+lTopography", "y+lmeters"]) fig.show() + +Terrain correction in spherical coordinates +------------------------------------------- + +So far we computed the terrain effect by projecting the topography grid and +the observation points to plain Cartesian coordinates and approximating the +topographic masses with rectangular prisms. +On regional to global scales the curvature of the Earth cannot be neglected: +the projection distorts the geometry of the topographic masses and the +computed terrain effect accumulates errors. +In such cases we can forward model the topographic masses directly in +geocentric spherical coordinates using tesseroids (spherical prisms), which +take the curvature of the Earth into account. + +We can build a model of the topographic masses through the +:func:`harmonica.tesseroid_layer` function. +Unlike :func:`harmonica.prism_layer`, its ``surface`` and ``reference`` +arguments must be passed as **radii** measured from the center of the Earth, +not as heights above a reference level. +We can obtain the radii of the surface of the reference ellipsoid at each +latitude with :meth:`boule.Ellipsoid.geocentric_radius` and add the +topographic heights to them: + +.. jupyter-execute:: + + import boule as bl + + ellipsoid = bl.WGS84 + + longitude, latitude = np.meshgrid(topography.longitude, topography.latitude) + reference = ellipsoid.geocentric_radius(latitude) + surface = reference + topography.values + +We will assign the same densities we used for the layer of prisms and define +the layer of tesseroids: + +.. jupyter-execute:: + + density = np.where(topography.values >= 0, 2670, 1040 - 2670) + + tesseroids = hm.tesseroid_layer( + coordinates=(topography.longitude, topography.latitude), + surface=surface, + reference=reference, + properties={"density": density}, + ) + tesseroids + +The radial coordinate of the observation points must be expressed in the same +way as the boundaries of the layer: as radii from the center of the Earth. +We will compute them the same way we defined the ``surface`` of the layer, by +adding the observation heights to the geocentric radius of the ellipsoid at +each latitude. +This keeps the observation points consistent with the model of the topographic +masses: + +.. jupyter-execute:: + + radius = ellipsoid.geocentric_radius(data.latitude) + data.height_geometric_m + +Tesseroid forward modelling requires every computation point to be located +outside of the tesseroids. +Since our observations were taken on the terrain surface, some of them can +fall slightly below the top boundary of the tesseroid that contains them: the +tops of the tesseroids are given by the topography grid, whose values don't +exactly coincide with the observation heights. +We can make sure every observation point is located on or above the top of its +tesseroid by clamping their radii. +We will clamp against the highest top among the neighboring tesseroids, so +that observation points falling exactly on the boundary between two tesseroids +are safely lifted as well: + +.. jupyter-execute:: + + top = ( + tesseroids.top.rolling(longitude=3, latitude=3, center=True, min_periods=1) + .max() + .sel( + longitude=xr.DataArray(data.longitude), + latitude=xr.DataArray(data.latitude), + method="nearest", + ) + .values + ) + radius = np.maximum(radius, top) + +Now we can compute the terrain effect through the +:meth:`harmonica.DatasetAccessorTesseroidLayer.gravity` method: + +.. jupyter-execute:: + + coordinates_sph = (data.longitude, data.latitude, radius) + terrain_effect_spherical = tesseroids.tesseroid_layer.gravity( + coordinates_sph, field="g_z" + ) + +And obtain a topography-free gravity disturbance that takes the curvature of +the Earth into account: + +.. jupyter-execute:: + + topo_free_disturbance_spherical = ( + data.gravity_disturbance_mgal - terrain_effect_spherical + ) + + cpt_lims = vd.minmax(topo_free_disturbance_spherical) + + fig = pygmt.Figure() + pygmt.makecpt(cmap="viridis", series=cpt_lims) + fig.plot( + x=data.longitude, + y=data.latitude, + fill=topo_free_disturbance_spherical, + cmap=True, + style="c3p", + projection="M15c", + frame=['ag', 'WSen+ggray'], + ) + fig.colorbar( + cmap=True, + frame=[ + "a50f25", + "x+lTopography-free gravity disturbance (tesseroids)", + "y+lmGal", + ], + ) + fig.show() + ---- .. grid:: 2 From 7b1cc91d4d8916eb85fc5148254815f57ccc34c8 Mon Sep 17 00:00:00 2001 From: "Jean Pierre Mandujano G." Date: Wed, 2 Sep 2026 19:34:49 -0500 Subject: [PATCH 2/2] Lower tesseroid tops instead of lifting stations; compare with prisms --- doc/user_guide/topographic_correction.rst | 119 ++++++++++++++++++---- 1 file changed, 99 insertions(+), 20 deletions(-) diff --git a/doc/user_guide/topographic_correction.rst b/doc/user_guide/topographic_correction.rst index 4bbc65c63..dea194b33 100644 --- a/doc/user_guide/topographic_correction.rst +++ b/doc/user_guide/topographic_correction.rst @@ -367,6 +367,16 @@ the layer of tesseroids: ) tesseroids +.. note:: + + We are using the geodetic latitude of the topography grid as the latitude + of the tesseroids, which live in geocentric spherical coordinates. + This assumes the difference between the two latitudes (up to 0.2 degrees) + has no significant effect on the terrain correction. + Converting the grid to geocentric spherical coordinates would avoid the + assumption, but a regular grid in geodetic coordinates is not regular in + spherical ones, so the topography would have to be regridded first. + The radial coordinate of the observation points must be expressed in the same way as the boundaries of the layer: as radii from the center of the Earth. We will compute them the same way we defined the ``surface`` of the layer, by @@ -381,29 +391,42 @@ masses: Tesseroid forward modelling requires every computation point to be located outside of the tesseroids. -Since our observations were taken on the terrain surface, some of them can -fall slightly below the top boundary of the tesseroid that contains them: the -tops of the tesseroids are given by the topography grid, whose values don't -exactly coincide with the observation heights. -We can make sure every observation point is located on or above the top of its -tesseroid by clamping their radii. -We will clamp against the highest top among the neighboring tesseroids, so -that observation points falling exactly on the boundary between two tesseroids -are safely lifted as well: +Since our observations were taken on the terrain surface, some of them fall +below the top of the tesseroid that contains them: the tops of the tesseroids +are given by the topography grid, which averages the terrain over each cell, +while the observation heights were measured at each station. +Rather than moving the observation points, we will trust the measured heights +and lower the top of every tesseroid that contains a station below it to the +radius of that station. +A station that sits exactly on the boundary between two tesseroids belongs to +both, so we look up the tesseroids on every side of each station (shifting its +coordinates by far less than their precision): .. jupyter-execute:: - top = ( - tesseroids.top.rolling(longitude=3, latitude=3, center=True, min_periods=1) - .max() - .sel( - longitude=xr.DataArray(data.longitude), - latitude=xr.DataArray(data.latitude), - method="nearest", - ) - .values - ) - radius = np.maximum(radius, top) + indices = np.arange(tesseroids.top.size).reshape(tesseroids.top.shape) + cells = tesseroids.top.copy(data=indices) + lowest_station = np.full(tesseroids.top.size, np.inf) + shift = 1e-9 + for shift_longitude in (-shift, shift): + for shift_latitude in (-shift, shift): + index = cells.sel( + longitude=xr.DataArray(data.longitude + shift_longitude), + latitude=xr.DataArray(data.latitude + shift_latitude), + method="nearest", + ) + np.minimum.at(lowest_station, index.values, radius) + lowest_station = lowest_station.reshape(tesseroids.top.shape) + + surface = np.minimum(surface, lowest_station) + tesseroids.tesseroid_layer.update_top_bottom(surface, reference) + +.. note:: + + The same situation arises with the layer of prisms, but the prism forward + model doesn't require the computation points to be outside of the prisms, + so it went unnoticed in the previous section: the mass above those + stations is still part of that model. Now we can compute the terrain effect through the :meth:`harmonica.DatasetAccessorTesseroidLayer.gravity` method: @@ -447,6 +470,62 @@ the Earth into account: ) fig.show() +Compare the terrain effects of prisms and tesseroids +---------------------------------------------------- + +Even though this region spans only a few degrees, the two models don't agree. +Let's plot the difference between the terrain effects computed with prisms and +with tesseroids: + +.. jupyter-execute:: + + difference = terrain_effect - terrain_effect_spherical + + cpt_lims = vd.minmax(difference, min_percentile=5, max_percentile=95) + + fig = pygmt.Figure() + pygmt.makecpt(cmap="viridis", series=cpt_lims) + fig.plot( + x=data.longitude, + y=data.latitude, + fill=difference, + cmap=True, + style="c3p", + projection="M15c", + frame=['ag', 'WSen+ggray'], + ) + fig.colorbar( + cmap=True, + frame=["af", "x+lTerrain effect difference (prisms - tesseroids)", "y+lmGal"], + ) + fig.show() + +The tesseroids produce a terrain effect that is systematically larger, by +about 4.5 mGal on average, than the one produced by the prisms. +This is the effect of the curvature of the Earth on the terrain correction: the +topographic masses far from an observation point lie below the plane that is +tangent to the Earth at that point, so they pull more strongly downwards than +the same masses laid flat in a Cartesian model. +Unlike the difference between the Bouguer and the topography-free +disturbances, this one is fairly uniform: it depends on how much topography +surrounds each station rather than on how rugged it is, so it grows with the +extent of the topography grid and shrinks towards its edges. +For this grid, which extends a few hundred kilometers around the observations, +it amounts to 3-4% of the terrain effect. +Whether that is negligible depends on the goal of the survey: it's comparable +to the differences we found above between the Bouguer and the topography-free +disturbances. + +.. hint:: + + This is the same effect that the classic Bullard B (curvature) correction + accounts for when applying a Bouguer correction with a spherical cap + instead of an infinite slab. + +The largest differences are found at stations that lie below the topography +grid. The prism model still has topographic mass above them, which pulls +upwards, while we removed it from the tesseroid model. + ---- .. grid:: 2