Discrepancies in results produced by v4.9.3 vs v5.0.1

[avalentino]

The problem has been use cartopy (http://scitools.org.uk/cartopy) with proj v5

When I try to perform he same operation with two different versions of proj I get different results.

Example 1

$ proj 
Rel. 4.9.3, 15 August 2016
usage: proj [ -bCeEfiIlormsStTvVwW [args] ] [ +opts[=arg] ] [ files ]

$  echo -e "-170.0 -80.0 \n +190.0 +84.0" | proj +ellps=WGS84 +proj=merc +lon_0=10.0 +lat_ts=0.0 +units=m +no_defs
-20037508.34	-15496570.74 
20037508.34	18764656.23

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

$ echo -e "-170.0 -80.0 \n +190.0 +84.0" | proj +ellps=WGS84 +proj=merc +lon_0=10.0 +lat_ts=0.0 +units=m +no_defs
-20037508.34	-15496570.74 
-20037508.34	18764656.23

Example 2

$ cs2cs
Rel. 4.9.3, 15 August 2016
usage: cs2cs [ -eEfIlrstvwW [args] ] [ +opts[=arg] ]
                   [+to [+opts[=arg] [ files ]

$ echo -6.826286 54.725116 | \
cs2cs +ellps=WGS84 +a=57.29577951308232 +proj=eqc +lon_0=0.0 +no_defs \
+to +ellps=WGS84 +a=6377340.189 +b=6356034.447938534 +proj=tmerc +lon_0=-8 +lat_0=53.5 +k=1.000035 +x_0=200000 +y_0=250000 +units=m +no_defs
275614.87	386984.15 0.00

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

$ echo -6.826286 54.725116 | \
cs2cs +ellps=WGS84 +a=57.29577951308232 +proj=eqc +lon_0=0.0 +no_defs \
+to +ellps=WGS84 +a=6377340.189 +b=6356034.447938534 +proj=tmerc +lon_0=-8 +lat_0=53.5 +k=1.000035 +x_0=200000 +y_0=250000 +units=m +no_defs
275614.27	386984.21 0.00

[kbevers]

The first one is a duplicate of #941. Add +over.

I can reproduce the second example but I am kind of confused as to what is going on here. Why are you firt setting the ellipsoid to WGS84 and then specifying the semimajor axis (+x) and for the to system both semimajor and semiminor axes (+a/+b)? I suspect the difference is caused by changes in how the earth shape is determined from the various ellipsoid parameters. If you could elaborate a bit about what it is you are trying to achieve it is easier for me to determine if this really is a bug or not. And if it is a bug, if it is in 4.9.3 or 5.0.1.

[avalentino]

Examples above have been written to track down an issue in the cartopy (http://scitools.org.uk/cartopy) test suite.
More in detail example 2 corresponds to the test test_transverse_mercator.TestOSNI.test_default [1].
Parameters for the cs2cs program have been obtained by dumping (with the pj_get_def C function) the projPJ object (src_crs.proj4 and self.proj4) passed to the pj_transform function in [2].

So sorry I can't elaborate on the intention of the original author, just reporting parameters passed to the proj4 functions.

[1] https://github.com/SciTools/cartopy/blob/3c379ea50459a6b538529523c2bafe14c1c8dfc7/lib/cartopy/tests/crs/test_transverse_mercator.py#L97
[2] https://github.com/SciTools/cartopy/blob/3c379ea50459a6b538529523c2bafe14c1c8dfc7/lib/cartopy/_crs.pyx#L295

[kbevers]

I have done a bit of git bisecting and narrowed example 2 down to this commit

I am not sure what to make of that, really. It is hard to tell what is right or wrong here since the two proj-strings in the example both define the shape of the earth twice (+ellps=WGS84 combined with +a and/or +b). I need a more well-defined test case to make any progress on this. Clearly there's a change but it is not so clear if it is the 4.9.3 or 5.0.1 answer that is correct. Or both or none of them.

Unless this problem can be reproduced with a more sane definition of coordinate systems I am inclined to close this issue on the basis of too little information.

[avalentino]

Thanks @kbevers.
Maybe @pelson can give some hint about the cartopy initialization process that leads to the above Exemple 2 (as describe in my previous comment the starting point is test_transverse_mercator.TestOSNI.test_default).

[pelson]

I agree that

+ellps=WGS84 +a=57.29577951308232 +proj=eqc +lon_0=0.0 +no_defs \
+to +ellps=WGS84 +a=6377340.189 +b=6356034.447938534 +proj=tmerc +lon_0=-8 +lat_0=53.5 +k=1.000035 +x_0=200000 +y_0=250000 +units=m +no_defs

is somewhat junk. The problem is a fundamental one with cartopy's Globe definition, which has failed to implement mutually exclusive keyword arguments - that essentially means it is possible to construct such proj4 inits.

I'm relieved to note that the new behaviour (v5) can be achieved in v4 by removing the ellps definition. This gives some weight to the idea that the new behaviour (even with somewhat junk input) is doing the right thing.

In summary, I don't think there is anything that needs to be done on the proj.4 side of things, and this issue can be closed from my perspective. Thanks for raising this @avalentino, and thanks @kbevers for the effort you've put into answering this and other questions. I've been super impressed with the way proj4 is going - to the whole dev team: keep up the awesome work! 👍

[avalentino]

Thanks a lot @kbevers and @pelson