使用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 則取出兩組數據相交的部分，結果以下圖所示：

繪圖以下：

能夠看到確實取出了相交的部分，至此咱們就獲得了想要的結果。

結束

只要是須要判斷兩組空間對象空間位置的都可以使用此函數，其他的諸如並集、取異等能夠自行試驗，或參考官方文檔。解決問題的途徑有不少，而最簡單最優美的解決方式老是無止境的，在解決某一實際問題時咱們無需過多的思考如何最佳，可是當閒暇時刻靜下心來的時候仍是應該想一想碰到的問題如何解決纔是最優的。