使用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（s），我查看shapefile.prj文件以确定使用了哪种投影。ESRI WKT与ESRI：102686对应，我在这里找到了它： http://epsg.io/102686 由于pyproj不认识ESRI：102686代码（会导致错误），因此我必须使用从相同网站（http://epsg.io/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

我做错了什么？

- musicformellons

Answer 1

解决了，添加了 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，则不会强制将地图投影坐标中的单位转换为米。请参见此处。

- musicformellons

---

太好了，这解决了我的问题！不过你知道为什么吗？我正在从EPSG 2230进行转换。 - ForgetfulFellow

如果可选关键字'preserve_units'为True，则地图投影坐标中的单位不会被强制为米。来源：https://jswhit.github.io/pyproj/pyproj.Proj-class.html - musicformellons

Answer 2

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

- snowman2