diff --git a/doc/user_guide/topographic_correction.rst b/doc/user_guide/topographic_correction.rst index c0cb0bce0..dea194b33 100644 --- a/doc/user_guide/topographic_correction.rst +++ b/doc/user_guide/topographic_correction.rst @@ -319,6 +319,213 @@ 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 + +.. 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 +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 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:: + + 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: + +.. 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() + +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