使用pyproj进行坐标转换

Question

我对地理信息系统（GIS）坐标转换不太熟悉，但我通过这个页面搞定了这个问题：http://all-geo.org/volcan01010/2012/11/change-coordinates-with-pyproj/，使用Python模块pyproj把shapefile的坐标从EPSG:28992转换成EPSG:4326，具体代码如下：

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

当我把这些坐标反向输入到谷歌地图时，它们能正确显示位置，所以这个方法是有效的。

现在我有另一个shapefile，我查看了shapefile.prj文件来确定使用了什么投影。文件里的ESRI WKT对应的是ESRI:102686，我在这里找到了相关信息：http://epsg.io/102686。由于pyproj不认识ESRI:102686这个代码（会报错），我得用proj4表示法，这个表示法我也是在同一个网站上找到的：

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)

我得到了例如这样的坐标，并把它们用在谷歌地图上：60.275122729462495, -61.873986125999316，这个位置在海洋里...

但我的结果应该是在美国的剑桥，马萨诸塞州，应该更接近于：41.00000, -71.5000000

我到底哪里出错了呢？

Answer 1

推荐使用 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)

Answer 2

问题解决了，我加上了 preserve_units = True，像这样：

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)

现在一切正常。如果可选的关键字 'preserve_units' 设置为 True，那么地图投影坐标的单位就不会被强制转换成米。详细信息可以查看 这里。