geopandas学习（七） 空间计算

1. 基于geopandas的矢量计算

geopandas中的矢量计算根据性质的不同可分为以下几类：

1.1 构造型方法

geopandas中的构造型方法（Constructive Methods）指的是从单个GeoSeries或GeoDataFrame中创建新的矢量数据的过程。

譬如早在系列第一篇文章数据结构篇中就介绍过的bounds、exterior、interiors、boundary、centroid、convex_hull、envelope等属性，就基于GeoSeries计算出对应的边界、内外轮廓线、重心等新的矢量数据，这些本文不再赘述，下面我们来学习geopandas中常用的其他构造方法。

- buffer()

geopandas中的buffer()方法源于shapely，用于缓冲区的创建，这里给非GIS专业的读者朋友解释一下什么是空间意义上的缓冲区：

缓冲区用于表示点、线、面等矢量数据的影响范围或服务范围，思想很简单，即为矢量数据拓展出一定宽度的边，图1展示了点、线以及面分别对应的缓冲区的示意：

而创建缓冲区时也需要遵循一定的参数，从而决定怎样向几何对象外进行缓冲，geopandas中buffer()和shapely中的buffer()方法参数一致，主要参数如下：

• distance：用于指定向外缓冲的距离，单位与矢量数据自带单位保持一致，在常见的投影坐标系如Web Mercator（EPSG:3857）下就是以米为单位，因此需要注意一定要先将矢量数据转换为合适的投影坐标系之后，再进行缓冲区分析才是合理有效的
• resolution：因为在创建缓冲区时，对于构成矢量对象的每一个点，都会以对应点为中心向外创建半径=缓冲区距离的圆，而Polygon类型始终是由有限个点所构成的，因此需要近似拼接出圆形的轮廓，resolution参数就用于决定每个四分之一圆弧上使用多少段连续的线段来近似拼接以表示圆的形状，默认参数值为16，足以近似模拟圆面积的99.8%

下面我们分别对点、线以及面绘制不同resolution参数取值下缓冲前后的对比图：

fig, axes = plt.subplots(3, 3, figsize=(9, 9))

axes = axes.flatten()

pt = shapely.geometry.Point([0, 0])
ls = shapely.geometry.LineString([(0, 0), (1, 0), (1, 1)])
polygon = shapely.geometry.Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])

axes[0] = gpd.GeoSeries([pt,
                         pt.buffer(distance=0.2, resolution=1)]) \
             .plot(ax=axes[0],
                   alpha=0.2)
axes[0].axis('off')
axes[0].set_title('点：resolution=3')

axes[1] = gpd.GeoSeries([pt,
                         pt.buffer(distance=0.2, resolution=4)]) \
             .plot(ax=axes[1],
                   alpha=0.2)
axes[1].axis('off')
axes[1].set_title('点：resolution=16')

axes[2] = gpd.GeoSeries([pt,
                         pt.buffer(distance=0.2, resolution=16)]) \
             .plot(ax=axes[2],
                   alpha=0.2)
axes[2].axis('off')
axes[2].set_title('点：resolution=16')

axes[3] = gpd.GeoSeries([ls,
                         ls.buffer(distance=0.2, resolution=1)]) \
             .plot(ax=axes[3],
                   alpha=0.2)
axes[3].axis('off')
axes[3].set_title('线：resolution=1')

axes[4] = gpd.GeoSeries([ls,
                         ls.buffer(distance=0.2, resolution=4)]) \
             .plot(ax=axes[4],
                   alpha=0.2)
axes[4].axis('off')
axes[4].set_title('线：resolution=4')

axes[5] = gpd.GeoSeries([ls,
                         ls.buffer(distance=0.2, resolution=16)]) \
             .plot(ax=axes[5],
                   alpha=0.2)
axes[5].axis('off')
axes[5].set_title('线：resolution=16')

axes[6] = gpd.GeoSeries([polygon,
                         polygon.buffer(distance=0.2, resolution=1)]) \
             .plot(ax=axes[6],
                   alpha=0.2)
