Sunday, 18 March 2018

coordinate system - Incorrect transformation when using pyproj


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]]

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()

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...