MGA94 and MGA2020 unexpected transformation

[jamiecook]

This is most likely a user error - but we've had two experienced GIS users scratching their heads on this one. We appear to be getting a bad transformation out of the MGA projections of Australia.

Code Sample (copy-pastable)

from pyproj import Proj, Transformer

wgs84,wgs84p,gda94,gda2020 = 'epsg:4326', 'epsg:3857', 'epsg:4348', 'epsg:7842'
mga2020 = 7856

def round_trip(point, epsg):
    p1,p2 = Proj(wgs84), Proj(epsg)
    t_forward = Transformer.from_proj(proj_from=p1, proj_to=p2, always_xy=True)
    t_back    = Transformer.from_proj(proj_from=p2, proj_to=p1, always_xy=True)
    point_t = t_forward.transform(*point)

    print(f"round trip from {p1.name} -> {p2.name} -> {p1.name}")
    print(f"{wgs84}   => {point}")
    print(f"{epsg}   => {point_t}")
    print(f"{wgs84}   => {t_back.transform(*point_t)}")
    print()

point = (93.41, -60.55)
point = (153, -27) # near Brisbane, Australia
round_trip(point, wgs84p)
round_trip(point, gda94)
round_trip(point, mga2020)

Problem description

So taking a point (which is in the defined bounds of the cartesian projection system) and going from WGS84, into that projection system and back out again should give (roughly) the same point.

This property holds for projected WGS84 (3857) and for Map Grid of Australia (MGA2020:7856), but does not hold for GDA projections (GDA94 'epsg:4348', GDA2020 'epsg:7842') both of which are geo-centric.

Expected Output

round trip from longlat -> merc -> longlat
epsg:4326   => (153, -27)
epsg:3857   => (17031882.091370855, -3123471.749104577)
epsg:4326   => (152.99999999999997, -27.0) # THIS IS EXPECTED

round trip from longlat -> geocent -> longlat
epsg:4326   => (153, -27)
epsg:4348   => (-5067052.800615102, 2581792.355845415)
epsg:4326   => (153.0, 0.0) # THIS IS NOT EXPECTED

round trip from longlat -> utm -> longlat
epsg:4326   => (153, -27)
7856   => (500000.0, 7013564.757530207)
epsg:4326   => (153.0, -26.999999999999993) # THIS IS EXPECTED

Environment Information

λ pyproj -v
pyproj info:
    pyproj: 3.6.1
      PROJ: 9.3.0
  data dir: /home/jamie/git/tmr/qfm/venv/lib/python3.10/site-packages/pyproj/proj_dir/share/proj
user_data_dir: /home/jamie/.local/share/proj
PROJ DATA (recommended version): 1.15
PROJ Database: 1.2
EPSG Database: v10.094 [2023-08-08]
ESRI Database: ArcGIS Pro 3.1 [2023-19-01]
IGNF Database: 3.1.0 [2019-05-24]

System:
    python: 3.10.12 (main, Jul 29 2024, 16:56:48) [GCC 11.4.0]
executable: /home/jamie/git/tmr/qfm/venv/bin/python
   machine: Linux-5.15.153.1-microsoft-standard-WSL2-x86_64-with-glibc2.35

Python deps:
   certifi: 2023.7.22
    Cython: 3.0.5
setuptools: 68.2.2
       pip: 23.3.1

Installed with pip

[jjimenezshaw]

Geocentric systems need 3D coordinates. Try with (153, -27, 0) instead of (153, -27). That will produce a meaningful 3D coordinate that makes sense as input for the geocentric system.

[jamiecook]

Right - i guess that makes sense.

I guess it's outside scope for this forum - but any suggestions on how that would be applied to a GeoDataFrame where all my geometries are 2D?

[jjimenezshaw]

Why do you want to use a geocentric system if your data is 2D? A geographic CRS is not enough?

[jamiecook]

I was looking for a cartesian system that had good coverage of the entire Australian continent - MGA2020 works okay within each zone - but i have an "all of Aus" dataset that I wanted to do some buffering operations with.

Now that you highlight it, there were probably other better options, but I was nerd-sniped by the unexpected behaviour here.

[jjimenezshaw]

a cartesian system that had good coverage of the entire Australian continent

Australia is not Austria (sorry, I had to make the joke). I mean, it is huge! A Cartesian 2D would be projected... with a lot of distortion due to the size of Australia. Ok, Geocentric is also Cartesian, but you are opening another Pandora's box: it is a 3D system, so all your coordinate management should support that.

[snowman2]

Possibly related: OSGeo/PROJ#2318

May be worth trying with a 3D geographic CRS (EPSG:4979).

[m-richards]

@jjimenezshaw is there a reason why pyproj will let you transform 2d data via a geocentric crs if 3d inputs are required?

After some experimenting I understand in this case the issue is because (153, -27) is converted to (-5067052.801, 2581792.356), when it really should be (-5067052.801, 2581792.356, -2878215.576) in the first transformation, and it is the missing z coordinate (which is treated as zero) responsible for this failing to round trip in the second transformation.

If there's a reason for pyproj to support this, then perhaps there's a case to prevent/ warn on these transformations e.g. at a geopandas level in this case - if there's an easy way to classify a crs / transformation as non-applicable from 2d to 3d (I imagine geocentric transformations might not be the only kind in this situation?)