使用Python实现子区域数据分类统计

前言

将近两年前，我写过一篇同名文章（见使用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 以下：

经过画图能够看出两个对象的位置关系：

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 则取出两组数据相交的部分，结果以下图所示：

绘图以下：

能够看到确实取出了相交的部分，至此咱们就获得了想要的结果。

结束

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