前言

將近兩年前，我寫過一篇同名文章（見使用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 則取出兩組資料相交的部分，結果如下圖所示：

繪圖如下：

可以看到確實取出了相交的部分，至此我們就得到了想要的結果。

結束

只要是需要判斷兩組空間物件空間位置的均可以使用此函式，其餘的諸如並集、取異等可以自行試驗，或參考官方文件。解決問題的途徑有很多，而最簡單最優美的解決方式總是無止境的，在解決某一實際問題時我們無需過多的思考如何最佳，但是當閒暇時刻靜下心來的時候還是應該想想碰到的問題如何解決才是最優的。