（資料科學學習手札84）基於geopandas的空間資料分析——空間計算篇（上）

費弗裡發表於2020-05-22

本文示例程式碼已上傳至我的Github倉庫https://github.com/CNFeffery/DataScienceStudyNotes

1 簡介

　　在本系列之前的文章中我們主要討論了geopandas及其相關庫在資料視覺化方面的應用，各個案例涉及的資料預處理過程也僅僅涉及到基礎的向量資料處理。在實際的空間資料分析過程中，資料視覺化只是對最終分析結果的釋出與展示，在此之前，根據實際任務的不同，需要銜接很多較為進階的空間操作，本文就將對geopandas中的部分空間計算進行介紹。

　　本文是基於geopandas的空間資料分析系列文章的第8篇，通過本文你將學習到geopandas中的空間計算（由於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專業的讀者朋友解釋一下什麼是空間意義上的緩衝區，緩衝區用於表示點、線、面等向量資料的影響範圍或服務範圍，思想很簡單，即為向量資料擴充出一定寬度的邊，圖1展示了點、線以及面分別對應的緩衝區的示意：

（資料科學學習手札84）基於geopandas的空間資料分析——空間計算篇（上）

圖1

　　而建立緩衝區時也需要遵循一定的引數，從而決定怎樣向幾何物件外進行緩衝，geopandas中buffer()和shapely中的buffer()方法引數一致，主要引數如下：

distance：用於指定向外緩衝的距離，單位與向量資料自帶單位保持一致，在常見的投影座標系如Web Mercator（EPSG:3857）下就是以米為單位，因此需要注意一定要先將向量資料轉換為合適的投影座標系之後，再進行緩衝區分析才是合理有效的

resolution：因為在建立緩衝區時，對於構成向量物件的每一個點，都會以對應點為中心向外建立半徑=緩衝區距離的圓，而Polygon型別始終是由有限個點所構成的，因此需要近似拼接出圓形的輪廓，resolution引數就用於決定每個四分之一圓弧上使用多少段連續的線段來近似拼接以表示圓的形狀，預設引數值為16，足以近似模擬圓面積的99.8%

　　下面我們分別對點、線以及面繪製不同resolution引數取值下緩衝前後的對比圖：

圖2

　　可以看出，resolution引數對最終形成的緩衝區形態影響較大，但預設16的引數下已經可以較準確地逼近圓形，且緩衝距離還可以設定為負數，即幾何物件向內收縮：

# 分別繪製多邊形、多邊形正向緩衝區、多邊形負向緩衝區
ax = gpd.GeoSeries([polygon,
                    polygon.buffer(distance=1),
                    polygon.buffer(distance=-0.25)]) \
        .plot(alpha=0.2)
ax.axis('off')

plt.savefig('圖3.png', dpi=300, bbox_inches='tight', pad_inches=0)

圖3

　　在本系列文章第一篇中介紹過shapely對向量資料格式的合法性有一定規定，如多邊形不能自交叉，可以通過is_valid()方法判斷幾何物件是否合法，而buffer()有一個隱藏功能就是其可以通過對非法的幾何物件建立距離為0的緩衝區來修正構成向量物件的點的不合理連線順序，從而使得向量物件變為合法的：

圖4

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')

圖5

simplify()

　　當原始的向量資料因為形狀複雜，包含的點較多時，會導致其檔案體積較大，如果我們需要在線上地圖上疊加它們，太大體積的向量資料不僅會拖慢網路傳輸速度，也會給圖形的渲染帶來更大的壓力，這時對向量資料進行簡化就非常有必要，geopandas中沿用shapely中的simplify()方法，幫助我們對過於複雜的線和麵進行簡化，和QGIS中簡化向量的方法一樣，simplify()使用了科學的Douglas-Peucker演算法，基於預先設定的閾值\(\epsilon\)，在遞迴判斷的過程中刪掉所有小於\(\epsilon\)的點，其過程示意如圖6：

圖6

　　譬如我們這裡基於-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)

圖7

　　可以看到在預設的閾值下，對應simplify()中的引數tolerance=0.5，折線得到有效地簡化，這在搭建web GIS平臺要渲染向量資料時非常有效，有效簡化後的向量資料可以在不損失太多視覺感知到的準確度的同時，帶來巨大的效能提升。

unary_union

　　我們都知道，不管是GeoSeries還是GeoDataFrame，其每一行資料都代表獨立的shapely向量要素，而通過unary_union屬性，我們可以將一整列向量合併為單獨的一個shapely向量物件，從而方便我們進行一些其他的操作：

圖8

　　並且如果原始資料中存在互相存在重疊的向量物件，通過unary_union之後，返回的shapely物件會自動對存在重疊的向量物件進行融合，這一點可以方便我們的很多日常操作：

圖9

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('圖10.png', dpi=300, bbox_inches='tight', pad_inches=0)

圖10

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)

圖11

translate()

　　translate()用於實現向量物件的平移操作，其主要引數有xoff和yoff，分別控制在x維度和y維度上的平移距離（與對應的投影單位保持一致）:

圖12

2.3 疊加分析

　　geopandas基於shapely中的overlay()，為GeoDataFrame賦予了同樣的可以作用到整個向量列的overlay()，使得我們可以對兩個GeoDataFrame中全部的向量物件兩兩之間進行基於集合關係的疊加分析（如圖13）:

圖13

　　overlay()中的主要引數如下：

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)

圖14

　　接下來我們將其中紅色部分對應的GeoDataFrame作為df1，灰色部分作為df2，來比較overlay()中不同引數對應的效果：

how='union'

　　首先我們設定how='union'，對polygon1與polygon2進行疊加分析：

overlay_result = gpd.overlay(df1=polygon1,
                           df2=polygon2,
                           how='union')
overlay_result

　　得到的結果如圖15：

圖15

　　可以發現，有些行存在缺失值而有些行又是完整的，我們分別繪製出這兩類記錄行：

# 存在缺失值的行
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('圖16.png', dpi=300, bbox_inches='tight', pad_inches=0)

圖16

　　在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

圖17

　　這時返回的結果中不再帶有缺失值，因為intersection只保留df1和df2彼此相交的部分：

ax = overlay_result.plot()

ax.set_xlim((-1, 10))
ax.set_ylim((0, 10))
plt.savefig('圖18.png', dpi=300, bbox_inches='tight', pad_inches=0)

圖18

how='difference'

　　在how='difference'引數下對polygon1和polygon2進行疊加分析：

overlay_result = gpd.overlay(df1=polygon1,
                             df2=polygon2,
                             how='difference')
overlay_result

圖19

　　這時返回的結果中不再有value2欄位，結合圖13可以知曉在how='difference'下的返回結果與Arcgis中的擦除功能一樣，返回的是df1中不與df2相交的部分，且以Multi的形式保留被切割開來的碎片向量：

ax = overlay_result.plot()

ax.set_xlim((-1, 10))
ax.set_ylim((0, 10))
plt.savefig('圖20.png', dpi=300, bbox_inches='tight', pad_inches=0)

圖20

how='symmetric_difference'

　　'symmetric的意思是對稱的，而在how='symmetric_difference'條件下，與Arcgis中的交集取反功能相同，兩個df中不相交的部分都會被返回：

overlay_result = gpd.overlay(df1=polygon1,
                             df2=polygon2,
                             how='symmetric_difference')
overlay_result

圖21

圖22

how='identity'

　　最後我們來看看how='identity'，對應Arcgis中的標識功能：

overlay_result = gpd.overlay(df1=polygon1,
                             df2=polygon2,
                             how='identity')
overlay_result

圖23

　　為了方便理解，我們同樣分別繪製出存在缺失值的行以及不存在缺失值的行：

# 存在缺失值的行
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('圖24.png', dpi=300, bbox_inches='tight', pad_inches=0)

圖24

　　從圖24中可以看出，在how='identity'條件下，所有df1中不與df2相交的部分，以及兩者相交的部分作為返回結果，且每個相交的部分都變為單獨的要素帶上所有涉及的屬性欄位，而df1中不涉及相交的部分則仍然以Multi的形式被返回。

keep_geom_type

　　有些時候我們需要做的不僅僅是面與面之間的疊加分析。比如在計算路網相關的指標時，我們可能會需要與目標區域存在疊置關係的部分路網，這就存在面與線之間的疊加分析。引數keep_geom_type就用於設定最終返回的向量資料型別是否必須與df1對應的型別相同，下面我們構造示例資料來學習keep_geom_type引數的作用：

圖25

　　True和False下結果如圖26所示：

圖26

　　其中GeometryCollection型別代表多型別要素集合，比如這裡疊加分析的結果包含了一條線和一個點：

圖27

　　在實際工作中，可以根據具體需要來選擇使用對應的引數組合來進行疊加分析。

2.4 空間融合與拆分

　　有時候我們希望對向量資料按照某些欄位進行分組，再分別對非向量列與向量列進行聚合及合併，類似於pandas中的groupby.agg()；而有些時候我們希望把向量型別為Multi-xxx的記錄行拆分成獨立要素行，譬如將Multi-Polygon拆分成Polygon組成的若干行。通過geopandass中的dissolve()和explode()方法，我們就可以實現上述功能：

dissolve()

　　dissolve()用於對向量資料進行融合，可以理解為對向量資料進行groupby+agg操作，即指定的單個或多個欄位值相等的分到一組，對非向量欄位進行指定規則的聚合計算，對向量列進行融合，其主要引數如下：

by：用於指定分組所依據的欄位，單個欄位傳入列名字串，多個欄位傳入列名列表

aggfunc：對分組欄位外的其他非向量列採取的聚合方式，與pandas中的agg一致，預設為first，也可以像agg那樣傳入欄位和函式一一對應的字典來分別聚合不同的列

as_index：bool型，用於設定是否在返回的結果中將分組依據列作為索引，預設為True

　　我們以world資料集為例，為了方便演示我們首先新增欄位less_than_median_gdp，用於判斷對應的國家GDP是否小於世界中位數水平：