Incorrect propagation of +R parameter in transform pipeline

[scottmcnab]

Example of problem

There is an inconsistency in behaviour between PROJ4 and PROJ9 for the transform from EPSG:4978 (ECEF) to ESRI:104199.

Here are the definitions of the two CRS:

$ projinfo epsg:4978 -o PROJ
PROJ.4 string:
+proj=geocent +datum=WGS84 +units=m +no_defs +type=crs

$ projinfo esri:104199 -o PROJ
Warning: object is deprecated

PROJ.4 string:
+proj=longlat +R=6378137 +no_defs +type=crs

When transforming a point using PROJ4 (version 4.9.3), I get the expected result:

$ echo "4057449.87922115 1236185.10760693 4747287.50008112" | cs2cs +proj=geocent +datum=WGS84 +units=m +no_defs +type=crs +to +proj=longlat +R=6378137 +no_defs +type=crs -f "%.10f"
16.9444298151   48.4111246681 -74.6128296684

However, when using PROJ9 (version 9.1.0), I get a different (incorrect) result:

$ echo "4057449.87922115 1236185.10760693 4747287.50008112" | cs2cs +proj=geocent +datum=WGS84 +units=m +no_defs +type=crs +to +proj=longlat +R=6378137 +no_defs +type=crs -f "%.10f"
16.9444298151   48.2199854103 -11992.8176014163

Even using Authority Codes with PROJ9 results in the same incorrect result:

$ echo "4057449.87922115 1236185.10760693 4747287.50008112" | cs2cs epsg:4978 esri:104199 -f "%.10f"
48.2199854103   16.9444298151 -11992.8176014163

Problem description

Inspecting the PROJ9 transform pipeline shows:

$ projinfo -s epsg:4978 -t esri:104199 -o PROJ
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Ballpark geocentric translation from WGS 84 to GCS_WGS_1984_Major_Auxiliary_Sphere (geocentric) + Conversion from GCS_WGS_1984_Major_Auxiliary_Sphere (geocentric) to GCS_WGS_1984_Major_Auxiliary_Sphere, unknown accuracy, World, has ballpark transformation

PROJ string:
+proj=pipeline
  +step +inv +proj=cart +R=6378137
  +step +proj=longlat +R=6378137
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

I think the problem is the first step incorrectly contains an unexpected +R parameter:
+step +inv +proj=cart +R=6378137

This transform step should use the ellipsoid definition from EPSG:4978 (WGS84), and not include the +R value, which should only apply to the ESRI:104199 transform step. I expected this first step to be something like:
+step +inv +proj=cart +ellps=WGS84

Expected Output

I would expect to see a transform pipeline similar to:

+proj=pipeline
  +step +inv +proj=cart +ellps=WGS94
  +step +proj=longlat +R=6378137
  +step +proj=unitconvert +xy_in=rad +xy_out=deg
  +step +proj=axisswap +order=2,1

In fact, if I manually construct a transform with this first step corrected (and remove the axisswap), then I get an output that exactly matches the expected output from PROJ4:

$ echo "4057449.87922115 1236185.10760693 4747287.50008112" | cct -d 10 +proj=pipeline +step +inv +proj=cart +ellps=WGS84 +step +proj=longlat +R=6378137 +step +proj=unitconvert +xy_in=rad +xy_out=deg
 16.9444298151   48.4111246681  -74.6128296675           inf

It looks to me that there is some logic in the pipeline builder that is incorrectly propagating the +R parameter between steps, when it should not.

Is this behaviour expected? Or is there some other way to get behaviour more consistent with PROJ4?

Many thanks

Environment Information

• PROJ version (9.1.0)
• Operation System: Ubuntu Linux

Installation method

• Compiled from source

[rouault]

I've studied this and my conclusion is that PROJ 4.9 or PROJ9 are both correct & incorrect at the same time, or said otherwise provide arbitrary results. This is just a ill-defined problem given there is no known datum transformation between WGS 84 (used by EPSG:4978) and D_WGS_1984_Major_Auxiliary_Sphere.
The difference come from the fact that PROJ 4.9 and PROJ9 use totally different logic to compute coordinate operations.
In PROJ 4.9, a null-datum transformation was done in the geographic coordinate space, that is the input cartesian coordinates were converted into geographic ones (in WGS 84 datum), then those geographic coordinates were considered as the same in the target datum (WGS_1984_Major_Auxiliary_Sphere), and then converted to the output CRS (this step being a no-op here)
In current PROJ versions, the null-datum transformation in the case of geocentric <--> geographic of different datums is done in the geocentric coordinate space.
That said in #3361 I propose to get back to PROJ < 6 era behavior.