[GIS] Incorrect transformation when using pyproj

coordinate systemprojpyprojpyshp

I am using pyproj to work with shp files and am relatively new to GIS. I am tryign to visualize the shp data, but am having a problem as I am incorrectly transforming the coordionates.

Here is an example from my code:

import shapefile
shp = shapefile.Reader(filePath)  # .shp file
rec = shp.shapeRecords()[0]
x, y = rec.shape.points[0]

Now, the resulting x and y are in utm format. For plotting, I want to convert to (lat,lon). I am using the following:

import pyproj
p = pyproj.Proj(init='EPSG:3857', proj='utm', zone=36N, ellps='WGS84')
p(x,y,inverse=True)

But this gives me incorrect (lat,lon) coordinates. the point in question is close to the Syrian-Israeli border, and has the raw value of (3970162.65625, 3853896.40625). But running p(x,y,inverse=True) gives the (lat,lon) coordinates of 67.809, 29.709, which turns out to be around Finland.

Any ideas why this transformation is wrong? I checked for the zone in here. And I am pretty sure the other parameters are correct as well. The shp file was created by exporting from ArcGIS desktop's ArcMap.

Edit in response to comment:

Here is the .prj file:

PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",
GEOGCS["GCS_WGS_1984",
DATUM["D_WGS_1984",
SPHEROID["WGS_1984",
6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Mercator_Auxiliary_Sphere"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
PARAMETER["Standard_Parallel_1",0.0],
PARAMETER["Auxiliary_Sphere_Type",0.0],
UNIT["Meter",1.0]]

Best Answer

For future users, here is some of my code. It might be helpful for a jumpstart.

from matplotlib.patches import Polygon
import pyproj
import shapefile

shp = shapefile.Reader(filePath)
recs = shp.shapeRecords()

p1 = pyproj.Proj(init='epsg:3857') 
p2 = pyproj.Proj(init='epsg:4326')

for j in range(1,len(recs)):
    rec = recs[j]
    verts = np.array(rec.shape.points)

    xin = np.array(verts)[:,0]
    yin = np.array(verts)[:,1]

    y,x = pyproj.transform(p1,p2,x=xin,y=yin)
    verts =  tuple(zip(y,x))
    poly = Polygon(verts, facecolor='red', edgecolor='0')
    ax.add_patch(poly)  

ax.autoscale_view()

plt.show()
Related Question