[GIS] Pyproj transform returning incorrect values (EPSG 3875 to EPSG 4326)

pyprojpython

I am trying to convert from EPSG 3875 to EPSG 4326 and have this code:

from pyproj import Proj, transform

def convert_coordinates(input_format, output_format, x, y):

    input_format = "+init=epsg:" +str(input_format)
    output_format = "+init=epsg:" +str(output_format)

    inProj = Proj(input_format)
    outProj = Proj(output_format)

    return transform(inProj,outProj,x,y)

x, y = 4152144, 219039
x2,y2 = convert_coordinates("3875", "4326", x, y)
print x2, y2

This returns 112.264404642 154.392043378 which is out of bounds for EPSG 4326. But ArcGIS returns 37.2993441706517, 1.96727415792103 which is in bounds and looks better.

Is there anything obviously wrong with this code?

Best Answer

The correct procedure is

from pyproj import Proj, transform
inProj  = Proj("+init=EPSG:3857"))
outProj = Proj("+init=EPSG:4326")
x, y = 4152144, 219039
print transform(inProj,outProj,x,y)
(37.29934417065166, 1.96727415792103)

Your error is convert_coordinates("3875", "4326", x, y): 3875 and 4326 are simple strings and not PyProj projections

Look at pyProj Coordinate Transformation Incorrect