python基于geopandas的空间数据分析之二：坐标参考系

AndersonHuang

已于 2023-11-01 16:24:32 修改

阅读量470

点赞数 1

分类专栏：数据挖掘数据分析 GIS 文章标签： python 数据分析开发语言

于 2023-10-30 09:45:59 首次发布

本文链接：https://blog.csdn.net/andersonhuang/article/details/134040250

版权

数据挖掘同时被 3 个专栏收录

22 篇文章 0 订阅

订阅专栏

数据分析

15 篇文章 0 订阅

订阅专栏

GIS

6 篇文章 0 订阅

订阅专栏

1 简介

在geopandas中，涉及到面积长度等计算的结果与所选择的投影坐标系关系密切，投影坐标系选择的不恰当会带来计算结果的偏差，直接关乎整个分析过程的有效与否。本文将探索geopandas对坐标参考系的管理内容。

2 坐标参考系

2.1 CRS的基本含义

在一个二维的平面中，我们可以使用下图所示的坐标系统，通过坐标 (x₀, y₀) 唯一确定点的位置：

图1 坐标系统

现实世界中的地球是一个球体，当我们想要用同样的方式利用坐标 (ϕ₀, λ₀) 来唯一确定地球球面上的某个位置时，需要一套适应球体形状的坐标系统。而当我们想要在纸面或电脑屏幕上绘制平面地图时，就又需要有一套将地球球面展平的方法。上述的这些用于在不同情况下定义对象位置信息的坐标系统，就称为 坐标参考系统（Coordinate Reference System），下文统称 CRS ：
图2 坐标参考系统

图2 坐标参考系统

CRS 可细分为地理坐标系和投影坐标系。

2.1.1 地理坐标系

以弧度制下度数为单位的 地理坐标系（Geographic Coordinate Systems） 帮助我们定位物体在地球球面上的具体位置以及绘制球体地图：
图3 WGS84地理坐标系

图3 WGS84地理坐标系

地理坐标系以地表上确定的某一个点为 原点(0, 0)，创建了包裹全球的网格。譬如WGS84，将本初子午线与赤道的交点作为原点：

图4 WGS84地理坐标系及其经纬网格

2.1.2 投影坐标系

地理坐标系 虽然解决了我们在地球球面上定位的问题，但纬度和经度位置没有使用统一的测量单位，因为经度不变的情况下，纬度每变化1单位因为是对固定弧长的映射，所以真实距离是固定不变的，纬度变化1度的真实距离恒等于：
$2Π \cdot 地球半径 /360 \approx 110.95 (千米)$
可是经度每变化1单位对应的真实距离要随着纬度的变化而变化，经度变化1度的真实距离为
$(2Π \cdot 赤道半径 /360) \cdot cos (当地纬度) \approx 111.314 cos (当地纬度) (千米)$

这就导致我们既不能直接在地理坐标系下精确度量几何对象的长度、面积，也无法直接用地理坐标系在平面上绘制出几何对象真实的形状。为了解决上述问题，各种各样的 投影坐标系（Projected Coordinate Systems） 被开发出来（图5，其中右下角为地理坐标系，其余均为投影坐标系）：

图5 各种CRS

投影坐标系 指的是从将3D球面展平为2D平面的一套数学计算方法，利用它可以优化 形状、比例/距离以及面积 的失真情况，但实际情况中没有在整个地球表面都能 “三全其美” 的投影坐标系，有些投影坐标系优化形状上的失真，有些投影坐标系优化距离上的失真，有些投影坐标系专门针对面积失真进行优化，而有些投影坐标系可以对局部区域进行三个方面上的优化。

图6 投影坐标系变换过程

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

图7 横轴墨卡托

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

图8 横轴墨卡托原点

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

图9 墨卡托投影法

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

图10 世界各国真实大小与墨卡托投影后差别

2.2 常用CRS格式

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

2.2.1 Proj4字符串

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

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

+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，中文名称为欧洲石油调查组织。EPSG（European Petroleum Survey Group）编码，是使用4或5位数字编码来唯一确定已存在的一种 CRS，可以在 http://spatialreference.org/ref/epsg/ 中查看和搜索所有已知的EPSG与CRS对应关系（图11）：

图11

或在QGIS软件中查看：

图12

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

