Add terrain correction in spherical coordinates to the User Guide - #700
Add terrain correction in spherical coordinates to the User Guide#700manduinca wants to merge 2 commits into
Conversation
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.
santisoler
left a comment
santisoler
left a comment
There was a problem hiding this comment.
Thanks @manduinca for taking the time to work on this. This tutorial is looking really good!
On one side I'm thinking if it would be better to move this section into its own .rst document, as a subsection of the Topographic Correction. But on the other side, that would require repeating some of the code in the Topographic Correction (like downloading the data and the DEM, etc). Also, I think it might be nice to include a plot showing the differences between the terrain correction performed with prisms and tesseroids.
For this particular case where the region spans only a few degrees, I'm not sure we'll be able to see a significant difference between them. That would lead me to consider if it is better to use a different example with a larger region, where the curvature of the Earth introduces a significant difference in the results.
I guess I'm thinking while I type... But in summary, could you add an extra plot with the differences between terrain correction using prisms and tesseroids?
I left a few comments below as well. Let me know what do you think. And thanks again for the effort!
| 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) |
There was a problem hiding this comment.
I'm a bit confused about this piece of code. I'm not fully understanding why it's necessary to perform the windowing operation. Getting the closest top for each observation point should be enough to build the top array, right?
I'm also thinking if isn't better to adjust the top surface of the tesseroids rather than shifting the observation points. In gravity surveys the height measurements are usually very accurate, carried out through differential GNSS. While the DEMs are models that average the terrain elevation within a region given by the resolution of the model. What do you think?
There was a problem hiding this comment.
Update: I see your point on observation points that fall in the boundaries of tesseroids. My guess is that in those cases the points will never be internal to the tesseroid: they'll always lie on a face or an edge. So I'm not sure if those are going to be problematic.
| 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: | ||
|
|
There was a problem hiding this comment.
I would add an admonition here explaining that by assigning the ellipsoidal latitude as coordinates for the tesseroids, we are assuming that there's not significant difference between them and the geocentric spherical latitude.
One alternative way is to convert the ellipsoidal coordinates into geocentric spherical and build the tesseroids from them. For this particular tutorial I'd avoid doing so, because that introduces issues like regridding the topography in geocentric spherical, since a regular grid in ellipsoidal coordinates does not translate into a regular grid in geocentric spherical coordinates.
Let me know if you need any help with this.
|
Thanks for the detailed review @santisoler. All three points addressed; here is what I found while doing it. Comparison plot. Added, and the difference is visible even in this region: the tesseroids give a terrain effect about 4.5 mGal larger than the prisms, quite uniformly (rms 5.2 mGal over effects of 12 to 242 mGal, so 3-4%). To make sure it was the curvature and not the Mercator projection, I ran a check with constant topography of 1500 m on the same 1 arc-minute grid, one station at the center, growing the extent of the grid, with prisms projected both with Mercator and with a local azimuthal equidistant projection:
(mGal; the infinite Bouguer slab is 167.95.) The prisms converge to the slab, the tesseroids grow past it linearly with the extent of the grid, and the projection makes no difference. Interpolated to the Hayford radius (166.7 km) the difference is +1.55 mGal, which is LaFehr's Bullard B value for 1500 m, so I mention that in a hint. Given that, I kept this region rather than a larger example: the effect is already there and it avoids downloading a second dataset. Latitude note. Added as a Clamp. You were right to question it. The nearest top is enough for almost all stations, but not for the ones sitting exactly on a cell boundary, and this dataset has many of those: 40 stations in longitude and 36 in latitude have coordinates that are round fractions of a degree, right on the edges of the 1 arc-minute grid. The problem is floating point rather than geometry: the boundaries of neighboring tesseroids differ by ~7e-15 degrees, so a station at longitude 31.575 is inside the tesseroid whose east edge is 31.574999999999996 while I went with your suggestion instead: the measured heights stay, and the top of every tesseroid containing a station below it is lowered to the radius of that station (629 tesseroids, median 15 m, worst 392 m). To cover the boundary stations I look up the cell on either side of each station, shifting its coordinates by 1e-9 degrees, and lower both. No station is moved and the forward runs without errors. A note points out that the prism model has the same stations below the DEM, just without complaining. On splitting the section into its own page: I left it here because of the duplication you mentioned (download, cropping, projection), but happy to move it if you prefer. The whole page runs end-to-end for me (with PyGMT stubbed out, since I don't have GMT locally), about a minute in total. |
2 participants
Closes #353
Adds a section to the topographic correction guide showing how to forward model the topographic masses with a tesseroid layer in geocentric spherical coordinates, mirroring the existing Cartesian workflow: build the layer from the topography grid using radii, express the observation points consistently with the layer (geodetic latitude +
geocentric_radius(latitude) + height), compute the topography-free gravity disturbance, and compare the terrain effects of both models.Since the observations sit on the terrain, 619 of the 3877 stations fall below the top of their tesseroid (the DEM averages each cell; median deficit 15 m, worst 392 m), and
tesseroid_gravityrequires computation points outside the tesseroids. Following the review, the section keeps the measured heights and lowers the top of the tesseroids that contain those stations instead of moving the stations. A note explains the geodetic vs. geocentric latitude assumption.The comparison shows the tesseroids give a terrain effect about 4.5 mGal (3-4%) larger than the prisms, uniformly across the region. A check with constant topography (details in the comments) confirms it's the curvature of the Earth, not the projection: the difference grows linearly with the extent of the grid and matches the Bullard B correction.
I couldn't build the docs locally (missing GMT binary for PyGMT), but the code of the whole page runs end-to-end (~1 min, tesseroid forward included).