利用Pyproj进行坐标转换
作者：郜科科
两个坐标系统的参考椭球不同，实地一个点的不同坐标系的值是不同的，不同的部门采用的坐标系统经常是不一致，所以要转换后才能相互利用。例如目前使用的北京市观测站点位置根据GPS的定位而来，GPS使用的地理坐标系为GCS_WGS_1984，所以其坐标的地理坐标系也为GCS_WGS_1984，而假如需要将这些点显示在Web端的地图上，Web端的投影坐标系WGS_1984_Web_Mercator_Auxiliary_Sphere，就需要进行地理坐标系转换为投影坐标系的操作。
关于地理坐标系投影坐标系的关系这里不再赘述，这里介绍一下WKID。WKID的英文全称是Well Known ID，即众所周知的编号。这个编号是大家坐下来一起讨论、约定和认同的，具体有唯一性。众多的坐标系统有了自己的WKID，就像每个人都有自己的身份证号一样，从出生就定了，即使是名字改了，还是可 能通过身份证号确定，这为空间数据的使用、转换、共享等起到关键作用。
例如下表为GCS_WGS_1984的WKID格式及某点的投影文件的具体实例：
|
WKID |
4326 |
|
名称 |
GCS_WGS_1984 |
|
参数 |
GEOGCS["GCS_WGS_1984",DATUM ["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]] |
需要转换为的投影坐标同样有其WKID，例如下边为WGS_1984_Web_Mercator_Auxiliary_Sphere的具体实例：
|
WKID |
3857 |
|
名称 |
WGS_1984_Web_Mercator_Auxiliary_Sphere |
|
参数 |
PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS ["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]], PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"], PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0], PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]] |
如果对自己需要进行转换的坐标系的WKID不了解，可以从以下两个网站进行查询：
地理坐标系WKID：https://developers.arcgis.com/javascript/3/jshelp/gcs.htm
投影坐标系WKID：https://developers.arcgis.com/javascript/3/jshelp/pcs.htm
下面列举一些我国常用的地理坐标系的WKID：
|
坐标系 |
WKID |
备注 |
|
地理坐标 |
4214 |
GCS_Beijing_1954 |
|
地理坐标 |
4326 |
GCS_WGS_1984 |
|
地理坐标 |
4490 |
GCS_China_Geodetic_Coordinate_System_2000 |
|
地理坐标 |
4610 |
GCS_Xian_1980 |
进行坐标系的转换有很多工具，其中比较常用的又Proj.4相关库，如果我们使用Python进坐标转换的话，有高级的Pyproj第三方库可以使用。其文档地址如下：
http://jswhit.github.io/pyproj/
这个库非常简单，我们只需要掌握其中的一个主要函数就可以了：
transform(p1, p2, x, y, z=None, radians=False)
示例：x2, y2, z2 =transform(p1, p2, x1, y1, z1, radians=False)
这个函数表示在p1坐标系和p2坐标系之间进行坐标转换，x1，y1，z1是由p1坐标系定义的坐标，z为高度单位是米。X2，y2，z3是由p2坐标系定义的坐标，它是经过转换过后返回的，默认z1=none。Radians参数表示是否用弧度返回值。
下面我们进行一下北京市观测站点的坐标转换，如下所示为转换的代码：
|
# 投影变换 def proj_trans(): # 读取经纬度 data = pd.read_excel(u"D:/Visualization/python/file/location.xlsx") lon = data.lon.values lat = data.lat.values print lon, lat p1 = pyproj.Proj(init="epsg:4326") # 定义数据地理坐标系 p2 = pyproj.Proj(init="epsg:3857") # 定义转换投影坐标系 x1, y1 = p1(lon, lat) x2, y2 = pyproj.transform(p1, p2, x1, y1, radians=True) print x2, y2 |
在上述代码中需要注意需要在转换前首先定义数据的地理坐标系和转换后的投影坐标系，这样才能进行有目的性的转换。
转换前后的结果如下所示：
|
name |
lon |
lat |
x |
y |
|
万寿西宫 |
116.366 |
39.8673 |
12953803.87 |
4846677.374 |
|
定陵 |
116.17 |
40.2865 |
12931985.25 |
4907663.441 |
|
东四 |
116.434 |
39.9522 |
12961373.59 |
4858998.543 |
|
天坛 |
116.434 |
39.8745 |
12961373.59 |
4847721.686 |
|
农展馆 |
116.473 |
39.9716 |
12965715.05 |
4861816.127 |
|
官园 |
116.361 |
39.9425 |
12953247.27 |
4857590.051 |
|
万柳 |
116.315 |
39.9934 |
12948126.57 |
4864983.232 |
|
顺义 |
116.72 |
40.1438 |
12993210.97 |
4886860.96 |
|
怀柔 |
116.644 |
40.3937 |
12984750.68 |
4923319.714 |
|
昌平 |
116.23 |
40.1952 |
12938664.41 |
4894348.889 |
|
奥体中心 |
116.407 |
40.0031 |
12958367.96 |
4866392.773 |
|
古城 |
116.225 |
39.9279 |
12938107.82 |
4855470.429 |
下面我们使用ArcMap中的Project工具进行实验，并对照一下实验结果，如下所示为在ArcMap中进行投影转换后的结果，与上述结果基本相同。