Inconsistency in 2D component between 2D and 3D transformation

Open ktheijs opened this issue 5 months ago • 1 comments

Example of problem

C:\OSGeo4W>echo 39.62996151 -77.77042607 | cs2cs EPSG:4326 EPSG:6488 -f "%.3f"
1095322.536     715954.909 0.000

C:\OSGeo4W>echo 39.62996151 -77.77042607 438.851 | cs2cs EPSG:4979 EPSG:6488+EPSG:6360 -f "%.3f"
1095323.285     715952.045 1556.137

Problem description

When doing a 3D transformation I expect the resulting 2D component to be the same as when I do a 2D transformation where I strip the vertical part from the transformation.

Expected Output

I expect that the X and Y value for the 3D transformation would also be 1095322.536 and 715954.909.

Environment Information

PROJ Rel. 9.6.2, June 4th, 2025
Windows 11

Installation method

I'm using command line on OSGeo4W GDAL 3.11.3 "Eganville", released 2025/07/12.

Sep 08 '25 10:09 ktheijs

I don't think this is a bug. The fundamental issue is that EPSG:4326 / EPSG:4979 is a datum ensemble of time-dependent WGS84 realizations, without a very precise relationship with NAD83(2011). So there are only approximate transformations available. Depending if you take the 2D or 3D path, the available transformations aren't the same. In both cases the result is geodeticaly questionable, but both answers are within the 2 m accuracy margin advertized. In the 2D case, the transformations available consider that WGS 84 ~= NAD83(2011) with a 2 m accuracy. In the 3D case, PROJ takes a lot of intermediate steps, so this is also questionable. I believe we must admit this is a limitation of the database that PROJ uses and of its algorithm.

$ projinfo EPSG:4326 EPSG:6488 --spatial-test intersects :
[...]
unknown id, Inverse of NAD83(2011) to WGS 84 (1) + SPCS83 Maryland zone (US survey foot), 2 m, Puerto Rico - onshore and offshore. United States (USA) onshore and offshore - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands - onshore and offshore.

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=lcc +lat_0=37.6666666666667 +lon_0=-77 +lat_1=39.45 +lat_2=38.3
        +x_0=399999.9998984 +y_0=0 +ellps=GRS80
  +step +proj=unitconvert +xy_in=m +xy_out=us-ft

$ projinfo EPSG:4979 EPSG:6488+6360 --spatial-test intersects
[...]
unknown id, Inverse of NAD83(HARN) to WGS 84 (3) + NAD83(HARN) to NAD83(FBN) (1) + NAD83(FBN) to NAVD88 height (1) + Inverse of Conversion from NAVD88 height (ftUS) to NAVD88 height + NAD83(FBN) to NAD83(2011) (NADCON5, CONUS) + SPCS83 Maryland zone (US survey foot), 1.14 m, United States (USA) - CONUS onshore - Alabama; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming., at least one grid missing

PROJ string:
+proj=pipeline
  +step +proj=axisswap +order=2,1
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +proj=cart +ellps=WGS84
  +step +inv +proj=helmert +x=-0.991 +y=1.9072 +z=0.5129 +rx=-0.0257899075194932
        +ry=-0.0096500989602704 +rz=-0.0116599432323421 +s=0
        +convention=coordinate_frame
  +step +inv +proj=cart +ellps=GRS80
  +step +proj=gridshift +grids=us_noaa_nadcon5_nad83_harn_nad83_fbn_conus.tif
  +step +inv +proj=vgridshift +grids=us_noaa_geoid03_conus.tif +multiplier=1
  +step +proj=unitconvert +z_in=m +z_out=us-ft
  +step +proj=gridshift +no_z_transform
        +grids=us_noaa_nadcon5_nad83_fbn_nad83_2007_conus.tif
  +step +proj=gridshift +no_z_transform
        +grids=us_noaa_nadcon5_nad83_2007_nad83_2011_conus.tif
  +step +proj=lcc +lat_0=37.6666666666667 +lon_0=-77 +lat_1=39.45 +lat_2=38.3
        +x_0=399999.9998984 +y_0=0 +ellps=GRS80
  +step +proj=unitconvert +xy_in=m +xy_out=us-ft

Sep 08 '25 11:09 rouault