python基于geopandas的空间数据分析之二：空间计算_python 空间计算

1 前言

 在实际的空间数据分析过程中，数据可视化只是对最终分析结果的发布与展示，在此之前，根据实际任务的不同，需要衔接很多较为进阶的空间操作，本文就将对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专业的读者朋友解释一下什么是空间意义上的缓冲区，缓冲区用于表示点、线、面等矢量数据的影响范围或服务范围，思想很简单，即为矢量数据拓展出一定宽度的边，下图展示了点、线以及面分别对应的缓冲区的示意：

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

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

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

 可以看出，resolution参数对最终形成的缓冲区形态影响较大，但默认16的参数下已经可以较准确地逼近圆形，且缓冲距离还可以设置为负数，即几何对象向内收缩：

# 分别绘制多边形、多边形正向缓冲区、多边形负向缓冲区
gs = gpd.GeoSeries([polygon,
                    polygon.buffer(distance=1),
                    polygon.buffer(distance=-0.25)])      
ax = gs.plot(alpha=0.2)
ax.axis('off')
plt.savefig('图3.png', dpi=300, bbox_inches='tight', pad_inches=0)
1
2
3
4
5
6
7

 在本系列文章第一篇中介绍过shapely对矢量数据格式的合法性有一定规定，如多边形不能自交叉，可以通过查看要素的is_valid属性判断几何对象是否合法，而buffer()有一个隐藏功能就是其针对存在自交叉的非法面要素，对shapely.validation.make_valid()修复后产生的多部件面要素，创建一个非常小的缓冲区，即可在不损失多少数据精度的情况下，实现修复：

• 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')
1
2
3
4
5
6
7
8
9
10
11

• simplify

 当原始的矢量数据因为形状复杂，包含的点较多时，会导致其文件体积较大，如果我们需要在在线地图上叠加它们，太大体积的矢量数据不仅会拖慢网络传输速度，也会给图形的渲染带来更大的压力，这时对矢量数据进行简化就非常有必要，geopandas中沿用shapely中的simplify()方法，帮助我们对过于复杂的线和面进行简化，和QGIS中简化矢量的方法一样，simplify()使用了科学的Douglas-Peucker算法，基于预先设定的阈值ϵ，在递归判断的过程中删掉所有小于ϵ的点，其过程示意如图：

 譬如我们这里基于-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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

 可以看到在预设的阈值下，对应simplify()中的参数tolerance=0.5，折线得到有效地简化，这在搭建web GIS平台要渲染矢量数据时非常有效，有效简化后的矢量数据可以在不损失太多视觉感知到的准确度的同时，带来巨大的性能提升。

• unary_union

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

 并且如果原始数据中存在互相存在重叠的矢量对象，通过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
3
4
5
6

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('图11.png', dpi=300, bbox_inches='tight', pad_inches=0)
1
2
3
4

• 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)
1
2
3
4

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

2.3 叠加分析

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

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)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16

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

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

overlay_result = gpd.overlay(df1=polygon1,
                             df2=polygon2,
                             how='union')
overlay_result
1
2
3
4

 可以发现，有些行存在缺失值而有些行又是完整的，我们分别绘制出这两类记录行：

# 存在缺失值的行
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('图17.png', dpi=300, bbox_inches='tight', pad_inches=0)
1
2
3
4
5
6
7
8
9
10
11

 在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
1
2
3
4

 这时返回的结果中不再带有缺失值，因为intersection只保留df1和df2彼此相交的部分：

ax = overlay_result.plot()

ax.set_xlim((-1, 10))
ax.set_ylim((0, 10))
plt.savefig('图19.png', dpi=300, bbox_inches='tight', pad_inches=0)
1
2
3
4
5

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

overlay_result = gpd.overlay(df1=polygon1,
                             df2=polygon2,
                             how='difference')
overlay_result
1
2
3
4

 这时返回的结果中不再有value2字段，结合图13可以知晓在how='difference’下的返回结果与Arcgis中的擦除功能一样，返回的是df1中不与df2相交的部分，且以Multi的形式保留被切割开来的碎片矢量：

ax = overlay_result.plot()

ax.set_xlim((-1, 10))
ax.set_ylim((0, 10))
plt.savefig('图21.png', dpi=300, bbox_inches='tight', pad_inches=0)
1
2
3
4
5

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