axes[6].axis('off')
axes[6].set_title('面：resolution=1')

axes[7] = gpd.GeoSeries([polygon,
                         polygon.buffer(distance=0.2, resolution=4)]) \
             .plot(ax=axes[7],
                   alpha=0.2)
axes[7].axis('off')
axes[7].set_title('面：resolution=4')

axes[8] = gpd.GeoSeries([polygon,
                         polygon.buffer(distance=0.2, resolution=16)]) \
             .plot(ax=axes[8],
                   alpha=0.2)
axes[8].axis('off')
axes[8].set_title('面：resolution=16')

plt.savefig('图2.png', dpi=300, bbox_inches='tight', pad_inches=0)

可以看出，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)

在本系列文章第一篇中介绍过shapely对矢量数据格式的合法性有一定规定，如多边形不能自交叉，可以通过is_valid()方法判断几何对象是否合法。

而buffer()有一个隐藏功能就是其可以通过对非法的几何对象创建距离为0的缓冲区来修正构成矢量对象的点的不合理连接顺序，从而使得矢量对象变为合法的：

polygon = shapely.geometry.Polygon([(0, 0), (1, 1), (1, 0), (0, 1)])
print(polygon.is_valid)
polygon

print(polygon.buffer(0).is_valid)
polygon

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

plt.savefig('图5.png', dpi=300, bbox_inches='tight', pad_inches=0)

- simplify()

当原始的矢量数据因为形状复杂，包含的点较多时，会导致其文件体积较大，如果我们需要在在线地图上叠加它们，太大体积的矢量数据不仅会拖慢网络传输速度，也会给图形的渲染带来更大的压力。

这时对矢量数据进行简化就非常有必要，geopandas中沿用shapely中的simplify()方法，帮助我们对过于复杂的线和面进行简化，和QGIS中简化矢量的方法一样，simplify()使用了科学的Douglas-Peucker算法，基于预先设定的阈值，在递归判断的过程中删掉所有小于的点，其过程示意如图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)

- unary_union

我们都知道，不管是GeoSeries还是GeoDataFrame，其每一行数据都代表独立的shapely矢量要素，而通过unary_union属性，我们可以将一整列矢量合并为单独的一个shapely矢量对象，从而方便我们进行一些其他的操作：

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.geometry.unary_union

type(world.geometry.unary_union)

并且如果原始数据中存在互相存在重叠的矢量对象，通过unary_union之后，返回的shapely对象会自动对存在重叠的矢量对象进行融合，这一点可以方便我们的很多日常操作：

polygons = gpd.GeoSeries([shapely.geometry.Polygon([(0, 0), (1, 0), (1, 1), (0, 1)]),
                          shapely.geometry.Polygon([(-0.5, 0), (0.5, 0), (0.5, 1), (-0.5, 1)])])
ax = polygons.plot(alpha=0.2)
ax.axis('off')

1.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)

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

- translate()

translate()用于实现矢量对象的平移操作，其主要参数有xoff和yoff，分别控制在x维度和y维度上的平移距离（与对应的投影单位保持一致）:

gpd.GeoSeries([shapely.geometry.Point([0, 0])]) \
   .translate(xoff=1, yoff=1)

1.3 叠加分析

geopandas基于shapely中的overlay()，为GeoDataFrame赋予了同样的可以作用到整个矢量列的overlay()，使得我们可以对两个GeoDataFrame中全部的矢量对象两两之间进行基于集合关系的叠加分析（如图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)

接下来我们将其中红色部分对应的GeoDataFrame作为df1，灰色部分作为df2，来比较overlay()中不同参数对应的效果：

- how=‘union’

首先我们设置how=‘union’，对polygon1与polygon2进行叠加分析：

overlay_result = gpd.overlay(df1=polygon1,
                             df2=polygon2,
                             how='union')
overlay_result

# 存在缺失值的行
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)

在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

这时返回的结果中不再带有缺失值，因为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)

