[PROJ] 7.2.1 RC1 ob_tran changes break existing logic

Roger Bivand Roger.Bivand at nhh.no
Mon Dec 28 03:41:20 PST 2020 

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

>
> /Kristian
>
>> On 28 Dec 2020, at 11:26, Kristian Evers <kristianevers at gmail.com> wrote:
>>
>> I can reproduce the pipeline you expect using the below call to projinfo, but it is unclear to me if your code is doing something similar behind the scenes. That is, using the projstring as a CRS definition in a transformation from the ob_tran “CRS” to a latlong CRS. Note that I’ve added “+type=crs” arguments to both projstrings to let projinfo know that these are legacy PROJ.4 CRS definitions.
>>
>>
>> $ PROJ_LIB=./data/ ./src/projinfo -s "+proj=ob_tran +o_proj=longlat +o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +type=crs" -t "proj=longlat +type=crs" -o PROJ
>> Candidate operations found: 1
>> -------------------------------------
>> Operation No. 1:
>>
>> unknown id, Inverse of unknown + Ballpark geographic offset from unknown to unknown, unknown accuracy, World, has ballpark transformation
>>
>> PROJ string:
>> +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
>>
>>
>> /Kristian
>>
>>> On 28 Dec 2020, at 11:10, Kristian Evers <kristianevers at gmail.com> wrote:
>>>
>>> Roger,
>>>
>>> Thanks for reporting this. We certainly need to find out what’s going on here. It is possible that this change shouldn’t have been included in 7.2.1.
>>>
>>> I can’t fully understand what you are trying to do from your code snippet below, can you provide an example that only uses PROJ tools or API?
>>>
>>> /Kristian
>>>
>>>> On 28 Dec 2020, at 10:18, Roger Bivand <Roger.Bivand at nhh.no> wrote:
>>>>
>>>> News item:
>>>>
>>>> createOperation(): make it work properly when one of the CRS is a BoundCRS of a DerivedGeographicCRS (+proj=ob_tran +o_proj=lonlat +towgs84=....) (#2441)
>>>>
>>>> appear to break:
>>>>
>>>> (crds <- matrix(data=c(9.05, 48.52), ncol=2))
>>>> (a <- project(crds, paste("+proj=ob_tran +o_proj=longlat",
>>>> "+o_lon_p=-162 +o_lat_p=39.25 +lon_0=180 +ellps=sphere +no_defs"),
>>>> use_ob_tran=TRUE, verbose=TRUE))
>>>> stopifnot(isTRUE(all.equal(a, matrix(c(-5.917698, -1.87195), ncol=2),
>>>> tolerance=.Machine$double.eps ^ 0.25)))
>>>>
>>>> completely in the R rgdal package. ob_tran has always needed special casing, because of the reverse invert direction, but this change is a major regression. The operation on input 9.05 48.52 should return -5.917698 -1.87195 but now chooses coordOp: +proj=noop +ellps=GRS80. In
>>>> 7.2.0, the chosen coordOp is (and probably should be):
>>>>
>>>> +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
>>>>
>>>> The NEWS item does not lead one to expect this regression.
>>>>
>>>> Hope some light can be thrown on this and that release will not be rushed.
>>>>
>>>> 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
>>>> _______________________________________________
>>>> PROJ mailing list
>>>> PROJ at lists.osgeo.org
>>>> https://lists.osgeo.org/mailman/listinfo/proj
>>>
>>> _______________________________________________
>>> PROJ mailing list
>>> PROJ at lists.osgeo.org
>>> https://lists.osgeo.org/mailman/listinfo/proj
>>
>> _______________________________________________
>> PROJ mailing list
>> PROJ at lists.osgeo.org
>> https://lists.osgeo.org/mailman/listinfo/proj
>
>

-- 
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

More information about the PROJ mailing list