Boundary oddities in Robinson projection

[QuLogic]

In Cartopy, the boundary for a (warped rectangular) map is calculated by making a rectangle spanning (-90, 90) latitude and (-180 + lon_0, 180 + lon_0) longitude, and then projecting it. For the Robinson projection, this used to work in 4.9.3, but produces strange results in 5+. This seems to be reproducible with just proj.

With 4.9.3:

$ echo -e '180 0\n-180 0' | proj +ellps=WGS84 +proj=robin +lon_0=0 +no_defs
17005833.33	-0.00
-17005833.33	-0.00
$ echo -e '182 0\n-178 0' | proj +ellps=WGS84 +proj=robin +lon_0=2 +no_defs
17005833.33	-0.00
-17005833.33	-0.00
$ echo -e '178 0\n-182 0' | proj +ellps=WGS84 +proj=robin +lon_0=-2 +no_defs
17005833.33	-0.00
-17005833.33	-0.00

With 5.0.0:

$ echo -e '180 0\n-180 0' | ./src/proj +ellps=WGS84 +proj=robin +lon_0=0 +no_defs
17005833.33	-0.00
-17005833.33	-0.00
$ echo -e '182 0\n-178 0' | ./src/proj +ellps=WGS84 +proj=robin +lon_0=2 +no_defs
-17005833.33	-0.00
-17005833.33	-0.00
$ echo -e '178 0\n-182 0' | ./src/proj +ellps=WGS84 +proj=robin +lon_0=-2 +no_defs
17005833.33	-0.00
17005833.33	-0.00

Here, if central longitude is positive, the positive boundary wraps around, and if the central longitude is negative, the negative boundary wraps around. The same occurs on master.

[QuLogic]

Bisect tells me that 90968ff is the first bad commit.

commit 90968ff934a348b02f982080f8b7ccf17e037a46
Author: Thomas Knudsen <thokn@sdfe.dk>
Date:   Fri Jan 5 12:03:21 2018 +0100

    Introduce compatibility for cs2cs-style proj-strings into the 4D API.
    
    Parameters such as towgs84, nadgrids and geoidgrids was previously only
    handled by pj_transform(). This commit add a compatibility layer in
    proj_create() by calling the pj_cs2cs_emulation_setup() function. This
    function sets up a handful of predefined transformation objects on the
    PJ object that is being created. Each of these transformation objects
    are related to the cs2cs-style parameters we are trying to emulate in
    the 4D API. That is, if the +towgs84 parameters is used we create
    P->helmert with the parameters specified in +towgs84. Similarly for
    +axis, +nadgrids and +geoidgrids. When these transformation objects
    exists we use them in the prepare and finalize functions in pj_fwd/
    pj_inv. If no cs2cs-style parametes are specified we skip those
    parts of the prepare and finalizing steps.
    
    Co-authored-by:Thomas Knudsen <thokn@sdfe.dk>
    Co-authored-by:Kristian Evers <kristianevers@gmail.com>

:040000 040000 c9a447b89baf0191b49885c976a79df2d899b60f 9f7bb44eeffa820616a4a739b31d014fc0341965 M	src
:040000 040000 036bea3f6d5f9d1874c77956c16734ef6fd43d5d c0c2cd8531941d248c056b90fb98c32c72fa86ea M	test

[QuLogic]

Possibly the same thing happening in Mercator?

# 4.9.3
$ echo -e '190 0\n-170 0' | proj +ellps=WGS84 +proj=merc +lon_0=10 +x_0=0.0 +y_0=0.0 +units=m +no_defs
20037508.34	0.00
-20037508.34	0.00

# 90968ff
$ echo -e '190 0\n-170 0' | ./src/proj +ellps=WGS84 +proj=merc +lon_0=10 +x_0=0.0 +y_0=0.0 +units=m +no_defs
-20037508.34	0.00
-20037508.34	0.00

[rouault]

This is a duplicate of #941
This is documented at https://proj4.org/usage/differences.html#longitude-wrapping-when-using-custom-central-meridian

[QuLogic]

Bleh, not only did I already know about +over, but I've also read that ticket before. Guess that's what happens when you leave things alone for four months. I don't quite understand why that change only seems to affect Mercator and Robinson though.

In any case, adding +over to those two fixes most of the Cartopy tests. There's still one broken Robinson plot that I need to investigate a bit still.

[QuLogic]

OK, the trouble with +over is that it seems to disable wrapping, and anything past 180 disappears. It just wasn't obvious in other tests for some reason. For example, here is what Cartopy produces with central longitude of 135 (+ellps=WGS84 +proj=robin +lon_0=135 +over +no_defs)

I tried also adding +lon_wrap=135, but that had no effect.

I'm still confused how the boundary worked, but none of the interior content, so I will continue to investigate. I will note that the documentation plot scripts also do the same thing, but I'm not sure if they were written to understand non-default central longitude.

[QuLogic]

We worked around this change by subtracting a small epsilon (1e-10) from the opposite boundary. This works for both old and new versions.