- how=‘difference’

在how='difference’参数下对polygon1和polygon2进行叠加分析：

overlay_result = gpd.overlay(df1=polygon1,
                             df2=polygon2,
                             how='difference')
overlay_result

这时返回的结果中不再有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)

- how=‘symmetric_difference’

'symmetric的意思是对称的，而在how='symmetric_difference’条件下，与Arcgis中的交集取反功能相同，两个df中不相交的部分都会被返回：

overlay_result = gpd.overlay(df1=polygon1,
                             df2=polygon2,
                             how='symmetric_difference')
overlay_result

ax = overlay_result.plot()

ax.set_xlim((-1, 10))
ax.set_ylim((0, 10))
plt.savefig('图22.png', dpi=300, bbox_inches='tight', pad_inches=0)

- how=‘identity’

最后我们来看看how=‘identity’，对应Arcgis中的标识功能：

overlay_result = gpd.overlay(df1=polygon1,
                             df2=polygon2,
                             how='identity')
overlay_result

为了方便理解，我们同样分别绘制出存在缺失值的行以及不存在缺失值的行：

# 存在缺失值的行
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中可以看出，在how='identity’条件下，所有df1中不与df2相交的部分，以及两者相交的部分作为返回结果，且每个相交的部分都变为单独的要素带上所有涉及的属性字段，而df1中不涉及相交的部分则仍然以Multi的形式被返回。

- keep_geom_type

有些时候我们需要做的不仅仅是面与面之间的叠加分析。

比如在计算路网相关的指标时，我们可能会需要与目标区域存在叠置关系的部分路网，这就存在面与线之间的叠加分析。

参数keep_geom_type就用于设定最终返回的矢量数据类型是否必须与df1对应的类型相同，下面我们构造示例数据来学习keep_geom_type参数的作用：

# 构造面
polygon = gpd.GeoDataFrame({
    'value1': [1],
    'geometry': [shapely.geometry.Polygon([(0, 0), (1, 0), (1, 1), (0, 1)])]
})

# 构造线
line = gpd.GeoDataFrame({
    'value2': [2],
    'geometry': [shapely.geometry.LineString([(1, 0.5), (-1, 0.5), (-1, 1), (0, 1)])]
})

ax = polygon.plot(alpha=0.2)
ax = line.plot(alpha=0.2, ax=ax)

True和False下结果如图26所示：

其中GeometryCollection类型代表多类型要素集合，比如这里叠加分析的结果包含了一条线和一个点：

# 默认参数：keep_geom_type=False
gpd.overlay(df1=polygon,
            df2=line,
            keep_geom_type=False).unary_union

1.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是否小于世界中位数水平：

median = world['gdp_md_est'].median()
world['less_than_median_gdp'] = world['gdp_md_est'].apply(lambda x: 1 if x < median else 0)
world.head()

接着我们以国家对应大洲列continent为分组依据，并对人口和GDP列进行求和，如图29所示，在非矢量列得到对应的聚合计算之后，矢量列也被融合为Multi-Polygon：

world.dissolve(by='continent', aggfunc={'pop_est': 'sum',
                                        'gdp_md_est': 'sum'})

- explode()

explode()功能与dissolve()相反，用于将Multi-xxx或Geometry-Collection类型的数据从一行拆分到多行，如下面的例子，非矢量字段会自动填充到每一行：

gpd.GeoDataFrame({
    'value1': [1],
    'value2': ['a'],
    'geometry': [shapely.geometry.MultiPoint([(0, 0), (1, 1), (2, 3)])]
})

gpd.GeoDataFrame({
    'value1': [1],
    'value2': ['a'],
    'geometry': [shapely.geometry.MultiPoint([(0, 0), (1, 1), (2, 3)])]
}).explode()

1.5 空间连接

类比常规表格数据的连接操作，在空间数据分析中也存在类似表连接的操作，譬如我们手头有一张包含设施点数据的矢量表，以及另一张包含行政区划面数据的矢量表，当我们想要通过某些操作来统计出每个行政区划面内部的设施点信息时，空间连接就可以非常方便快捷地实现这类需求。

