OSGeo / PROJ

PROJ - Cartographic Projections and Coordinate Transformations Library
https://proj.org
Other
1.75k stars 788 forks source link

Improving geographic <--> vertical ballpark transformations by using a "generic" geoid model ? #1743

Open rouault opened 4 years ago

rouault commented 4 years ago

Something I've had in a corner of my head for some time...

We have from the EPSG dataset a number of registered tranformations between a geographic 3D CRS and a vertical CRS (NAD83 to NAVD88, etc...), but there are far more vertical CRS in the EPSG dataset for which there is no transformation available. Currently if converting between a compound CRS using such vertical CRS and a geographic 3D CRS, we will just return the Z value unchanged, using the "Ballpark vertical transformation". I'm wondering if it wouldn't make more sense to use a default model, like a EGM 96 or EGM 2008 (only for Earth-related transformations...), while still advertizing it as "Ballpark vertical transformation").

There are pros/cons: current approach gives really wrong result, but it is quite obvious when the user compares its input and output that no appropriate vertical done, so he knows that is wrong. If we apply a default model, we should get results that have a precision of one meter or so (when looking at all the geoid models they seem to only defer by the order of one meter, often less, sometimes a bit more, compared to EGM 2008), but the user of our "black box" APIs ( cs2cs, proj_create_crs_to_crs() ) will have a hard time figuring out it is just a best guess operation that has been applied.

One option would be to have a flag to proj_create_operations() to say "use default geoid model if there is nothing best available", but that wouldn't benefit to black box API style uses.

Ah, I should also mention that recently, I did a move in that direction for a particular (but common) case. Let me explain with an example. Let's say the EPSG dataset has a transformation "NAD83(NSRS2007) to NAVD88 height", but no "NAD83 to NAVD88" (it uses to have ones, until a recent update of the EPSG dataset, that have marked them as deprecated. here it is an instance where deprecated in EPSG doesn't mean broken... Not sure why they didn't just use the supersession mechanism). When transforming a compound CRS using NAVD88 height to a geographic 3D CRS, and none of source and target geodetic datums are the ones used by the vertical transformation, like NAD83 or WGS84, I will still use the transformation "something else to NAVD88 height" as a likely reasonable approximation for the actual desired transformation. Example between "NAD83 / UTM zone 11N + NAVD88 height" and "NAD83":

$ projinfo -s "EPSG:26911+5703" -t EPSG:4269 --spatial-test intersects --summary
Candidate operations found: 70
[...]
unknown id, Inverse of UTM zone 11N + NAD83 to NAD83(HARN) (4) + NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83(NSRS2007) to NAVD88 height (1) + Inverse of NAD83(HARN) to NAD83(NSRS2007) (1) + Inverse of NAD83 to NAD83(HARN) (4), 0.35 m, USA - California - south of 36.5°N, at least one grid missing
[...]

One can see that it does first a horizontal transformation from "NAD83 / UTM zone 11N + NAVD88 height" to NAD83(NSRS2007) (using NAD833(HARN) as an intermediate...), to be able to interpolate point from that grid, does the vertical adjustment, and then reconverts the horizontal part back to NAD83

kbevers commented 4 years ago

Just so I am sure that I undestand this correctly, this sitatution only comes up when doing the "augmented" EPSG codes like EPSG:26911+5703, right? That is, when users effectively create their own CRS made up of two official CRS's.

On one hand I can see how this is useful and gives users what they ask for in a very pleasing manner. On the other hand this is borderline geodetically unsound since we will be combining CRS's that probably/possibly was meant to never be used together and we will be claiming that a height is referred to a particular reference surface but in reality it is referring to something else. If we go down that route we need to put a significant uncertainty on the transformations as a whole.

kbevers commented 4 years ago

@busstoptaktik your input here would be appreciated

rouault commented 4 years ago

when doing the "augmented" EPSG codes like EPSG:26911+5703, right? That is, when users effectively create their own CRS made up of two official CRS's.

