coordinates conversion with pyproj

5
votes

I am not experienced in gis coordinates conversion but managed, using this page: http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/ to convert shapefile coordinates from EPSG:28992 to EPSG:4326 using the python module pyproj using these statements: 

wgs84=pyproj.Proj("+init=EPSG:4326")
epsg28992=pyproj.Proj("+init=EPSG:28992")
pyproj.transform(epsg28992, wgs84,x,y)

When I reverse and enter these coordinates in google maps they give me correct locations. So this is working fine.

Now I have another shapefile(s) and I look at the shapefile.prj file to determine what projection was used. The ESRI WKT corresponds with ESRI:102686 which I find here: http://epsg.io/102686 As the ESRI:102686 code is not known by pyproj (gives error), I have to use proj4 notation which I got from the same site (http://epsg.io/102686):

wgs84=pyproj.Proj("+init=EPSG:4326") 
esri102686=pyproj.Proj("+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333 +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000.0000000001 +datum=NAD83 +units=us-ft +no_defs")
pyproj.transform(esri102686, wgs84,x,y)

I get e.g. coordinates and use these in google maps: 60.275122729462495, -61.873986125999316 which is somewhere in the ocean...

But my results should be in Cambridge, MA in the US, so more around: 41.00000, -71,5000000

What am I doing wrong?

2
pythongisshapefile

2 Answers

5
votes

Solved, added preserve_units = True, like this:

esri102686 = pyproj.Proj("+proj=lcc +lat_1=41.71666666666667 +lat_2=42.68333333333333                +lat_0=41 +lon_0=-71.5 +x_0=200000 +y_0=750000.0000000001 +datum=NAD83 +units=us-ft     +no_defs",preserve_units= True)

Now it works fine. If the optional keyword 'preserve_units' is True, the units in map projection coordinates are not forced to be meters. See here.

0
votes

Recommended method with pyproj 2.x:

>>> import pyproj
>>> pyproj.__version__
'2.2.0'
>>> from pyproj import Transformer
>>> transformer = Transformer.from_crs("ESRI:102686", "EPSG:4326")
>>> transformer
<Concatenated Operation Transformer: pipeline>
Inverse of NAD_1983_StatePlane_Massachusetts_Mainland_FIPS_2001_Feet + NAD83 to WGS 84 (1)
>>> transformer.transform(794207.7467209417, 2461029.5953208264)
(40.999999999999844, -71.0)