我们都清楚常规表格数据的连接，是按照设定的连接方式，将每张表中指定的某列或某些列数值相等的记录行合并为同一行，最后汇整成连接结果表返回：

而空间连接不同于常规表连接，其合并同一行的依据不是检查指定的列数值是否相等，而是基于不同矢量表其矢量列之间的空间拓扑关系，譬如相交、包含等。

在geopandas中我们利用sjoin函数来实现空间连接，其使用方式类似pandas中的merge接近，主要参数如下：

• left_df：GeoDataFrame，传入空间连接对应的左表
• right_df：GeoDataFrame，传入空间连接对应的右表
• how：字符型，用于决定连接方式，'inner’表示内连接，且连接结果表中的矢量列来自左表；'left’表示左连接，且结果表中的矢量列来自左表；'right’表示右连接，最终结果表中的矢量列来自右表
• op：字符型，用于设定拓扑判断的规则，'intersects’代表相交，即几何对象之间存在共有的边或内部点；'contains’代表包含，即一个几何对象至少有一个点位于另一个几何对象内部，且其本身没有任何点落在另一个结几何对象的外部；‘within’表示在内部，是’contains’的相反情况，即左表被右表矢量’contains’
• lsuffix：字符型，代表当左右表连接之后存在重名列时，为左表重名的列添加的后缀，默认为’left’
• rsuffix：字符型，意义类似lsuffix，默认为’right’

了解过sjoin()中的核心参数后，我们来通过实际例子理解它们的具体作用，how的作用与pandas中效果的一致，这里不多解读，我们来重点学习op各参数的不同效果：

• 参数op

intersects是空间连接中最常使用的模式，即相比较的两个几何对象有至少1个公共点就会被匹配上，下面我们以柏林公交站点数据为例，首先我们先读入柏林行政区划面数据，其中字段Gemeinde_n是每个行政区划的名称：

# 读入柏林行政区划面文件
Berlin = gpd.read_file('Berlin/Bezirke__Berlin.shp')
Berlin.head() # Gemeinde_n代表镇，即Berlin中每个面文件对应的行政区划名称

接着再读入柏林全部交通车站数据，其中fclass列代表对应车站的类别：

Berlin_transport = gpd.read_file('Berlin/gis_osm_transport_free_1.shp')
Berlin_transport.head()

对站点的空间分布进行简单的可视化：

fig, ax = plt.subplots(figsize=(8, 8))

# 投影选择web mercator 
ax.axis('off')
# 区划图层
ax = Berlin.to_crs('EPSG:3857') \
           .plot(ax=ax,
                 facecolor='none',
                 edgecolor='grey')
# 站点图层
ax = Berlin_transport.to_crs('EPSG:3857') \
                     .plot(ax=ax,
                           facecolor='lightgrey',
                           alpha=0.4,
                           edgecolor='grey',
                           markersize=6)

fig.savefig('图5.png', dpi=300, bbox_inches='tight', pad_inches=0)

接着我们就利用sjoin()将区划面作为左表，站点作为右表，在op='intersects’参数设置下进行空间连接，再衔接groupby，以统计出各区划面内部的公交站点数量：

gpd.sjoin(left_df=Berlin,
          right_df=Berlin_transport.query("fclass=='bus_stop'"),
          op='intersects') \
   .groupby('Gemeinde_n') \
   .size()

再设置op=‘contains’，因为进行连接的对象是左表面要素，右表点要素，所以这里的效果等价于op=‘intersects’：

gpd.sjoin(left_df=Berlin,
          right_df=Berlin_transport.query("fclass=='bus_stop'"),
          op='contains') \
   .groupby('Gemeinde_n') \
   .size()

但当op='within’时，按照拓扑规则，如果依旧是左表面要素，右表点要素，得到的结果就会为空，反过来则正常：

