2

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 ?

5

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)
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.