本文示例代码已上传至我的 GitHub 仓库 https://github.com/CNFeffery/DataScienceStudyNotes
1 简介
在本系列之前的文章中我们主要讨论了 geopandas 及其相关库在数据可视化方面的应用, 各个案例涉及的数据预处理过程也仅仅涉及到基础的矢量数据处理. 在实际的空间数据分析过程中, 数据可视化只是对最终分析结果的发布与展示, 在此之前, 根据实际任务的不同, 需要衔接很多较为进阶的空间操作, 本文就将对 geopandas 中的部分空间计算进行介绍.
本文是基于 geopandas 的空间数据分析系列文章的第 8 篇, 通过本文你将学习到 geopandas 中的空间计算(由于 geopandas 中的空间计算内容较多, 故拆分成上下两篇发出, 本文是上篇).
2 基于 geopandas 的矢量计算
geopandas 中的矢量计算根据性质的不同可分为以下几类:
2.1 构造型方法
geopandas 中的构造型方法 (Constructive Methods) 指的是从单个 GeoSeries 或 GeoDataFrame 中创建新的矢量数据的过程, 譬如早在系列第一篇文章数据结构篇中就介绍过的 bounds,exterior,interiors,boundary,centroid,convex_hull,envelope 等属性就基于 GeoSeries 计算出对应的边界, 内外轮廓线, 重心等新的矢量数据, 这些本文不再赘述, 下面我们来学习 geopandas 中常用的其他构造方法.
buffer()
geopandas 中的 buffer()方法源于 shapely, 用于缓冲区的创建, 这里给非 GIS 专业的读者朋友解释一下什么是空间意义上的缓冲区, 缓冲区用于表示点, 线, 面等矢量数据的影响范围或服务范围, 思想很简单, 即为矢量数据拓展出一定宽度的边, 图 1 展示了点, 线以及面分别对应的缓冲区的示意:
图 1
而创建缓冲区时也需要遵循一定的参数, 从而决定怎样向几何对象外进行缓冲, geopandas 中 buffer()和 shapely 中的 buffer()方法参数一致, 主要参数如下:
distance: 用于指定向外缓冲的距离, 单位与矢量数据自带单位保持一致, 在常见的投影坐标系如 web Mercator(EPSG:3857)下就是以米为单位, 因此需要注意一定要先将矢量数据转换为合适的投影坐标系之后, 再进行缓冲区分析才是合理有效的
resolution: 因为在创建缓冲区时, 对于构成矢量对象的每一个点, 都会以对应点为中心向外创建半径 = 缓冲区距离的圆, 而 Polygon 类型始终是由有限个点所构成的, 因此需要近似拼接出圆形的轮廓, resolution 参数就用于决定每个四分之一圆弧上使用多少段连续的线段来近似拼接以表示圆的形状, 默认参数值为 16, 足以近似模拟圆面积的 99.8%
下面我们分别对点, 线以及面绘制不同 resolution 参数取值下缓冲前后的对比图:
图 2
可以看出, resolution 参数对最终形成的缓冲区形态影响较大, 但默认 16 的参数下已经可以较准确地逼近圆形, 且缓冲距离还可以设置为负数, 即几何对象向内收缩:
- # 分别绘制多边形, 多边形正向缓冲区, 多边形负向缓冲区
- ax = gpd.GeoSeries([polygon,
- polygon.buffer(distance=1),
- polygon.buffer(distance=-0.25)]) \
- .plot(alpha=0.2)
- ax.axis('off')
- plt.savefig('图 3.png', dpi=300, bbox_inches='tight', pad_inches=0)
图 3
在本系列文章第一篇中介绍过 shapely 对矢量数据格式的合法性有一定规定, 如多边形不能自交叉, 可以通过 is_valid()方法判断几何对象是否合法, 而 buffer()有一个隐藏功能就是其可以通过对非法的几何对象创建距离为 0 的缓冲区来修正构成矢量对象的点的不合理连接顺序, 从而使得矢量对象变为合法的:
图 4
total_bounds
total_bounds 你应该不会感到陌生, 在前面很多篇文章中我们都使用到它来限定图像的画幅范围, 其返回依次记录了整列矢量数据所在最小矩形区域左下角 x, 左下角 y, 右上角 x 以及右上角 y 的 numpy 数组:
- geom = gpd.GeoSeries([shapely.geometry.Point([0, 0]),
- shapely.geometry.Point([0, 1]),
- shapely.geometry.Polygon([(1, 1), (1.5, 1), (1.25, 1.25)])])
- ax = geom.plot(alpha=0.4)
- # 绘制 total_bounds 范围
- ax = gpd.GeoSeries([shapely.geometry.box(*geom.total_bounds.tolist())]) \
- .plot(ax=ax,
- alpha=0.1,
- color='red')
- ax.axis('off')
图 5
simplify()
当原始的矢量数据因为形状复杂, 包含的点较多时, 会导致其文件体积较大, 如果我们需要在在线地图上叠加它们, 太大体积的矢量数据不仅会拖慢网络传输速度, 也会给图形的渲染带来更大的压力, 这时对矢量数据进行简化就非常有必要, geopandas 中沿用 shapely 中的 simplify()方法, 帮助我们对过于复杂的线和面进行简化, 和 QGIS 中简化矢量的方法一样, simplify()使用了科学的 Douglas-Peucker 算法, 基于预先设定的阈值 \(\epsilon\), 在递归判断的过程中删掉所有小于 \(\epsilon\)的点, 其过程示意如图 6:
图 6
譬如我们这里基于 - 1 到 1 之间的均匀分布, 创建一条上下波动的折线, 再用 simplify()来简化它:
- import numpy as np
- import matplotlib.patches as mpatches
- np.random.seed(10) # 固定随机数种子
- # 创建线
- line = shapely.geometry.LineString([(_, np.random.uniform(-1, 1)) for _ in range(10)])
- # 绘制简化前
- ax = gpd.GeoSeries([line]).plot(color='red')
- # 绘制简化后
- ax = gpd.GeoSeries([line]).simplify(tolerance=0.5).plot(color='blue',
- ax=ax,
- linestyle='--')
- # 制作图例映射对象列表
- LegendElement = [plt.Line2D([], [], color='red', label='简化前'),
- plt.Line2D([], [], color='blue', linestyle='--', label='简化后')]
- # 将制作好的图例映射对象列表导入 legend()中, 并配置相关参数
- ax.legend(handles = LegendElement,
- loc='lower left',
- fontsize=10)
- ax.set_ylim((-2.5, 1))
- ax.axis('off')
- plt.savefig('图 7.png', dpi=300, bbox_inches='tight', pad_inches=0)
图 7
可以看到在预设的阈值下, 对应 simplify()中的参数 tolerance=0.5, 折线得到有效地简化, 这在搭建 Web GIS 平台要渲染矢量数据时非常有效, 有效简化后的矢量数据可以在不损失太多视觉感知到的准确度的同时, 带来巨大的性能提升.
unary_union
我们都知道, 不管是 GeoSeries 还是 GeoDataFrame, 其每一行数据都代表独立的 shapely 矢量要素, 而通过 unary_union 属性, 我们可以将一整列矢量合并为单独的一个 shapely 矢量对象, 从而方便我们进行一些其他的操作:
图 8
并且如果原始数据中存在互相存在重叠的矢量对象, 通过 unary_union 之后, 返回的 shapely 对象会自动对存在重叠的矢量对象进行融合, 这一点可以方便我们的很多日常操作:
图 9
2.2 仿射变换
geopandas 中封装了几种常见的仿射变换操作, 如旋转等:
rotate()
rotate()对矢量列中的每个要素分别进行旋转操作, 其主要参数如下:
angle: 数值型, 用于指定需要旋转的角度
origin: 用于指定旋转操作的中心, 默认为 center, 是矢量对象 bbox 矩形范围的中心, centroid 表示矢量对象的重心, 或者也可以传入格式如 (x0, y0) 的坐标元组来自定义旋转中心
要注意的是 rotate()旋转方向是逆时针, 如下面的例子, 红色是旋转 90 度之后的美国:
- ax = world.query("iso_a3 =='USA'").plot(color='blue',
- alpha=0.4)
- ax = world.query("iso_a3 =='USA'").rotate(angle=90,
- origin='center') \
- .plot(color='red',
- ax=ax,
- alpha=0.4)
- plt.savefig('图 10.png', dpi=300, bbox_inches='tight', pad_inches=0)
图 10
scale()
scale()方法用于对矢量对象进行各个维度上的放缩操作, 其主要参数如下:
xfact: 数值型, 表示对 x 维度上进行放缩的因子, 默认为 1 即不放缩, 小于 1 则缩放, 大于 1 则放大
yfact: 同 xfact, 控制 y 维度上的放缩因子
origin: 同 scale()中的 origin, 用于确定缩放中心
如下面的例子, 我们在 x 以及 y 维度上对美国进行 0.5 倍放缩, 红色代表缩放之后:
- ax = world.query("iso_a3 =='USA'").plot(alpha=0.4)
- ax = world.query("iso_a3 =='USA'").scale(xfact=0.5, yfact=0.5) \
- .plot(color='red',
- alpha=0.5,
- ax=ax)
- plt.savefig('图 11.png', dpi=300, bbox_inches='tight', pad_inches=0)
图 11
translate()
translate()用于实现矢量对象的平移操作, 其主要参数有 xoff 和 yoff, 分别控制在 x 维度和 y 维度上的平移距离(与对应的投影单位保持一致):
图 12
2.3 叠加分析
geopandas 基于 shapely 中的 overlay(), 为 GeoDataFrame 赋予了同样的可以作用到整个矢量列的 overlay(), 使得我们可以对两个 GeoDataFrame 中全部的矢量对象两两之间进行基于集合关系的叠加分析(如图 13):
图 13
overlay()中的主要参数如下:
df1:GeoDataFrame, 作为输入的第一个矢量数据集
df2:GeoDataFrame, 作为输入的第二个矢量数据集
how: 字符型, 用于声明空间叠加的类型, 对应图 13, 有'intersection','union','symmetric_difference','difference', 以及额外的'identity', 他们之间的区别下文会进行详细介绍
keep_geom_type:bool 型, 当 df1 与 df2 矢量类型不同时(譬如面与线数据之间进行叠加分析), 用于决定在叠加分析产生结果中, 是否只保留与 df1 矢量类型相同的记录, 默认为 True
首先我们构造示例矢量数据, 以方便演示 overlay()不同参数下结果的区别:
- polygon1 = gpd.GeoDataFrame({
- 'value1': [1, 2],
- 'geometry': [shapely.geometry.Polygon([(1, 0), (3, 0), (3, 10), (1, 10)]),
- shapely.geometry.Polygon([(6, 0), (8, 0), (8, 10), (6, 10)])]
- })
- polygon2 = gpd.GeoDataFrame({
- 'value2': [3, 4],
- 'geometry': [shapely.geometry.Polygon([(-1, 3), (-1, 5), (10, 5), (10, 3)]),
- shapely.geometry.Polygon([(-1, 6), (-1, 8), (10, 8), (10, 6)])]
- })
- ax = polygon1.plot(color='red', alpha=0.4)
- ax = polygon2.plot(color='grey', alpha=0.4, ax=ax)
- ax.axis('off')
- plt.savefig('图 14.png', dpi=300, bbox_inches='tight', pad_inches=0)
图 14
接下来我们将其中红色部分对应的 GeoDataFrame 作为 df1, 灰色部分作为 df2, 来比较 overlay()中不同参数对应的效果:
how='union'
首先我们设置 how='union', 对 polygon1 与 polygon2 进行叠加分析:
- overlay_result = gpd.overlay(df1=polygon1,
- df2=polygon2,
- how='union')
- overlay_result
得到的结果如图 15:
图 15
可以发现, 有些行存在缺失值而有些行又是完整的, 我们分别绘制出这两类记录行:
- # 存在缺失值的行
- ax = overlay_result[overlay_result.isna().any(axis=1)] \
- .plot(color='grey')
- # 不存在缺失值的行
- ax = overlay_result[~overlay_result.isna().any(axis=1)] \
- .plot(color='blue',
- ax=ax)
- ax.set_xlim((-1, 10))
- ax.set_ylim((0, 10))
- plt.savefig('图 16.png', dpi=300, bbox_inches='tight', pad_inches=0)
图 16
在 how='union'下, 叠加分析的结果会包含所有存在相交的部分, 以及 df1 与 df2 各自剩下的不相交的部分, 如图中蓝色部分即为 df1 与 df2 相交从而不存在缺失值的部分, 而剩余的灰色部分因为没有相交, 无法获得来自另一个 GeoDataFrame 的属性值, 所以返回出来的结果会在对应的字段下填充为缺失值.
how='intersection'
在 how='intersection'参数下对 polygon1 和 polygon2 进行叠加分析:
- overlay_result = gpd.overlay(df1=polygon1,
- df2=polygon2,
- how='intersection')
- overlay_result
图 17
这时返回的结果中不再带有缺失值, 因为 intersection 只保留 df1 和 df2 彼此相交的部分:
- ax = overlay_result.plot()
- ax.set_xlim((-1, 10))
- ax.set_ylim((0, 10))
- plt.savefig('图 18.png', dpi=300, bbox_inches='tight', pad_inches=0)
图 18
how='difference'
在 how='difference'参数下对 polygon1 和 polygon2 进行叠加分析:
- overlay_result = gpd.overlay(df1=polygon1,
- df2=polygon2,
- how='difference')
- overlay_result
图 19
这时返回的结果中不再有 value2 字段, 结合图 13 可以知晓在 how='difference'下的返回结果与 Arcgis 中的擦除功能一样, 返回的是 df1 中不与 df2 相交的部分, 且以 Multi 的形式保留被切割开来的碎片矢量:
- ax = overlay_result.plot()
- ax.set_xlim((-1, 10))
- ax.set_ylim((0, 10))
- plt.savefig('图 20.png', dpi=300, bbox_inches='tight', pad_inches=0)
图 20
how='symmetric_difference'
'symmetric 的意思是对称的, 而在 how='symmetric_difference'条件下, 与 Arcgis 中的交集取反功能相同, 两个 df 中不相交的部分都会被返回:
- overlay_result = gpd.overlay(df1=polygon1,
- df2=polygon2,
- how='symmetric_difference')
- overlay_result
图 21
图 22
how='identity'
最后我们来看看 how='identity', 对应 Arcgis 中的标识功能:
- overlay_result = gpd.overlay(df1=polygon1,
- df2=polygon2,
- how='identity')
- overlay_result
图 23
为了方便理解, 我们同样分别绘制出存在缺失值的行以及不存在缺失值的行:
- # 存在缺失值的行
- ax = overlay_result[overlay_result.isna().any(axis=1)] \
- .plot(color='grey')
- # 不存在缺失值的行
- ax = overlay_result[~overlay_result.isna().any(axis=1)] \
- .plot(color='blue',
- ax=ax)
- ax.set_xlim((-1, 10))
- ax.set_ylim((0, 10))
- plt.savefig('图 24.png', dpi=300, bbox_inches='tight', pad_inches=0)
图 24
从图 24 中可以看出, 在 how='identity'条件下, 所有 df1 中不与 df2 相交的部分, 以及两者相交的部分作为返回结果, 且每个相交的部分都变为单独的要素带上所有涉及的属性字段, 而 df1 中不涉及相交的部分则仍然以 Multi 的形式被返回.
keep_geom_type
有些时候我们需要做的不仅仅是面与面之间的叠加分析. 比如在计算路网相关的指标时, 我们可能会需要与目标区域存在叠置关系的部分路网, 这就存在面与线之间的叠加分析. 参数 keep_geom_type 就用于设定最终返回的矢量数据类型是否必须与 df1 对应的类型相同, 下面我们构造示例数据来学习 keep_geom_type 参数的作用:
图 25
True 和 False 下结果如图 26 所示:
图 26
其中 GeometryCollection 类型代表多类型要素集合, 比如这里叠加分析的结果包含了一条线和一个点:
图 27
在实际工作中, 可以根据具体需要来选择使用对应的参数组合来进行叠加分析.
2.4 空间融合与拆分
有时候我们希望对矢量数据按照某些字段进行分组, 再分别对非矢量列与矢量列进行聚合及合并, 类似于 pandas 中的 groupby.agg(); 而有些时候我们希望把矢量类型为 Multi-xxx 的记录行拆分成独立要素行, 譬如将 Multi-Polygon 拆分成 Polygon 组成的若干行. 通过 geopandass 中的 dissolve()和 explode()方法, 我们就可以实现上述功能:
dissolve()
dissolve()用于对矢量数据进行融合, 可以理解为对矢量数据进行 groupby+agg 操作, 即指定的单个或多个字段值相等的分到一组, 对非矢量字段进行指定规则的聚合计算, 对矢量列进行融合, 其主要参数如下:
by: 用于指定分组所依据的字段, 单个字段传入列名字符串, 多个字段传入列名列表
aggfunc: 对分组字段外的其他非矢量列采取的聚合方式, 与 pandas 中的 agg 一致, 默认为 first, 也可以像 agg 那样传入字段和函数一一对应的字典来分别聚合不同的列
as_index:bool 型, 用于设定是否在返回的结果中将分组依据列作为索引, 默认为 True
我们以 world 数据集为例, 为了方便演示我们首先新增字段 less_than_median_gdp, 用于判断对应的国家 GDP 是否小于世界中位数水平:
图 28
接着我们以国家对应大洲列 continent 为分组依据, 并对人口和 GDP 列进行求和, 如图 29 所示, 在非矢量列得到对应的聚合计算之后, 矢量列也被融合为 Multi-Polygon:
图 29
explode()
explode()功能与 dissolve()相反, 用于将 Multi-xxx 或 Geometry-Collection 类型的数据从一行拆分到多行, 如下面的例子, 非矢量字段会自动填充到每一行:
图 30
来源: https://www.cnblogs.com/feffery/p/12909284.html