图13
对应的EPSG编码为2381。

2.2.3 EPSG编码的选用

非常重要！！！ 地理空间分析中的大多数错误来自于为所需操作选择了错误的CRS。常见的CRS误区：

混合坐标系：当结合数据集时，空间对象必须具有相同的参考系统。请确保将所有的东西转换为相同的CRS。我们在下面向你展示如何进行这种转换。
计算面积：在测量一个形状的面积之前，使用一个等面积的CRS。
计算距离：在计算物体之间的距离时，使用等距CRS。

常用EPSG编码

EPSG:4326 (WGS84)
在世界地图方面，EPSG:4326是比较著名的一个，加密前的高德、百度用的也是WGS84，因为由美国主导的GPS系统就是在用它，它还有一个名气更大的别名叫作WGS84，WGS(World Geodetic System)是世界大地测量系统的意思，由于是1984年定义的，所以叫WGS84，之前的版本还有WGS72、WGS66、WGS60。
EPSG:3857(墨卡托投影)
另一个比较知名的编码是EPSG:3857，这也是一张世界地图，目前主要是各大互联网地图公司以它为基准，例如Google地图，Microsoft地图都在用它。
GCJ-02坐标系
GCJ-02是由中国国家测绘局（G表示Guojia国家，C表示Cehui测绘，J表示Ju局）制订的地理信息系统的坐标系统。它是一种对经纬度数据的加密算法，即加入随机的偏差。国内出版的各种地图系统（包括电子形式），必须至少采用GCJ-02对地理位置进行首次加密。
EPSG：4490（CGCS2000）
我国的GPS系统-北斗导航系统以及国家发行的“天地图”，用的是这一套地理坐标系统，中文名“中国国家2000地理坐标系统”，英文全称翻译名“中国大地坐标系2000”。英文名 China Geodetic Coordinate System 2000。
北京54、西安80
是我国已经逐渐停止使用的两个地理坐标系统。北京54坐标系统WKID是4214，西安80坐标系统的WKID是4610。

2.2.4 小结

上述坐标系的关系可以用下面的脑图来大概表示
图25

google的摩卡托坐标，也就是我们经常看到的 EPSG:3857 坐标系。
EPSG:3857 的数据一般是这种的。[12914838.35,4814529.9]，看上去相对数值较大。不利于存储，比较占内存。
4326 WGS-84：是国际标准，GPS坐标（Google Earth使用、或者GPS模块）
EPSG:4326 的数据一般是这种的。[22.37，114.05]。利于存储，可读性高
所以我们常常看到和用到的坐标系数据往往不是墨卡托坐标，而是EPSG:4326坐标系下的坐标数据。

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')

图14

EPSG

pyproj.CRS.from_user_input(2381) # 传入编码数字
#pyproj.CRS.from_user_input("EPSG:2381") # 传入字符串

图15

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

pyproj.CRS.from_user_input(2381).to_proj4()

图16

3.1 CRS的设置与再投影

在现实的空间分析计算任务中，创建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下将其绘制出来：

# 查看CRS
china.crs

图18

%matplotlib widget
china.plot(color='red', alpha=0.8)

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

%matplotlib widget
china.to_crs("EPSG:2381").plot(color='red', alpha=0.8)

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

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

from shapely import geometry
import matplotlib.pyplot as plt

cq = gpd.GeoSeries(data=[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计算面积结果如下（大家注意这里的计算其实是毫无意义的，因为EPSG:4326是地理坐标系，单位是度，无任何面积计算意义！）：

图24
可以看出使用ESPG:2380计算出的面积比较接近大家记忆中的960万平方公里。

参考文章:
1 （数据科学学习手札75）基于geopandas的空间数据分析——坐标参考系篇

2 地理空间数据分析入门：GeoPandas 和 Shapely

3 geopandas–坐标系学习

4 地图坐标系总结

AndersonHuang

关注

1
点赞
踩
2

收藏

觉得还不错? 一键收藏
0
评论
python基于geopandas的空间数据分析之二：坐标参考系

在geopandas中，涉及到面积长度等计算的结果与所选择的投影坐标系关系密切，投影坐标系选择的不恰当会带来计算结果的偏差，直接关乎整个分析过程的有效与否。本文将探索geopandas对坐标参考系的管理内容。
复制链接

扫一扫