gpd.sjoin(left_df=Berlin_transport.query("fclass=='bus_stop'"),
          right_df=Berlin,
          op='within') \
   .groupby('Gemeinde_n') \
   .size()

类似的，其他类型几何对象之间的空间连接你也可以根据自己的需要进行操作，值得一提的是，利用sjoin()进行空间左、右、内连接时，因为结果表依旧是GeoDataFrame，所以只会保留一列矢量列，按照上文中参数介绍部分的描述，只有右连接时结果表中的矢量列才来自右表，但无论采取什么连接方式，结果表中未被保留的矢量列对应的index会被作为单独的一列保存下来，帮助我们可以按图索骥利用loc方式索引出需要的数据：

1.6 拓扑关系判断

geopandas中除了在上一篇文章中介绍的叠加分析以及上文介绍的空间连接中基于拓扑关系判断实现多表数据联动之外，还针对GeoSeries与GeoDataFrame设计了一系列方法，可以直接进行矢量数据之间的拓扑关系判断并返回对应的bool型判断结果，以contains()为例，在比较矢量数据之间拓扑关系时，矢量数据与待比较矢量数据之间主要有以下几种格式：

• 长度n与长度1进行比较

当主体矢量列长度为n，而输入待比较的矢量列长度为1时，返回的bool值是待比较矢量列与主题矢量列一一进行比较后的结果：

geo1 = gpd.GeoSeries([Polygon([(-1, -1), (1, -1), (1, 1), (-1, 1)]),
                      Polygon([(1, 1), (2, 1), (2, 2), (1, 2)])])
geo2 = gpd.GeoSeries([Point(0, 0)])

geo1.contains(geo2)

• 长度1与长度n进行比较

与前面一种情况类似，只不过这里是将主体矢量列与待比较矢量列一一比较之后的结果：

geo1 = gpd.GeoSeries([Polygon([(-1, -1), (1, -1), (1, 1), (-1, 1)])])
geo2 = gpd.GeoSeries([Point(0, 0), Point(2, 2)])

geo1.contains(geo2)

• 长度m与长度m-n（n>0）进行比较

这里所说的情况指主体矢量与待比较矢量长度都不为1，且主体矢量列的长度大于待比较矢量，这时返回的结果只会对主体矢量列前m-n个要素与待比较矢量对应位置一一比较，主体矢量被截断未能进行比较的部分默认返回False：

geo1 = gpd.GeoSeries([Polygon([(-1, -1), (1, -1), (1, 1), (-1, 1)]),
                      Polygon([(1.5, 1.5), (2.5, 1.5), (2.5, 2.5), (1.5, 2.5)]),
                      Polygon([(3, 3), (4, 3), (4, 4), (3, 4)])])
geo2 = gpd.GeoSeries([Point(0, 0), Point(2, 2)])

geo1.contains(geo2)

geopandas中进行拓扑关系判断的基本原则了解完了，下面罗列出常用的一些拓扑关系判断API，均为GeoSeries或GeoDataFrame的方法：

• intersects()：检查相交关系
• contains()：检查包含关系，即主体矢量完全包裹住待比较的矢量且它们的边界互不接触，譬如面对点的包含
• within()：检查主体矢量是否在待检查矢量的内部
• touches()：检查触碰关系，即两个矢量之间至少有一个1个公共点，但它们的内部无任何相交区域
• crosses()：检查交叉关系，常见如线与线之间的交叉
• disjoint()：检查不相交关系，即两个矢量之间没有任何接触
• geom_equals()：检查是否完全相同
• overlaps()：检查重叠关系

1.7 空间裁切

在空间数据分析中，裁切也是非常常用的操作，譬如我们想要获取某个公交站周围500米半径内部的路网矢量，就可以使用到裁切。

在geopandas中我们可以使用clip()函数来基于蒙版矢量对目标矢量进行裁切，其主要参数如下：

