Regression in +proj=merc with shifted central meridian

[strezen]

Hi,

It seems that proj 5 has a regression when computing Mercator projection with non-default +lat_0 option. Lets try the forward transform of the equatorial points with left- and rightmost longitudes:

$ proj +proj=merc +lon_0=0 +lat_ts=0 +datum=WGS84
-180 0
-20037508.34	0.00
180 0
20037508.34	0.00

Seems to be Ok and it works the same way using the proj 4.9.3 and the latest one.
Now let's shift the central meridian with the +lat_0 option. proj 5 gives the following:

$ proj +proj=merc +lon_0=110 +lat_ts=0 +datum=WGS84
-70 0
-20037508.34	0.00
290 0
-20037508.34	0.00

which seems to be wrong. The proj 4.9.3 produces correct result:

$ proj +proj=merc +lon_0=110 +lat_ts=0 +datum=WGS84
-70 0
-20037508.34	0.00
290 0
20037508.34	0.00

It seems to be related to much more complicated logic of projection parameter preparation in pj_fwd.c but I am failing to find where it goes wrong.

Best regards,
Andrey

[kbevers]

I would argue that this was a bug in PROJ 4.9.3. -70 deg and 290 deg are equivalent, e.g. sin(-70 deg) == sin(290 deg) and cos(-70 deg) == cos(290 deg). So of course you get the same result for both. If you add the +over flag you get what you want:

$ proj
Rel. 5.0.1, April 1st, 2018
usage: proj [ -beEfiIlormsStTvVwW [args] ] [ +opts[=arg] ] [ files ]

$ ./bin/proj +proj=merc +lon_0=110 +lat_ts=0 +datum=WGS84 +over
290 0
20037508.34     0.00

[strezen]

This sounds reasonable, but the real problem here is in backward compatibility. It may kick in various places. I will adopt to new behavior though, so feel free to close the issue if you believe it is not a problem.

Thank you,
Andrey

[pelson]

Following on from the conversation in #985, and given the example above, wouldn't it be consistent for the following example to produce the same result without the +over flag?

$ proj +proj=merc +lon_0=0 +lat_ts=0 +datum=WGS84
-180 0
-20037508.34	0.00
180 0
20037508.34	0.00

This is indicating that the +over flag is necessary for lon_0<>0, but not for lon_0=0.

[kbevers]

This is indicating that the +over flag is necessary for lon_0<>0, but not for lon_0=0.

Really it is just indicating that a map has two edges and that you can control which on a coordinate end up at. As you can see from the output below it is nothing new. Adding +over doesn't change anything, since you don't step outside the normal range of -180..180.

(base) C:\dev>proj
Rel. 5.0.1, April 1st, 2018
usage: proj [ -beEfiIlormsStTvVwW [args] ] [ +opts[=arg] ] [ files ]

(base) C:\dev>proj +proj=merc +lon_0=0 +lat_ts=0 +datum=WGS84
180 12
20037508.34     1336830.21
-180 12
-20037508.34    1336830.21

(base) C:\dev>proj +proj=merc +lon_0=0 +lat_ts=0 +datum=WGS84 +over
180 12
20037508.34     1336830.21
-180 12
-20037508.34    1336830.21

(base) C:\dev>proj
Rel. 4.9.3, 15 August 2016
usage: proj.exe [ -bCeEfiIlormsStTvVwW [args] ] [ +opts[=arg] ] [ files ]

(base) C:\dev>proj +proj=merc +lon_0=0 +lat_ts=0 +datum=WGS84
pj_open_lib(proj_def.dat): call fopen(C:\Users\b012349\AppData\Local\Continuum\anaconda3\Library\share\proj_def.dat) - succeeded

180 12
20037508.34     1336830.21
-180 12
-20037508.34    1336830.21

(base) C:\dev>proj +proj=merc +lon_0=0 +lat_ts=0 +datum=WGS84 +over

180 12
20037508.34     1336830.21
-180 12
-20037508.34    1336830.21

[strezen]

The actual problem here is that proj 5 is not a drop-in replacement for proj 4. This change in behavior is breaking and should be listed somewhere as an incompatibility (or am I just missed it?).

[kbevers]

