Inverse projection gives different results for pyproj 1.9.6 and 2.4.0

[suneeta-mall]

Code Sample, a copy-pastable example if possible

This python code is refered as ./pypro_bug.py in the descrption of issue.

#!/usr/bin/env python

from pyproj import Proj
EPSG_3857_PROJ4 = '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'

mercator_projection = Proj(EPSG_3857_PROJ4)
metres_x, metres_y = 3640945.346538976, -9059984.893837836
lon, lat = mercator_projection(metres_x, metres_y, inverse=True)
print(f"Projected long lat are:  ({lon}, {lat}) ")

Problem description

pyproj gives different results across two versions of pyproj (1.9.6 and 2.4.0).Please see this image and following commands for more details.

conda create --quiet --yes --channel conda-forge --name kidmenot python=3.6 pyproj=1.9.6
conda activate kidmenot
# content of ./pypro_bug.py is detailed above
./pypro_bug.py
 conda install --yes -q pyproj=2.4.0
./pypro_bug.py

Expected Output

I expected the projected long, lat to be same with pyproj 1.9.6 or 2.4.0.

Environment Information

• Output from: python -m pyproj -v
  pyproj info:
  pyproj: 2.4.0
  PROJ: 6.2.0
  data dir: /Users/suneeta.mall/miniconda3/envs/kidmenot/lib/python3.6/site-packages/pyproj/proj_dir/share/proj

System:
python: 3.6.7 | packaged by conda-forge | (default, Jul 2 2019, 02:07:37) [GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]
executable: /Users/suneeta.mall/miniconda3/envs/kidmenot/bin/python
machine: Darwin-19.0.0-x86_64-i386-64bit

Python deps:
pip: 19.3.1
setuptools: 41.4.0
Cython: None

• pyproj version (python -c "import pyproj; print(pyproj.__version__)")
• PROJ version (python -c "import pyproj; print(pyproj.proj_version_str)")
• PROJ data directory (python -c "import pyproj; print(pyproj.datadir.get_data_dir())")
• Python version (python -c "import sys; print(sys.version.replace('\n', ' '))")
• Operation System Information (python -c "import platform; print(platform.platform())")
  $python -c "import pyproj; print(pyproj.version)"
  1.9.6
  $python -c "import pyproj; print(pyproj.proj_version_str)"
  4.9.3
  $python -c "import pyproj; print(pyproj.datadir.get_data_dir())"
  /Users/suneeta.mall/miniconda3/envs/kidmenot/share/proj
  $python -c "import sys; print(sys.version.replace('\n', ' '))"
  3.6.7 | packaged by conda-forge | (default, Jul 2 2019, 02:07:37) [GCC 4.2.1 Compatible Clang 4.0.1 (tags/RELEASE_401/final)]
  $python -c "import platform; print(platform.platform())"
  Darwin-19.0.0-x86_64-i386-64bit

Installation method

• conda, pip wheel, from source, etc...
  conda

Conda environment information (if you installed with conda):

Environment (conda list):

$ conda list | grep -E "proj|aenum"

Details about conda and system ( conda info ):

$ conda info

[snowman2]

Hmmm, seems like something to look into. The output of proj and pyproj.Proj do differ.

$ echo 3640945.346538976 -9059984.893837836 | proj -d 8 -I +proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
32.70716853	-62.83532199
$ proj
Rel. 6.2.0, September 1st, 2019
usage: proj [-bdeEfiIlmorsStTvVwW [args]] [+opt[=arg] ...] [file ...]
$ python
Python 3.7.3 | packaged by conda-forge | (default, Jul  1 2019, 21:52:21) 
[GCC 7.3.0] :: Anaconda, Inc. on linux
Type "help", "copyright", "credits" or "license" for more information.
>>> from pyproj import Proj
>>> pj = Proj("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
>>> pj
Proj('+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +R=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs', preserve_units=True)
>>> pj(3640945.346538976, -9059984.893837836, inverse=True)
(32.70716853439809, -62.99092501209502)

[snowman2]

This gives the proper output:

>>> from pyproj import CRS, Transformer
>>> cc = CRS("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
>>> trans = Transformer.from_crs(cc, cc.geodetic_crs, always_xy=True)
>>> trans.transform(3640945.346538976, -9059984.893837836)
(32.70716853439809, -62.83532198512574)
>>> trans = Transformer.from_crs(cc.geodetic_crs, cc, always_xy=True)
>>> trans.transform(3640945.346538976, -9059984.893837836, direction="inverse")
(32.70716853439809, -62.83532198512574)

[snowman2]

@rouault has confirmed the current behavior is correct and gave a nice explanation here: OSGeo/PROJ#1699 (comment)

[suneeta-mall]

@snowman2 & @rouault Thanks for the detailed explanation.