[PROJ] 7.2.1 RC1 ob_tran changes break existing logic

[Kristian Evers]

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.

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>

/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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20201228/e560755c/attachment-0001.html>