3

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]]
Improve this question
10
  • 1
    EPSG 3857 isn't a UTM projection, Zone 36 north is EPSG:32636. The EPSG code should be all you need... proj, zone, ellips are all specified in the EPSG. Those numbers don't look right for UTM, are you sure it's UTM and not web mercator or Lamberts?
    – Michael Stimson
    Commented Apr 13, 2017 at 1:02
  • Thanks for pointing that out. I am not sure, is there perhaps a way to check this?
    – splinter
    Commented Apr 13, 2017 at 1:06
  • 2
    Where do the points come from? You say shapefile so there should be a .prj file with it, open in Notepad (or your favorite text editor) and copy/paste the whole string into your question. If there's not a prj file then you need to find out from the agency or people you got the data from what spatial reference they used to store the data in.
    – Michael Stimson
    Commented Apr 13, 2017 at 1:07
  • 1
    there is an example jswhit.github.io/pyproj in your case p1 = Proj(init='epsg:3857') p2 = Proj(init='epsg:4326') to project from Web Mercator to WGS84 geographic.
    – Michael Stimson
    Commented Apr 13, 2017 at 1:21
  • 1
    That would be another question. The rules here are one question per question, if you have another then please ask on a new thread. When you say 'cartesian' what do you mean? Mathematically cartesian coordinates are X, Y (optional Z) which can describe just about any spatial reference.
    – Michael Stimson
    Commented Apr 13, 2017 at 1:59

1 Answer 1

Reset to default
3

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()
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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