Yes. Note that users of PROJ don't necessarily always specifying themselves the CRS. Might come from some dataset they are reprojecting.

On the other hand this is borderline geodetically unsound since we will be combining CRS's that probably/possibly was meant to never be used together

Yes, as the "ballpark geographic transformation" we apply when converting between 2 geodetic datums that we can't bind together with a known transformation path.

If we go down that route we need to put a significant uncertainty on the transformations as a whole.

That would be the same as the current ballpark vertical transformation: set with a unknown accuracy, and sorted last among results

stale[bot] commented 4 years ago

This issue has been automatically marked as stale because it has not had recent activity. It will be closed if no further activity occurs. Thank you for your contributions.

busstoptaktik commented 4 years ago

pinned this, as the imminent arrival of EGM2020, and the prospect of a global vertical reference frame makes it worth taking this suggestion for another spin

rouault commented 4 years ago

There's an interesting exception to my initial comment. The recently integrated Austrian at_bev_GEOID_BESSEL_Oesterreich.tif grid between EPSG:9267 (MGI) and EPSG:9274 (EVRF2000 Austria height) is a "geoid-like" model, but the values in it are very different from all other geoid models on Austria.

See

$ gdallocationinfo -wgs84 at_bev_GEOID_BESSEL_Oesterreich.tif 13 47
Report:
  Location: (84P,83L)
  Band 1:
    Value: 2.36299991607666

vs

(geoid model between ETRS89 and EVRF2000 Austria height)

$ gdallocationinfo -wgs84 at_bev_GEOID_GRS80_Oesterreich.tif 13 47
Report:
  Location: (84P,83L)
  Band 1:
    Value: 50.1870002746582

or

(geoid model between ETRS89 and GHA Austria height)

Report:
  Location: (560P,492L)
  Band 1:
    Value: 50.2265625

or

$ gdallocationinfo -wgs84 us_nga_egm08_25.tif 13 47

Report:
  Location: (4632P,1032L)
  Band 1:
    Value: 50.3139915466309

This is presumably due to MGI being a ancient, pre-GPS era, datum, with an origin far from GRS80 center of mass, and so ellipsoidal heights on it are very different from the ones on WGS84 / ETRS89

busstoptaktik commented 2 years ago

What do you mean? Also it is not yet released, no?

@ValZapod, as I wrote, but you carefully edited out of your malquotation:

pinned this [issue], as the imminent arrival of EGM2020, and the prospect of a global vertical reference frame makes it worth taking this suggestion for another spin

What part of the words imminent and prospect are unclear here?

Improvements of software does not happen magically out of thin air - it helps to be prepared for the gifts the future may bring, so you can take advantage of them as early as possible.

kbevers commented 2 years ago

So are you going to implement XGM2019e ICGEM grid format

Probably not. You are welcome to write up a RFC that argues why it's a good idea to include in PROJ and of course implement it following a positive response. Adding support for another grid file format is another maintainance burden. We are not going to do that without a well-argued case in favour of it and someone that's willing to keep it alive for at least a number of years in the future.

rouault commented 2 years ago

Beyond the issue of handing a new format, I believe the challenge here with ICGEM format is that it is likely not ready for immediate/easy consumption by PROJ as it contains the spherical harmonic coefficients. There are likely computations to do as a post processing stage to get simple vertical shift values to apply. Similar exercice has been done I believe in https://geographiclib.sourceforge.io/ to derive the EGM96 & EGM2008 grids used by PROJ from their original delivery format.

MyqueWooMiddo commented 1 year ago

What do you mean? Also it is not yet released, no?

@ValZapod, as I wrote, but you carefully edited out of your malquotation:

pinned this [issue], as the imminent arrival of EGM2020, and the prospect of a global vertical reference frame makes it worth taking this suggestion for another spin

What part of the words imminent and prospect are unclear here?

Improvements of software does not happen magically out of thin air - it helps to be prepared for the gifts the future may bring, so you can take advantage of them as early as possible.

Where do you get EGM2020 resource ? Is that more accurate and resolution less than 1 minute?