• gdf：GeoDataFrame或GeoSeries，代表将要被裁切的矢量数据集
• mask：GeoDataFrame、GeoSeries或shapely中的Polygon、Multi-Polygon对象，代表蒙版矢量
• keep_geom_type：同叠加分析overlay中的同名参数

基于实际例子进行演示，我们读入数据berlin_footway_WGS84.shp，包含了柏林全部的步道路网线数据，并转换到适合柏林地区的投影EPSG:32633：

Berlin_footway = gpd.read_file('Berlin/berlin_footway_WGS84.shp')
Berlin_footway.to_crs('EPSG:32633', inplace=True) # 转换为合适的投影坐标系
Berlin_footway.plot()

接下来我们从上文中使用到的柏林车站点数据中筛选出租车站点，与步道路网数据统一坐标参考系，生成500米缓冲区，并利用上一篇文章中介绍过的unary_union来得到MultiPolygon对象：

# 整合所有出租车站点500米缓冲区
taxi_station_500_buffer = Berlin_transport.query("fclass=='taxi'") \
                                          .to_crs('EPSG:32633') \
                                          .buffer(distance=500) \
                                          .unary_union
taxi_station_500_buffer

万事俱备，接下来我们使用clip()来裁切所有出租车站点500米缓冲区内部的步行道路网：

# 裁切所有出租车站点500米缓冲区内部的路网线数据
taxi_station_500buffer_roads = gpd.clip(gdf=Berlin_footway,
                                        mask=taxi_station_500_buffer)

在交互模式下同时绘制出缓冲区以及裁切出的路网：

# 交互式绘图模式
%matplotlib widget 

ax = taxi_station_500buffer_roads.plot(edgecolor='lightgrey',
                                       linewidth=1)
gpd.GeoSeries([taxi_station_500_buffer]).plot(ax=ax, 
                                              facecolor='none',
                                              edgecolor='grey',
                                              linestyle='--')

ax.axis('off')

可以看出我们需要的道路网都被正确裁切出来。

• 与叠加分析进行对比

需要注意的是，clip()中的mask参数，即蒙版矢量，无论是GeoDataFrame还是GeoSeries亦或是纯粹的shapely矢量，在执行裁切时，都会被整合为一个矢量对象整体，因此与之前文章介绍过的overlay()叠加分析有着本质上的不同。

举个实际的例子，当我们想算出整个柏林被出租车站点500米缓冲区所覆盖的步道路网总长度时，可以在上文裁切计算结果的基础上直接求得：

# 求得整个柏林被出租车站点500米缓冲区所覆盖的步道路网总长度
taxi_station_500buffer_roads.length.sum()

但当我们想要针对每个站点求出各自500米缓冲区内部的步道路网长度时，就需要叠加分析，因为叠加分析的矢量叠置操作是在df1与df2各自行元素两两之间建立起的：

# 整合所有出租车站点500米缓冲区GeoSeries
taxi_station_500_buffer_ = Berlin_transport.query("fclass=='taxi'") \
                                           .to_crs('EPSG:32633')

taxi_station_500_buffer_['geometry'] = taxi_station_500_buffer_.buffer(500)

# 对柏林步道与出租车缓冲区进行叠加分析
taxi_station_500buffer_roads_ = gpd.overlay(Berlin_footway,
                                           taxi_station_500_buffer_)

# 求出每个站点周边500米半径内步道路网长度
taxi_station_500buffer_roads_.groupby('osm_id_2') \
                             .agg({'geometry': lambda s: s.length.sum()})

查看裁切与叠加分析分别结果表路网矢量总长度也可以看出叠加分析中的结果是针对每个站点分别计算的，因此对于彼此重叠的站点500米缓冲区就会出现重复重叠的路段：

# 叠加分析路网结果表矢量总长度
taxi_station_500buffer_roads_.groupby('osm_id_2') \
                             .agg({'geometry': lambda s: s.length.sum()}) \
                             .geometry.sum()

# 裁切分析路网结果表矢量总长度
taxi_station_500buffer_roads.length.sum()