[GIS] pyProj Coordinate Transformation Incorrect

coordinate systempyprojpython

I'm using pyProj to transform extents between two different geographical systems. In a case when both my systems are the same I expect to get back the same extents as inputted but this isn't happening for the below mentioned case:

import pyproj
inProj = pyproj.Proj("+init=EPSG:{0}".format(3857)) # Mercator
outProj = pyproj.Proj("+init=EPSG:{0}".format(3857))
x =  (-20037507.0672, 20037509.6184)
y =  (-1467048.29156, 8625918.8737)
print pyproj.transform(inProj,outProj,x,y)

The output comes out to be : ((-20037507.0672, -20037507.067178484), (-1467048.2915600014, 8625918.873699998))

Which is incorrect as my X-Coordinate extents are same for xMax and xMin.

If I change my co-ordinate systems

import pyproj
inProj = pyproj.Proj("+init=EPSG:{0}".format(3857)) # Mercator
outProj = pyproj.Proj("+init=EPSG:{0}".format(4326)) # WGS84
x =  (-20037507.0672, 20037509.6184)
y =  (-1467048.29156, 8625918.8737)
print pyproj.transform(inProj,outProj,x,y)

My X-coordinate output comes out to be : (-179.99998854118687, -179.99998854099357)

Which is again incorrect. Is this due to a domain issue in pyProj or is there any other reason ?

Best Answer

You make a mistake with the pyproj parameters 

from pyproj import Proj, transform
inProj  = Proj("+init=EPSG:3857"))
outProj = Proj("+init=EPSG:4326")
x1,y1 =  (-20037507.0672, 20037509.6184)
x2,y2 =  (-1467048.29156, 8625918.8737)
# with the x,y coordinates of pointx
print transform(inProj,outProj,x1,y1)
(-179.99998854118687, 85.08398750388278)
# with the x,y coordinates of point pointy
print transform(inProj,outProj,x2,y2)
(-13.178719028497799, 61.16289280460126)

Control with GDAL (Python GDAL/OGR Cookbook: Projections)

from osgeo import osr
from_srs = osr.SpatialReference()
from_srs.ImportFromEPSG(3857)
to_srs =  osr.SpatialReference()
to_srs.ImportFromEPSG(4326)
transf = osr.CoordinateTransformation(from_srs,to_srs)
print transf.TransformPoint(x1,y1)[:2]
(-179.99998854118684, 85.05112976833757)
# with the x,y coordinates of point y
print transf.TransformPoint(x2,y2)[:2]
(-13.178719028497795, 61.000416666720845)

If you want to use lists of coordinates:

x =  (-20037507.0672,-1467048.29156)
y =  (20037509.6184, 8625918.8737)
for i,j in zip(x,y):
   print transform(inProj,outProj,i, j)
(-179.99998854118687, 85.08398750388278)
(-13.178719028497799, 61.16289280460126)
Related Question