overlay_result = gpd.overlay(df1=polygon1,
                             df2=polygon2,
                             how='symmetric_difference')
overlay_result
1
2
3
4

ax = overlay_result.plot()

ax.set_xlim((-1, 10))
ax.set_ylim((0, 10))
plt.savefig('图23.png', dpi=300, bbox_inches='tight', pad_inches=0)
1
2
3
4
5

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

overlay_result = gpd.overlay(df1=polygon1,
                             df2=polygon2,
                             how='identity')
overlay_result
1
2
3
4

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

# 存在缺失值的行
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('图25.png', dpi=300, bbox_inches='tight', pad_inches=0)
1
2
3
4
5
6
7
8
9
10
11

 从图中可以看出，在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({
	'value1':[1],
	'geometry':[shapely.geometry.LineString([(1,0.5),(-1,0.5),(-1,1),(0,1)])]
})

ax = polygon.plot(alpha=0.2)
ax = line.plot(ax=ax, alpha=0.2)
1
2
3
4
5
6
7
8
9
10
11
12
13

 True和False下结果如图所示：

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

 在实际工作中，可以根据具体需要来选择使用对应的参数组合来进行叠加分析。

2.4 空间融合与拆分

 有时候我们希望对矢量数据按照某些字段进行分组，再分别对非矢量列与矢量列进行聚合及合并，类似于pandas中的groupby.agg()；而有些时候我们希望把矢量类型为Multi-xxx的记录行拆分成独立要素行，譬如将Multi-Polygon拆分成Polygon组成的若干行。通过geopandass中的dissolve()和explode()方法，我们就可以实现上述功能：

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()
1
2
3

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

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

• dissolve()

 dissolve()用于对矢量数据进行融合，可以理解为对矢量数据进行groupby+agg操作，即指定的单个或多个字段值相等的分到一组，对非矢量字段进行指定规则的聚合计算，对矢量列进行融合，其主要参数如下：

by：用于指定分组所依据的字段，单个字段传入列名字符串，多个字段传入列名列表
aggfunc：对分组字段外的其他非矢量列采取的聚合方式，与pandas中的agg一致，默认为first，也可以像agg那样传入字段和函数一一对应的字典来分别聚合不同的列
as_index：bool型，用于设定是否在返回的结果中将分组依据列作为索引，默认为True

 我们以world数据集为例，为了方便演示我们首先新增字段less_than_median_gdp，用于判断对应的国家GDP是否小于世界中位数水平：

• explode()

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

gdf = gdp.GeoDataFrame({
	'value1':[1],
	'value2':['a'],
	'geometry':[shapely.geometry.MultiPoint([(0,0),(1,1),(2,3)])]
})
gdf
1
2
3
4
5
6

gdf.explode(index_parts=True)
1

2.5 空间连接

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

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

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

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

left_df：GeoDataFrame，传入空间连接对应的左表
right_df：GeoDataFrame，传入空间连接对应的右表
how：字符型，用于决定连接方式，'inner’表示内连接，且连接结果表中的矢量列来自左表；'left’表示左连接，且结果表中的矢量列来自左表；'right’表示右连接，最终结果表中的矢量列来自右表
op：（该参数在0.10版本以上修改为predicate）字符型，用于设定拓扑判断的规则，‘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中每个面文件对应的行政区划名称
1
2
3

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

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

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

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

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

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

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

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

• sjoin_nearest()

 有些时候我们需要判断的并不是左右两表中矢量列相交、包含等直接的拓扑关系，而是左右两表矢量列之间距离至多xx米内，彼此距离最近的成对匹配结果，像这样的空间距离关系判断，这在旧版本的geopandas中，通常可以左右两边分别做缓冲区后进行常规空间连接，再分组进行距离计算，最后才能筛选出所需的结果，颇为麻烦。

 0.10版以上的geopandas新增的sjoin_nearest()就可以支持我们一步到位开展上述分析计算功能，它的主要参数有：

left_df：连接对应的左GeoDataFrame
right_df：连接对应的右GeoDataFrame
how：设置连接方式，可选的有’left’、‘right’及’inner’，默认为’inner’
max_distance：重要参数，用于设置最大搜索距离阈值，当矢量间的距离小于此阈值时才会进行连接
lsuffix：设置左表重名字段后缀文字，默认为’left’
rsuffix：设置右表重名字段后缀文字，默认为’right’
distance_col：设置连接结果表中记录对应矢量间距离的字段名称，默认不设置时不会在结果表中添加距离信息

 下面我们来通过一个简单的例子来体验这个功能：

