Option +over broken in PROJ 7.2

[mheiskan]

Example of problem

The following command line converts between normal and rotated longlat coordinates, and works fine with PROJ 4.8, the result is the configured location of the south pole:

> echo 0 -90 | cs2cs -f "%.7f" +proj=ob_tran +o_proj=longlat +o_lon_p=0 +lon_0=20 +o_lat_p=30 +to_meter=.0174532925199433 +R=6371220 +wktext +over +towgs84=0,0,0 +no_defs +to +proj=longlat +R=6371220 +towgs84=0,0,0 +no_defs
20.0000000  -30.0000000 0.0000000

However, with PROJ 7.2 one gets the expected result only by removing the +over option:

> echo 0 -90 | /usr/proj72/bin/cs2cs -f "%.7f" +proj=ob_tran +o_proj=longlat +o_lon_p=0 +lon_0=20 +o_lat_p=30 +to_meter=.0174532925199433 +R=6371220 +wktext +over +towgs84=0,0,0 +no_defs +to +proj=longlat +R=6371220 +towgs84=0,0,0 +no_defs
0.3490659   -0.9882081 0.0000007
> echo 0 -90 | /usr/proj72/bin/cs2cs -f "%.7f" +proj=ob_tran +o_proj=longlat +o_lon_p=0 +lon_0=20 +o_lat_p=30 +to_meter=.0174532925199433 +R=6371220 +wktext +towgs84=0,0,0 +no_defs +to +proj=longlat +R=6371220 +towgs84=0,0,0 +no_defs
20.0000000  -30.0000000 0.0000000

The unexpected result seems to have the longitude in radians, but the latitude is broken in some more complicated way.

Environment Information

• PROJ version 7.2
• RHEL 7.9
• RPM from pgdg-common

[kbevers]

Quite likely this is related to #941. If so, it is not a bug, rather a bug that was fixed.

Why do you need +over? You seem to be getting what you want without it.

[mheiskan]

This was a test simplified to a single coordinate. Sure, I will get the correct result I want when I remove +over for this test, but I have not gotten into our bigger tests yet which include polygon and raster data in longitude range 0...360, which is common for global weather models. Even if our remaining tests work too, this could still break somebody else's code. +over should not change the expected result of 20 degrees to radians, and do something more complicated for the latitude. As such I do not quite see how this could be related to ticket #941 you mentioned.

[kbevers]

As such I do not quite see how this could be related to ticket #941 you mentioned.

I mentioned that ticket since the +over behaviour changed from 4 to 5. You are making quiting a jump from 4.8.0 to 7.2, this was first on the list of likely issues. I agree that it probably isn't what's wrong here. As you may be aware of, PROJ 7 is an entirely different beast than PROJ.4 and one of the big differences is the way transformations are instantiated. The PROJ.4 way of using projstrings as CRS definitions is now deprecated and should be avoided where possible. Use EPSG-codes or WKT2 descriptions instead if at all possible.

By "modernizing" your projstrings a bit I can at least get around the radian issue:

echo 0 -90 | cs2cs -f "%.7f" "+proj=ob_tran +o_proj=longlat +o_lon_p=0 +lon_0=20 +o_lat_p=30 +to_meter=.0174532925199433 +R=6371220 +over +type=crs" +to "+proj=longlat +R=6371220 +type=crs"
20.0000000      -56.6201562 0.0000000

In the above I've removed +no_def, +towgs84=0,0,0 and +wktext. +wktext has never been used by PROJ (it's a GDAL thing), +no_def was removed in version 6 (I think it was) and the +towgs84=0,0,0 amounts to nothing since it shifts the coordinates by 0.0 meters. The +to_meters could be removed as well in this particular case. I've added +type=crs to make sure that PROJ treats the projstrings as CRS definitions.

To explain what's wrong with the latitude we can investigate the transformation with projinfo:

projinfo -s "+proj=ob_tran +o_proj=longlat +o_lon_p=0 +lon_0=20 +o_lat_p=30 +R=6371220 +over +type=crs" -t "+proj=longlat +R=6371220 +type=crs" -o PROJ
Candidate operations found: 1
-------------------------------------
Operation No. 1:

unknown id, Transformation from unknown to unknown, unknown accuracy, unknown domain of validity

PROJ string:
+proj=pipeline
  +step +inv +proj=ob_tran +o_proj=longlat +o_lon_p=0 +lon_0=20 +o_lat_p=30
        +R=6371220 +over
  +step +proj=longlat +R=6371220
  +step +proj=unitconvert +xy_in=rad +xy_out=deg

From the resulting pipeline we can see that the input coordinate is not transformed to radians before being passed to ob_tran but the finishing step converts from radians to degrees. This is definitely a bug. The resulting pipeline should be

+proj=pipeline
  +step +proj=unitconvert +xy_in=deg +xy_out=rad
  +step +inv +proj=ob_tran +o_proj=longlat +o_lon_p=0 +lon_0=20 +o_lat_p=30
        +R=6371220 +over
  +step +proj=longlat +R=6371220
  +step +proj=unitconvert +xy_in=rad +xy_out=deg

The problem here is that ob_tran can digest both projected and angular coordinates and the logic doesn't account for that it would seem. Internally PROJ uses radians to handle angular coordinates so passing degrees instead of radians will obviously produce the wrong results.

Depending on your use case and surrounding setup we may be able to come up with a work around.