Unitconvert from degrees to radians not working with cct

[trey-stafford]

Example of problem

The latest version of proj does not appear to correctly perform unit conversion from degrees to radians.

In the following example, I use proj v9.5.0 cct to construct a pipeline that should convert from degrees to radians. The output appears to be in degrees.

$ docker run -it osgeo/proj:9.5.0 bash
root@e260fb4c9363:/# echo "-117.61748591770593, 59.4967490316935, 329.39612, 1990.5" | cct -v -d 14 +proj=pipeline +ellps=WGS84 +step +proj=unitconvert +xy_in=deg +xy_out=rad
-117.61748591770591  59.49674903169350  329.39612000000000     1990.5000

The transformation works as expected using proj v5.2:

$ docker run -it osgeo/proj:5.2 bash
root@2c3e5a11ed10:/# echo "-117.61748591770593, 59.4967490316935, 329.39612, 1990.5" | cct -v -d 14 +proj=pipeline +ellps=WGS84 +step +proj=unitconvert +xy_in=deg +xy_out=rad
-2.05281238718203  1.03841416483580  329.39612000000000     1990.5000

Problem description

I expect the output of unitconvert +xy_in=deg +xy_out=rad to be in radians. Instead the output is in degrees.

Expected Output

I expect the output to be in radians.

Environment Information

• PROJ version (proj): v9.5.0
• Operation System Information: Ubuntu/Docker

Installation method

• Docker. I have also tested this against conda-installed proj/pyproj.

[rouault]

I believe this is a side effect of a change of many years ago of edb0684#diff-6ddde4f3f187062a5e8b59f634a5c067d4585ce0d9779227e215725515538bcbR482 where the output of +proj=unitconvert +xy_in=deg +xy_out=rad is rightfully qualified as PJ_IO_UNITS_RADIANS instead of PJ_IO_UNITS_WHATEVER. But cct considers since its very beginnings that when proj_angular_output() == true (which means the output unit is radians), then the result is output as degree (the rational is that for example if the last step is for example something like "+inv +step +proj=utm ....", which outputs radians, you generally want to see the result as degree)

@kbevers any thoughs on what we should do here:

• do nothing but document that oddity ?
• add a hack in cct to detect if there's an explicit last step to convert to radians, and give up the post conversion to degree ?

[trey-stafford]

Hi @rouault , thanks for your timely response!

But cct considers since its very beginnings that when proj_angular_output() == true (which means the output unit is radians), then the result is output as degree (the rational is that for example if the last step is for example something like "+inv +step +proj=utm ....", which outputs radians, you generally want to see the result as degree)

This makes sense to me - I was wondering if something like this might be happening under the hood. Thanks for confirming! At a minimum, I think a note about this in the docs would be great!

[kbevers]

I'm surprised this hasn't surfaced earlier! I suppose it's because internally we more or less always couple unitconvert with an axisswap in the previous step. If it had been done the other way round we would have noticed this right away.

add a hack in cct to detect if there's an explicit last step to convert to radians, and give up the post conversion to degree ?

My immediate thought is that cct it should do exactly what it is asked to do. I can accept an exception when cct is only fed a projection, as I think the output would be too surprising and cause more confusion than is necessary. If I'm not mistaken there's a similar problem with the pre-conversion to radians. This gie test highlights a few problems:

<gie>

operation +proj=unitconvert +xy_in=rad +xy_out=deg

accept 1 1 1 1
expect 57.2957795131   57.2957795131        1.0000        1.0000
roundtrip 1

operation +proj=unitconvert +xy_in=deg +xy_out=rad

accept 57.2957795131   57.2957795131        1.0000        1.0000
expect 1 1 1 1
roundtrip 1

operation +proj=pipeline 
          +step +proj=unitconvert +x_in=rad +xy_out=deg 
          +step +proj=noop 
          +step +proj=unitconvert +xy_in=deg +xy_out=rad

accept 1 1 1 1
expect 1 1 1 1
roundtrip 1

</gie>

results:

gie unitconvert.gie
-------------------------------------------------------------------------------
Reading file 'unitconvert.gie'
     -----
     FAILURE in unitconvert.gie(6):
     expected: 57.2957795131 57.2957795131 1.0000 1.0000
     got:      1.000000000000   1.000000000000   1.000000000   1.000000000
     deviation:  79614.254892 mm,  expected:  0.500000 mm
     -----
     FAILURE in unitconvert.gie(12):
     expected: 1 1 1 1
     got:      57.295779513100   57.295779513100   1.000000000   1.000000000
     deviation:  999999999.999000 mm,  expected:  0.500000 mm
     -----
     FAILURE in unitconvert.gie(21):
     expected: 1 1 1 1
     got:      57.295779513082   57.295779513082   1.000000000   1.000000000
     deviation:  999999999.999000 mm,  expected:  0.500000 mm
-------------------------------------------------------------------------------
total:  3 tests succeeded,  0 tests skipped,  3 tests FAILED!
-------------------------------------------------------------------------------

The tests that passes are the inverse operation started with the roundtrip instruction, which isn't too surprising.

If we can make the above test pass without breaking any other tests I think we're good.

[rouault]

This gie test highlights a few problems:

I'm not sure how gie is relevant here. The post-unit conversion from radians to degree is a cct specific code path

[kbevers]

I'm not sure how gie is relevant here. The post-unit conversion from radians to degree is a cct specific code path

Maybe I've misunderstood the problem, but to my eyes it looks like it's the same issue in gie, so perhaps it can be solved at a higher level than in the individual apps? In any case, I didn't see any connection to cct (or gie) in the changeset you linked, so didn't think that it was specific to either of those applications

[rouault]

I didn't see any connection to cct

yes sorry my analysis was a bit terse. This is related to

PROJ/src/apps/cct.cpp

Line 527 in 01194d8

if (proj_angular_output(P, direction) ||

which has been there since the beginnings of cct. But the fact that we reports the radians unit as the output unit of +proj=unitconvert +xy_in=deg +xy_out=rad is a later change.

gie has a similar logic in

PROJ/src/apps/gie.cpp

Line 1052 in 01194d8

T.b = proj_angular_output(T.P, T.dir) ? todeg_coord(T.P, T.dir, co) : co;

to convert from radians to degrees

[zyizhang]

I believe this is a side effect of a change of many years ago of edb0684#diff-6ddde4f3f187062a5e8b59f634a5c067d4585ce0d9779227e215725515538bcbR482 where the output of +proj=unitconvert +xy_in=deg +xy_out=rad is rightfully qualified as PJ_IO_UNITS_RADIANS instead of PJ_IO_UNITS_WHATEVER. But cct considers since its very beginnings that when proj_angular_output() == true (which means the output unit is radians), then the result is output as degree (the rational is that for example if the last step is for example something like "+inv +step +proj=utm ....", which outputs radians, you generally want to see the result as degree)

@kbevers any thoughs on what we should do here:

• do nothing but document that oddity ?
• add a hack in cct to detect if there's an explicit last step to convert to radians, and give up the post conversion to degree ?

Thanks for clarification. I also met this issues. And in some cases, I am wondring if this action will destroy the whole pipeline? For example, after ' +step +inv +proj=utm', we want to do ' +step +proj=hgridshift +grids=ch_swisstopo_CHENyx06a.tif ' , for which the inputs should be in rad but not degree (or always taking degrees as inputs?)? I am not sure about the solutions.
Thanks in advance.