[PROJ] WGS84 to GDA94 and GDA2020

[Even Rouault]

Nigel,

Le 11/10/2024 à 12:27, Nigel Stewart via PROJ a écrit :
> > First of all I think you are being way too unrealistic about what 
> you can use ChatGPT for.
> >If you'd tried the command it's suggested to you, you would know that 
> it doesn't work.
> > It has given you an ill-formed CRS, as well as suggesting you use 
> deprecated syntax.
>
> It didn't smell good, I must admit.
>
> > It would be helpful if you can be more explicit about what you are 
> trying to do.
> > What EPSG-codes are you trying to get a transformation between?
>
> In a nutshell:
>
> # to GDA94 / MGA zone 56
>
> $ echo "0 0 0" | cs2cs +proj=tmerc +lat_0=-34.0 +lon_0=150.0 +k=1 
> +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs +to 'EPSG:28356'
> 222908.70       6233785.28 0.00
>
> # to GDA2020 / MGA zone 56
>
> $ echo "0 0 0" | cs2cs +proj=tmerc +lat_0=-34.0 +lon_0=150.0 +k=1 
> +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs +to 'EPSG:7856'
> 222908.70       6233785.28 0.00
>
> # to GDA2020 / MGA zone 56
>
> $ echo "222908.70       6233785.28 0.00" | cs2cs 'EPSG:28356' +to 
> 'EPSG:7856'
> 222909.17       6233786.71 0.00
>
> Our expectation is that GDA2020 is within a meter or two of GDA94 due 
> to the tectonic drift of the Australian plate towards the north east.
> Our concern is that from these custom transverse mercator projections, 
> this difference isn't reflected.
> We can see the Helmert transformation being applied in the third case, 
> but not in the first two.

This is fully expected.

Basically GDA94 ~= WGS84 at epoch 1994.0 and GDA2020 ~= WGS84 at epoch 
2020.0

So when you ask to transform between GDA94 and WGS84 a null datum 
transformation is a perfectly valid answer since they matched in 1994.0. 
As well as a null one between GDA2020 and WGS84 since they matched in 2020.0

There are also alternate transformations available in the EPSG dataset 
where they consider that when transforming between GDA94 and WGS84, then 
the user could mean they assume this WGS84 to be a recent one, so you 
actually want to use the transformation between GDA94 and GDA2020


$ projinfo -s "+proj=tmerc +lat_0=-34.0 +lon_0=150.0 +k=1 +x_0=0 +y_0=0 
+datum=WGS84 +type=crs" -t EPSG:28356 --spatial-test intersects -o PROJ
Candidate operations found: 3
-------------------------------------
Operation No. 1:

unknown id, Inverse of unknown + axis order change (2D) + Inverse of 
GDA94 to WGS 84 (1) + Map Grid of Australia zone 56, 3 m, Australia 
including Lord Howe Island, Macquarie Island, Ashmore and Cartier 
Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All 
onshore and offshore.

PROJ string:
+proj=pipeline
   +step +inv +proj=tmerc +lat_0=-34 +lon_0=150 +k=1 +x_0=0 +y_0=0 
+ellps=WGS84
   +step +proj=utm +zone=56 +south +ellps=GRS80

-------------------------------------
Operation No. 2:

unknown id, Inverse of unknown + axis order change (2D) + Inverse of 
GDA94 to WGS 84 (2) + Map Grid of Australia zone 56, 3 m, Australia 
including Lord Howe Island, Macquarie Island, Ashmore and Cartier 
Islands, Christmas Island, Cocos (Keeling) Islands, Norfolk Island. All 
onshore and offshore.

PROJ string:
+proj=pipeline
   +step +inv +proj=tmerc +lat_0=-34 +lon_0=150 +k=1 +x_0=0 +y_0=0 
+ellps=WGS84
   +step +proj=push +v_3
   +step +proj=cart +ellps=WGS84
   +step +inv +proj=helmert +x=0.06155 +y=-0.01087 +z=-0.04019 
+rx=-0.0394924
         +ry=-0.0327221 +rz=-0.0328979 +s=-0.009994 
+convention=coordinate_frame
   +step +inv +proj=cart +ellps=GRS80
   +step +proj=pop +v_3
   +step +proj=utm +zone=56 +south +ellps=GRS80

-------------------------------------
Operation No. 3:

unknown id, Inverse of unknown + axis order change (2D) + Inverse of 
GDA94 to WGS 84 (3) + Map Grid of Australia zone 56, 3 m, Australia - 
Australian Capital Territory; New South Wales; Northern Territory; 
Queensland; South Australia; Tasmania; Western Australia; Victoria., at 
least one grid missing

PROJ string:
+proj=pipeline
   +step +inv +proj=tmerc +lat_0=-34 +lon_0=150 +k=1 +x_0=0 +y_0=0 
+ellps=WGS84
   +step +inv +proj=hgridshift
         +grids=au_icsm_GDA94_GDA2020_conformal_and_distortion.tif
   +step +proj=utm +zone=56 +south +ellps=GRS80

Grid au_icsm_GDA94_GDA2020_conformal_and_distortion.tif needed but not 
found on the system. Can be obtained at 
https://cdn.proj.org/au_icsm_GDA94_GDA2020_conformal_and_distortion.tif


You'll notice all those 3 transformations are advertized to have a 3 
metre accuracy. The first one (note that the heuristics to determine the 
default one has changed among PROJ versions), which is the one used by 
default by cs2cs, uses the WGS84 = GDA94 hypothesis. The second and 
third ones us WGS84 = GDA2020

If you want to use the second one, you need to pass its pipeline to cct:

$ echo 0 0 0 | cct +proj=pipeline +step +inv +proj=tmerc +lat_0=-34 
+lon_0=150 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +step +proj=push +v_3 +step 
+proj=cart +ellps=WGS84 +step +inv +proj=helmert +x=0.06155 +y=-0.01087 
+z=-0.04019 +rx=-0.0394924 +ry=-0.0327221 +rz=-0.0328979 +s=-0.009994 
+convention=coordinate_frame +step +inv +proj=cart +ellps=GRS80 +step 
+proj=pop +v_3 +step +proj=utm +zone=56 +south +ellps=GRS80
   222908.2299   6233783.8497        0.0000           inf

Rule of thumb: as soon as one of your CRS is WGS84 related, you will 
generally go into troubles if you are interested in geodetic grade 
transformations.

Even

-- 
http://www.spatialys.com
My software is free, but my time generally not.
Mostly a fool, sometimes acting like a professional (or vice versa)
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20241011/10616cb5/attachment.htm>