Python地信专题 | 基于geopandas的空间数据分析-坐标参考系篇

可是经度每变化1单位对应的真实距离要随着纬度的变化而变化，经度变化1度的真实距离为：

这就导致我们既不能直接在地理坐标系下精确度量几何对象的长度、面积，也无法直接用地理坐标系在平面上绘制出几何对象真实的形状。

为了解决上述问题，各种各样的投影坐标系（Projected Coordinate Systems）被开发出来（图5，其中右下角为地理坐标系，其余均为投影坐标系）：

投影坐标系指的是从将3D球面展平为2D平面的一套数学计算方法，利用它可以优化形状、比例/距离以及面积的失真情况。

但实际情况中没有在整个地球表面都能“三全其美”的投影坐标系，有些投影坐标系优化形状上的失真，有些投影坐标系优化距离上的失真，有些投影坐标系专门针对面积失真进行优化，而有些投影坐标系可以对局部区域进行三个方面上的优化。

常用的投影坐标系如横轴墨卡托（Universal Transverse Mercator，简称UTM），基于经度将全球等分为编号0-60的区域，且每个区域又进一步细分为南半球区域或北半球区域，譬如图7所示为美国本土跨过的区域：

划分出的每个区域，其原点位于左下角顶点，距离区域中轴线500千米（图8）：

针对这样划分出的独立区域利用墨卡托投影法创建各自独立的坐标网格，这个过程可以通俗地理解为用圆筒包裹地球球体，从球心发散出的光穿过球体上每个位置点投射到外部圆筒内壁从而完成3D向2D的变换：

当然，这样做的后果是越靠近极点的几何对象被拉伸形变得越严重（图10），这也就是为什么俄罗斯疆域看起来如此庞大的原因：

2.2 常用CRS格式

通过前文我们了解到什么是CRS，而在计算机系统中要使用CRS，需要将其文档化，下面我们来了解CRS两种常见的文档存储格式。

2.2.1 Proj4

Proj4字符串是一种识别空间或坐标参考系统的简洁方法，通过其定义的语法规则，将想要定义的CRS全部参数信息保存到一条字符串中。

Proj4字符串包含了一种CRS全部元素信息，用+连接每个元素定义部分，如下面的例子记录了横轴墨卡托北11区CRS对应的Proj4字符串：

+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0

它记录了如下信息：

上述例子记录了投影坐标系的Proj4，下面我们再来看看地理坐标系对应的Proj4，如下例：

+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0

它记录了如下信息：

与投影坐标系相比，没有单位units的信息，因为地理坐标系通常单位为十进制度数；而上述两个示例中都带有towgs84=0,0,0，这是一个转换因子，在需要进行数据转换时使用。

2.2.2 EPSG编码

EPSG（European Petroleum Survey Group）编码，使用4或5位数字编码来唯一确定已存在的一种CRS，可以在http://spatialreference.org/ref/epsg/中查看和搜索所有已知的EPSG与CRS对应关系（图11）：

或在QGIS中查看：

譬如对于重庆，因为地跨东经105°11~110°11，中轴线距离108E更近，常用如下投影：

对应的EPSG编码为2381。

3 geopandas中的坐标参考系管理

至此，我们已经对CRS有了较为全面的了解，打好了基础，接下来我们来正式学习geopandas中的坐标参考系管理。

geopandas调用pyproj作为CRS管理的后端，因此所有可以被pyproj.CRS.from_user_input()接受的合法输入同样可以被geopandas识别，譬如针对上文所说的应用于重庆区域绘图的Xian 1980 / 3-degree Gauss-Kruger CM 108E：

import pyproj

pyproj.CRS.from_user_input('+proj=tmerc +lat_0=0 +lon_0=108 +k=1 +x_0=500000 +y_0=0 +ellps=IAU76 +units=m +no_defs')

pyproj.CRS.from_user_input(2381)

直接传入字符串格式的EPSG亦可：

查看对应的Proj4信息，关键参数与前面Proj4一致，只是以Proj4形式传入时系统会视作创建未知CRS一样，因此相对于官方CRS缺少了一些无关紧要的其他信息：

3.1 CRS的设置与再投影

在上一篇文章（数据科学学习手札74）基于geopandas的空间数据分析——数据结构篇中我们介绍了创建GeoSeries和GeoDataFrame的方法。

实际上，现实的空间分析计算任务中，必须要为数据设置合适的CRS，在geopandas.GeoSeries()和geopandas.GeoDataFrame()中就包含参数crs。

下面我们举例说明，还是先用到geopandas自带的世界国家地区数据，我们从中选择中国（坚持一个中国，我们将台湾地区组合进国土中）：

import geopandas as gpd

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
# 利用name字段选择中国区域
china = world.loc[world['name'].isin(['China', 'Taiwan'])]
china

查看其crs属性即为其对应CRS，为WGS84对应的EPSG:4326，在当前的CRS下将其绘制出来：

利用to_crs()将其再投影到EPSG:2381并进行绘制：

通过比较可以发现，再投影之后的中国形变失真情况得到缓解，且坐标系单位范围也发生了变化（EPSG:2381单位：米），接下来我们参考谷歌地图上点击出的重庆渝中区某地坐标：

基于此创建只包含一个点的GeoSeries，尝试将其与EPSG:2381下的中国地图一同绘制：

from shapely import geometry
import matplotlib.pyplot as plt

cq = gpd.GeoSeries([geometry.Point([106.561203, 29.558078])],
              crs='EPSG:4326')

fig, ax = plt.subplots()
china.to_crs(crs='EPSG:2381').plot(ax=ax, color='red', alpha=0.8)
cq.plot(ax=ax, color='orange', markersize=100, marker='x')
plt.xticks(rotation=20)

可以看出我们创建在重庆境内的点并没有绘制在正确的位置，接下来我们对cq进行再投影，再尝试将其与EPSG:2381下的中国绘制在一起：

fig, ax = plt.subplots()
china.to_crs(crs='EPSG:2381').plot(ax=ax, color='red', alpha=0.8)
# 先再投影到EPSG:2381
cq.to_crs(crs='EPSG:2381').plot(ax=ax, color='orange', markersize=100, marker='x')
plt.xticks(rotation=20)