[PROJ] 7.2.1 RC1 ob_tran changes break existing logic

[Roger Bivand]

Thanks!

On Mon, 28 Dec 2020, Kristian Evers wrote:

> Thanks, I now understand what you are doing. Basically you are doing 
> your best to emulate the legacy behaviour of PROJ.4. And that is the 
> root cause of this problem. At least partially. That is not saying that 
> it shouldn’t work, though. From the examples below it looks like 7.2.0 
> has two bugs that in conjunction returns the expected result. The two 
> bugs are
>
> 1. proj_create_crs_to_crs_from_pj() does not return a +proj=noop
>    transformation when two identical +proj=ob_tran CRS’s are used as
>    source and target.
> 2. proj_crs_get_geodetic_crs() returning a CRS that is identical to the
>    +proj=ob_tran CRS that it is fed. Should probably return
>    +proj=latlong +ellps=sphere instead.
>
> In 7.2.1 1) is fixed but not 2). I assume everything would work as 
> expected if 2) is also fixed. Unfortunately I am not familiar enough 
> with this part of the code base to provide a fix. Hopefully Even can 
> find a bit of time to look into that.

Yes, that seems a concise summary. And yes, the cases considered here were 
from legacy files, IIRC NetCDF weather/climate model output in ob_tran. If 
need be, I could force the proj_crs_get_geodetic_crs() output to be 
overridden in the > 7.2.0 case, but a proper fix would cover everybodies' 
needs.

>
> For reference, the project_ng_coordOp() function can be seen here: 
> https://github.com/cran/rgdal/blob/87d0d538a5f54681acb41de1b4647fc5a3e271f1/src/proj6.cpp#L656-L719 
> <https://github.com/cran/rgdal/blob/87d0d538a5f54681acb41de1b4647fc5a3e271f1/src/proj6.cpp#L656-L719>
>

The source of project_ng_coordOp() is at:

https://r-forge.r-project.org/scm/viewvc.php/pkg/src/proj6.cpp?root=rgdal&view=log

lines 767-855. The github version is an unauthorised copy and is pretty 
out of date (maintained on R-Forge/SVN because rgdal is maintenance-only, 
the sf package is the development route that we are trying to move users 
over to).

Best wishes,

Roger

> /Kristian
>
>
>> On 28 Dec 2020, at 12:41, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>>
>> On Mon, 28 Dec 2020, Kristian Evers wrote:
>>
>>> … and with cs2cs I get your expected result:
>>>
>>> $ echo 9.05 48.52 |PROJ_LIB=./data/ ./src/cs2cs -f %.6f "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs +type=crs" +to "proj=longlat +type=crs"
>>> -5.917698	-1.871950 0.000000
>>>
>>> I reckon there’s a bit of fine-tuning needed in your test suite. My best guess is adding “+type=crs” to the projstrings.
>>
>> Thanks for responding so quickly. The input to the rgdal C function project_ng_coordOp() is +proj=ob_tran +o_proj=longlat +o_lon_p=-162 \
>> +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs +type=crs, and inv=TRUE. In that function, the PROJ string is passed through proj_create(), then proj_crs_get_geodetic_crs() is used on that "target" to get the implied "source". If inv=TRUE, the coordOp is found by proj_create_crs_to_crs_from_pj() from "target" to "source". Why Proj4 inverted ob_tran is beyond me, but it has always done so. The proj=noop ellps=GRS80 transformation is also found without inversion.
>>
>> The values passed to proj_create_crs_to_crs_from_pj() with 7.2.1 are:
>>
>> from (proj_as_proj_string()): +proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs +type=crs
>>
>> to (proj_as_proj_string()): +proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs +type=crs
>>
>> output (proj_pj_info(pj_transform).definition): proj=noop ellps=GRS80
>>
>> output of proj_normalize_for_visualization() (proj_pj_info(pj_transform).definition): proj=noop ellps=GRS80
>>
>> With 7.2.0:
>>
>> from (proj_as_proj_string()): +proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs +type=crs
>>
>> to (proj_as_proj_string()): +proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs +type=crs
>>
>> output (proj_pj_info(pj_transform).definition): proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad step inv proj=ob_tran o_proj=longlat o_lon_p=-162 o_lat_p=39.25 lon_0=180 ellps=sphere step proj=unitconvert xy_in=rad xy_out=deg
>>
>> output of proj_normalize_for_visualization()
>> (proj_pj_info(pj_transform).definition): proj=pipeline step proj=unitconvert xy_in=deg xy_out=rad step inv proj=ob_tran o_proj=longlat o_lon_p=-162 o_lat_p=39.25 lon_0=180 ellps=sphere step proj=unitconvert xy_in=rad xy_out=deg
>>
>> I don't see what has happened. I can grasp the stupidity of from and to being equal PROJ strings in 7.2.0 but giving the valid pipeline; I can also  see that  proj_crs_get_geodetic_crs() might be puzzled by ob_tran and not return something like +proj=longlat +ellps=sphere.
>>
>> Any advice gratefully received.
>>
>> Best wishes,
>>
>> Roger
>
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; e-mail: Roger.Bivand at nhh.no
https://orcid.org/0000-0003-2392-6140
https://scholar.google.no/citations?user=AWeghB0AAAAJ&hl=en