coordinates conversion with pyproj

Question

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?

Answer 1

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.