It isn't listed anywhere, no. This was a bug that was fixed as part of major refactoring where the projection logic was adapted to fit in a more generic transformation setting. As such this change was never planned, but a consequence of re-doing the logic and it was first discovered after release. I don't think we were aware that the original behaviour was ill-defined before this issue was created. As you may know, PROJ has historically been badly documented and has had a similarly bad test-suite. Many of these problems may have slipped by unnoticed. Hence why it isn't mentioned in the release notes. I can do that retro-actively if you think it is beneficial.

[strezen]

Please mention it somewhere (eqdc and merc are affected, maybe others). I am not arguing it was a bug which was fixed, that's all right, but taking in account a long history of PROJ a lot of software depends on it and its bugs. Any additional info could save some time to people involved.

Thanks!

[kbevers]

Please mention it somewhere

I'll do that. No one reads release notes so maybe there's a better place to put. Any suggestions?

(eqdc and merc are affected, maybe others)

More or less all projections, apart from a few odd ones that override longitude wrapping.

[rouault]

No one reads release notes so maybe there's a better place to put. Any suggestions?

Another possibility is to have a Backward incompatibility/behaviour changes page where we list for each release such changes

[kbevers]

Good idea. In what chapter of the docs do you suggest to put such a page?

[rouault]

Might make sense are a sub-chapter of https://proj4.org/download.html
Hum actually I see we do have https://proj4.org/development/migration.html but that one is mostly about the API. The new section about Backward incompatibility/behaviour changes could potentially point to https://proj4.org/development/migration.html for the 4.9.3 -> 5.0 transition

[kbevers]

Yeah, the migration page would likely be missed by non-API users of PROJ. Download section might work. Alternatively as a section under "Using PROJ" called "Differences between versions" or something like that?

[strezen]

a section under "Using PROJ" called "Differences between versions"

Looks good. And a link pointing to that section designed in the same way as "Attention!" blocks pointing to API changes. When the one opens https://proj4.org there are those API warnings and there should be additional one "See the 'Differences between versions' section".

[kbevers]

See #1012

[pelson]

Following on from my comment, I figured I must be confused with the description. Hopefully this analysis is helpful to others, but I compared the results of proj v4 & v5 with and without over for the two different rotations (i.e. with and without +lon_0=110):

In all cases, I'm calling proj v4 and v5 with exactly the same arguments, with v5 immediately following v4 calls:

$ echo -180 0 | proj -E +datum=WGS84 +proj=merc +lat_ts=0
v4: -180 0 -20037508.34 0.00
v5: -180 0 -20037508.34 0.00

$ echo 180 0 | proj -E +datum=WGS84 +proj=merc +lat_ts=0
v4: 180 0 20037508.34 0.00
v5: 180 0 20037508.34 0.00

Adding the +over flag does nothing, as it seems the projection is already handling the over in this case:

$ echo -180 0 | proj -E +datum=WGS84 +proj=merc +lat_ts=0 +over
v4: -180 0 -20037508.34 0.00
v5: -180 0 -20037508.34 0.00

$ echo 180 0 | proj -E +datum=WGS84 +proj=merc +lat_ts=0 +over
v4: 180 0 20037508.34 0.00
v5: 180 0 20037508.34 0.00

However, for the +lon_0=110 case:

$ echo -70 0 | proj -E +datum=WGS84 +proj=merc +lon_0=110 +lat_ts=0
v4: -70 0 -20037508.34 0.00
v5: -70 0 -20037508.34 0.00

$ echo 290 0 | proj -E +datum=WGS84 +proj=merc +lon_0=110 +lat_ts=0
v4: 290 0 20037508.34 0.00
v5: 290 0 -20037508.34 0.00

There is no magic "+over" functionality like there was for a standard +lon_0. Adding the +over flag does make the versions consistent though:

$ echo -70 0 | proj -E +datum=WGS84 +proj=merc +lon_0=110 +lat_ts=0 +over
v4: -70 0 -20037508.34 0.00
v5: -70 0 -20037508.34 0.00

$ echo 290 0 | proj -E +datum=WGS84 +proj=merc +lon_0=110 +lat_ts=0 +over
v4: 290 0 20037508.34 0.00
v5: 290 0 20037508.34 0.00

In summary, there seems to be special treatment when +lon_0 is not set for a longitude of +180 as stated in my comment. I'm completely OK with using +over - it is the right solution. I just wanted to point out the special case, and make it clear that anybody relying on this special behaviour should be aware that this is ripe for being fixed.

Use +over. Over and out. 😄