利用Pyproj进行地理投影坐标系转换

利用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中进行投影转换后的结果，与上述结果基本相同。