import geopandas as gpd
from shapely.geometry import Point

# 构造示例点要素表1
gdf1 = gpd.GeoDataFrame(
    {
        'id1': list('abc'),
        'geometry': [
            Point(0, 0),
            Point(1, 0),
            Point(-1, 0)
        ]
    }
)

# 构造示例点要素表2
gdf2 = gpd.GeoDataFrame(
    {
        'id2': list('def'),
        'geometry': [
            Point(0.4, 0),
            Point(1.2, 0),
            Point(-1.3, 0)
        ]
    }
)

ax = gdf1.plot(color='red')
ax = gdf2.plot(color='green', ax=ax)
ax.axis('equal')
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30

 颜色即用来区分我们的左右表对应矢量点位置，下面直接运用sjoin_nearest()进行空间最近连接，设置的距离阈值为0.35：

gpd.sjoin_nearest(gdf1, gdf2, max_distance=0.35, distance_col='对应距离')
1

2.6 拓扑关系判断

 geopandas中除了 叠加分析 以及 空间连接中基于拓扑关系判断实现多表数据联动 之外，还针对GeoSeries与GeoDataFrame设计了一系列方法，可以直接进行矢量数据之间的拓扑关系判断并返回对应的bool型判断结果，下面罗列出常用的一些拓扑关系判断API，均为GeoSeries或GeoDataFrame的方法：

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

 以contains()为例，在比较矢量数据之间拓扑关系时，矢量数据与待比较矢量数据之间主要有以下几种格式：

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

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

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

• 长度m-n（n>0）与长度n进行比较
  这时的情况就与前面一种类似，即从头开始两两位置匹配上的要素才会进行比较及结果的输出，多出的得不到匹配的要素会自动返回False：
  

 geopandas中进行拓扑关系判断的基本原则了解完了。

2.7 空间裁切

 在空间数据分析中，裁切也是非常常用的操作，譬如我们想要获取某个公交站周围500米半径内部的路网矢量，就可以使用到裁切。在geopandas中我们可以使用clip()函数来基于蒙版矢量对目标矢量进行裁切，其主要参数如下：

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

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

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

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

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

 在交互模式下同时绘制出缓冲区以及裁切出的路网：
 可以看出我们需要的道路网都被正确裁切出来。

• 与叠加分析进行对比

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

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

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

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

2.8 新版本（0.10以上）增加的特征

• sjoin()、sjoin_nearest()、overlay()和clip()亦可作为GeoDataFrame的方法来使用
  在以前的版本中，我们只能使用gpd.XXX()的方式来使用sjoin()、overlay()、clip()等方法，而在这次新版本更新中，我们可以像pandas里的merge()、join()那样作为方法使用，以上文介绍的sjoin_nearest()为例，只需向sjoin_nearest()方法中传入右表即可：

gdf1.sjoin_nearest(gdf2, 
                   max_distance=0.35, 
                   distance_col='对应距离')
1
2
3

• GeoSeries新增批量XY转点方法from_xy()
  新版本中为GeoSeries对象新增了from_xy()方法来快速实现坐标转点，下面与gpd.points_from_xy()的效果进行对比：

gpd.points_from_xy(x=range(10), y=range(10))

gpd.GeoSeries.from_xy(x=range(10), y=range(10))
1
2
3

• 支持对矢量数据自动推断合适的横轴墨卡托坐标参考系
  其实这个特性在0.9版本中就已加入，但是还有一些小问题，而新版本中这个功能更加完善，效果如下：
  

• to_file()方法在driver参数缺省时可自动识别导出文件类型
  在新版本中，若未在to_file()中指定driver参数，geopandas会自动根据文件后缀名来自动推断要导出的矢量文件类型：
  

参考文章：
1.（数据科学学习手札84）基于geopandas的空间数据分析——空间计算篇（上）
2.（数据科学学习手札88）基于geopandas的空间数据分析——空间计算篇（下）
3.（数据科学学习手札129）geopandas 0.10版本重要新特性一览