GDA94 to GDA2020 horizontal shift fails when using EGM2008 & EGM96 vertical crs about proj

[rouault]

PROJ is behaving here as best as it can given its heuristics and available data. When using EGM96 or EGM2008, you are using geoid models referenced to WGS 84, so you're GDA94 <--> GDA2020 transformation, also involves going through WGS 84. But both GDA94 and GDA2020 are considered as being equivalent to WGS 84 (with a ~ 3m inaccuracy), hence no horizontal shift

$ echo "400500 6500000 60" | PROJ_NETWORK=ON PROJ_DEBUG=2 bin/cs2cs --3d EPSG:28350 +to EPSG:7850+3855 -f "%.3f"
Using coordinate operation Inverse of Map Grid of Australia zone 50 + WGS 84 to EGM2008 height (1) using GDA2020 to WGS 84 (2) + Map Grid of Australia zone 50

and "GDA2020 to WGS 84 (2)" (EPSG:8450) is a no-op transformation, that has the remark "Approximation at the 3m level assuming WGS 84 is equivalent to GDA2020"

For what you want to do, you would have to probably do in 3 separate steps:

• GDA94 to GDA2020 3D
• then lie a bit by chaining with another transformation from WGS 84 3D to WGS 84 + EGM2008
• lie a bit again by chaining with "GDA2020 + EGM2008 height" to "GDA2020 / MGA zone 50 + EGM2008 height"

$ echo "400500 6500000 60" | PROJ_NETWORK=ON bin/cs2cs -d 9 --3d EPSG:28350 +to GDA2020 | PROJ_NETWORK=ON bin/cs2cs -d 9 --3d "WGS 84" "WGS 84 + EGM2008 height" | bin/cs2cs -d 3 "GDA2020 + EGM2008 height" "GDA2020 / MGA zone 50 + EGM2008 height"
400500.998	6500001.495 91.965

[scottmcnab]

Thanks for the explanation. That makes sense given the peculiar way GDA2020 is defined with respect to GDA94.

However, it seems to me in many (most?) cases horizontal accuracy may be preferrable over vertical accuracy. I was wondering if (as a future option perhaps?) there could be some way to specify "prioritise horizontal accuracy" in the heuristics algorithm? This way PROJ could automatically construct a pipeline similar to the steps you described above.

I don't know how such an algorithm might work, but I'm guessing something like "find best horizontal transform available" which would find the GDA2020 shift, even though WGS 84 transforms are also defined. Then "find suitable vertical transforms" for each side of the horizontal transform. This could be implemented as a future command-line option similar to the --3d flag to preserve backwards compatibility.

While it is possible to work around it manually for this particular case by being aware of it, I wonder if there are other CRS combinations where similar horizontal shifts are being missed? It would be good to have an automated solution that catches all of them.

[rouault]

If you look at the output of "projinfo -s EPSG:28350 -t EPSG:7850+3855 --3d", you'll find alternative pipelines than the one selected by default which will give better results.

For example if you select the second pipeline it proposes, you'll get:

$ echo "400500 6500000 60" | PROJ_NETWORK=ON bin/cct +proj=pipeline +step +inv +proj=utm +zone=50 +south +ellps=GRS80 +step +inv +proj=vgridshift +grids=us_nga_egm08_25.tif +multiplier=1 +step +proj=hgridshift +grids=au_icsm_GDA94_GDA2020_conformal_and_distortion.tif +step +proj=utm +zone=50 +south +ellps=GRS80
  400501.0015   6500001.4911       92.0437           inf

All of the first pipelines listed have the same computed accuracy, hence the algorithm doesn't really know which one is "best". The algorithm to compute pipelies is extremely complicated. I've fully authored it, but I can't hardly predict the results it will give, and it is very difficult to further tune without breaking other use cases. You're here in a quite advanced / not-really-geodetic-safe use case where you mix many datums: GDA94, GDA2020 and WGS84 (by asking for a EGM2008 height which is defined against WGS84), and depending on context one might consider that GDA94==WGS84 and GDA2020==WGS84, or that GDA94<->WGS84 should actually be understood as GDA94<->GDA2020, since both options are actually listed in EPSG. This is messy