前言

将近两年前，我写过一篇同名文章（见使用Python实现子区域数据分类统计）。

当时是为了统计县域内的植被覆盖量，折腾了一段时间，解决了这个问题。最近，又碰到了一个类似的需求，也需要统计某个小范围内的数据。简单来说，这个需求是将两个 shp 文件的任意两个对象做相交判断，最后形成一个新的空间对象集合，最后对此集合进行简单统计分析即可。

解决方案

明白了这一点之后，再看之前的代码，就发现当时用了很笨的方法。写了两个循环，先是取出大范围的 shp 中的每一个对象，再读取小范围 shp 的每一个对象，将小范围的 shp 空间对象逐个与大的空间对象进行相交操作。循环操作，就能将大范围的 shp 对象完全切割完毕。这段话说的稍微有些绕，感兴趣的可以直接看一下上篇文章。

今天又一次碰到了这个问题，回头找到了原来的文章，但是总感觉写的很丑，难道必须要用这么难看的方法来解决这个问题吗？想了半天，有没有简单的方法能够解决呢？

思考半天，找到了答案，直接对两个 GeoDataFrame 对象做类似数据库的 join 操作不就可以了嘛，只是任意两个判断的时候用空间操作代替数据库的匹配操作。想到这，就开始翻看 geopandas 的用户手册，果然让我找到了。

解决路径

1. 包引用

from geopandas import *
from shapely.geometry import * 

2. 创建两个 GeoDataFrame 对象

geopandas 可以直接将 shp 文件读为 GeoDataFrame 对象，如下：

shpdata = GeoDataFrame.from_file(path)

此处，采用模拟的方式创建两个 GeoDataFrame 对象，如下：

p1 = Point([1, 2])
p2 = Point([1.5, 1.7])
p3 = Point([1.8, 1.5])
p4 = Point([1.4, 2.2])
gdf1 = geopandas .GeoSeries([p1, p2]).buffer(0.3)
gdf2 = geopandas .GeoSeries([p3, p4]).buffer(0.2)

首先创建4个点对象，使用前两个创建第一个 GeoSeries 对象，后两个创建第二个 GeoSeries 对象。 buffer 函数执行缓冲区分析，将点以一定的距离扩展成面。GeoSeries 简单的说是只包含空间属性的对象，不包含 GeoDataFrame 的其他字段，所以需要为其附加其他字段，为第一个添加 left 字段，为第二个添加 right 字段，并赋值，如下：

gdf1 = GeoDataFrame({"left": [1, 2]}, geometry=gdf1)
gdf2 = GeoDataFrame({"right":[3, 4]}, geometry=gdf2)

两个 GeoDataFrame 如下：

gdf1 gdf2

通过画图可以看出两个对象的位置关系：

ax = gdf1.plot(color='red')
gdf2.plot(ax=ax, color='green')

3. 两两相交

官网翻阅半天，找到了 overlay 函数，overlay 是覆盖的意思，从意思我们就能猜测出是对两个对象做覆盖的操作。

官网介绍如下：

When working with multiple spatial datasets – especially multiple polygon or line datasets – users often wish to create new shapes based on places where those datasets overlap (or don’t overlap). These manipulations are often referred using the language of sets – intersections, unions, and differences. These types of operations are made available in the geopandas library through the overlay function.

The basic idea is demonstrated by the graphic below but keep in mind that overlays operate at the DataFrame level, not on individual geometries, and the properties from both are retained. In effect, for every shape in the first GeoDataFrame, this operation is executed against every other shape in the other GeoDataFrame:

参考http://geopandas.org/set_operations.html

大意是说当执行两个空间对象的相交、合并、取异操作的时候就可以使用此函数。

此函数可以判断两个空间对象的交集、并集以及不同的部分，此处我们只需要取出交集就可以了。

intersection_data = geopandas.overlay(gdf1, gdf2, how='intersection')

参数 how 设置为 intersection 则取出两组数据相交的部分，结果如下图所示：

绘图如下：

可以看到确实取出了相交的部分，至此我们就得到了想要的结果。

结束

只要是需要判断两组空间对象空间位置的均可以使用此函数，其余的诸如并集、取异等可以自行试验，或参考官方文档。解决问题的途径有很多，而最简单最优美的解决方式总是无止境的，在解决某一实际问题时我们无需过多的思考如何最佳，但是当闲暇时刻静下心来的时候还是应该想想碰到的问题如何解决才是最优的。