（数据科学学习手札75）基于geopandas的空间数据分析——坐标参考系篇
（数据科学学习手札75）基于geopandas的空间数据分析——坐标参考系篇
本文对应代码已上传至我的
Github仓库https://github.com/CNFeffery/DataScienceStudyNotes
1 简介
在上一篇文章中我们对geopandas中的数据结构展开了较为全面的学习，其中涉及到面积长度等计算的过程中提到了具体的计算结果与所选择的投影坐标系关系密切，投影坐标系选择的不恰当会带来计算结果的偏差，直接关乎整个分析过程的有效与否。
作为基于geopandas的空间数据分析系列文章的第二篇，通过本文你将会学习到geopandas中的坐标参考系管理。
2 坐标参考系基础
2.1 CRS
在一个二维的平面中，我们可以使用如图1所示的坐标系统，通过坐标((x_{0},y_{0}))唯一确定点的位置：
图1
现实世界中的地球作为一个球体，当我们想要用同样的方式利用坐标((\phi_{0},\lambda_{0}))来唯一确定地球球面上的某个位置时，需要一套适应球体形状的坐标系统。而当我们想要在纸面或电脑屏幕上绘制平面地图时，就又需要有一套将地球球面展平的方法，上述的这些用于在不同情况下定义对象位置信息的坐标系统，就称为坐标参考系统（Coordinate Reference System，下文统称CRS）：
图2
CRS可细分为地理坐标系和投影坐标系。
2.1.1 地理坐标系
以弧度制下度数为单位的地理坐标系（Geographic Coordinate Systems）帮助我们定位物体在地球球面上的具体位置以及绘制球体地图：
图3 WGS84地理坐标系示意图
地理坐标系以地表上确定的某一个点为原点((0,0))，创建了包裹全球的网格，譬如WGS84，将本初子午线与赤道的交点作为原点（图5）：
图4 WGS84地理坐标系及其经纬网格
2.1.2 投影坐标系
地理坐标系虽然解决了我们在地球球面上定位的问题，但纬度和经度位置没有使用统一的测量单位，因为经度不变的情况下，纬度每变化1单位因为是对固定弧长的映射，所以真实距离是固定不变的，纬度变化1度的真实距离恒等于：
[ 2\pi\cdot地球极半径/360\approx110.95(千米) ] 可是经度每变化1单位对应的真实距离要随着纬度的变化而变化，经度变化1度的真实距离为： [ (2\pi\cdot赤道半径/360)\cdot\cos(当地纬度)\approx111.314\cos(当地纬度)(千米) ]
这就导致我们既不能直接在地理坐标系下精确度量几何对象的长度、面积，也无法直接用地理坐标系在平面上绘制出几何对象真实的形状。为了解决上述问题，各种各样的投影坐标系（Projected Coordinate Systems）被开发出来（图4，其中右下角为地理坐标系，其余均为投影坐标系）：
图5 各种CRS
投影坐标系指的是从将3D球面展平为2D平面的一套数学计算方法，利用它可以优化形状、比例/距离以及面积的失真情况，但实际情况中没有在整个地球表面都能“三全其美”的投影坐标系，有些投影坐标系优化形状上的失真，有些投影坐标系优化距离上的失真，有些投影坐标系专门针对面积失真进行优化，而有些投影坐标系可以对局部区域进行三个方面上的优化。
图6 投影坐标系变换过程示意
常用的投影坐标系如横轴墨卡托（Universal Transverse Mercator，简称UTM），基于经度将全球等分为编号0-60的区域，且每个区域又进一步细分为南半球区域或北半球区域，譬如图7所示为美国本土跨过的区域：
图7
划分出的每个区域，其原点((0,0))位于左下角顶点，距离区域中轴线500千米（图8）：
图8
针对这样划分出的独立区域利用墨卡托投影法创建各自独立的坐标网格，这个过程可以通俗地理解为用圆筒包裹地球球体，从球心发散出的光穿过球体上每个位置点投射到外部圆筒内壁从而完成3D向2D的变换：
图8
当然，这样做的后果是越靠近极点的几何对象被拉伸形变得越严重（图9），这也就是为什么俄罗斯疆域看起来如此庞大的原因：
图9 世界各国真实大小与墨卡托投影后差别
2.2 常用CRS格式
通过前文我们了解到什么是CRS，而在计算机系统中要使用CRS，需要将其文档化，下面我们来了解CRS两种常见的文档存储格式。
2.2.1 Proj4
Proj4字符串是一种识别空间或坐标参考系统的简洁方法，通过其定义的语法规则，将想要定义的CRS全部参数信息保存到一条字符串中。
- Proj4语法
Proj4字符串包含了一种CRS全部元素信息，用+连接每个元素定义部分，如下面的例子记录了横轴墨卡托北11区CRS对应的Proj4字符串：
+proj=utm +zone=11 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0它记录了如下信息：
proj=utm：声明投影方法为墨卡托 zone=11：声明对应北11区（因为这里是横轴墨卡托所以拥有独立分区，但并不是所有CRS都有分区，且在
Proj4中区号加S才为南半球分区如11S，否则默认为北半球分区） datum=WGS84：声明基准面为WGS84(基准面是椭球体用来逼近某地区用的，因此各个国家都有各自的基准面。国内常用的基准面有：BEIJING1954，XIAN1980，WGS84等) units=m：声明坐标系单位设置为米 ellps=WGS84：声明椭球面（如何计算地球的圆度）使用WGS84
上述例子记录了投影坐标系的Proj4，下面我们再来看看地理坐标系对应的Proj4，如下例：
+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0它记录了如下信息：
proj=longlat：声明这是一个地理坐标系 datum=WGS84：声明基准面为WGS84 ellps=WGS84：声明椭球面使用WGS84
与投影坐标系相比，没有单位units的信息，因为地理坐标系通常单位为十进制度数；而上述两个示例中都带有towgs84=0,0,0，这是一个转换因子，在需要进行数据转换时使用。
2.2.2 EPSG编码
EPSG（European Petroleum Survey Group）编码，使用4或5位数字编码来唯一确定已存在的一种CRS，可以在http://spatialreference.org/ref/epsg/中查看和搜索所有已知的EPSG与CRS对应关系（图10）：
图10
或在QGIS中查看：
图11
譬如对于重庆，因为地跨东经105°11~110°11，中轴线距离108E更近，常用如下投影：
图12
对应的EPSG编码为2381。
3 geopandas中的坐标参考系管理
至此，我们已经对CRS有了较为全面的了解，打好了基础，接下来我们来正式学习geopandas中的坐标参考系管理，geopandas调用pyproj作为CRS管理的后端，因此所有可以被pyproj.CRS.from_user_input()接受的合法输入同样可以被geopandas识别，譬如针对上文所说的应用于重庆区域绘图的Xian 1980 / 3-degree Gauss-Kruger CM 108E：
- Proj4
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')
图13
- EPSG
pyproj.CRS.from_user_input(2381)
图14
直接传入字符串格式的EPSG亦可：
图15
查看对应的Proj4信息，关键参数与前面Proj4一致，只是以Proj4形式传入时系统会视作创建未知CRS一样，因此相对于官方CRS缺少了一些无关紧要的其他信息：
图16
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
图17
查看其crs属性即为其对应CRS，为WGS84对应的EPSG:4326，在当前的CRS下将其绘制出来：
图18
利用to_crs()将其再投影到EPSG:2381并进行绘制：
图19
通过比较可以发现，再投影之后的中国形变失真情况得到缓解，且坐标系单位范围也发生了变化（EPSG:2381单位：米），接下来我们参考谷歌地图上点击出的重庆渝中区某地坐标：
图20
基于此创建只包含一个点的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)
图21
可以看出我们创建在重庆境内的点并没有绘制在正确的位置，接下来我们对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)
图22
这时我们定义的点被绘制到正确的位置。
同样地，可以在投影后计算更为准确的面积，这里举一个粗糙的例子（实际计算国土面积不会这样粗糙），以中国中轴线东经104.19度最靠近的105度经线对应的EPSG:2380为CRS计算面积：
图23
如果直接用原来的ESPG:4326计算面积结果如下：
图24
可以看出使用ESPG:2380计算出的面积比较接近大家记忆中的960万平方公里。