Transforming coordinates with pyproj

Question

I do not experience the conversion of gis coordinates, but I control using this page: http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/ to convert the coordinates of a shapefile from EPSG: 28992 to EPSG : 4326 using the pyproj python module using these operators:

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

When I go back and enter these coordinates on Google maps, they give me the correct locations. So it works fine.

Now I have another shapefile (s) and I look at the shapefile.prj file to determine which projection was used. ESRI WKT matches ESRI: 102686, which I find here: http://epsg.io/102686 Since the code ESRI: 102686 is not known to pyproj (gives an error), I have to use the proj4 notation 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, for example. coordinates and use them on google maps: 60.275122729462495, -61.873986125999316 which is located somewhere in the ocean ...

But my results should be in Cambridge, Massachusetts in the USA, so 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 great. If the optional preserve_units keyword is True, the units in the map projection coordinates must not be meters. See here .

Answer 2

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)