> 本文示例程式碼已上傳至我的`Github`倉庫[https://github.com/CNFeffery/DataScienceStudyNotes](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展示了點、線以及面分別對應的緩衝區的示意：

圖1 　　而建立緩衝區時也需要遵循一定的引數，從而決定怎樣向幾何物件外進行緩衝，`geopandas`中`buffer()`和`shapely`中的`buffer()`方法引數一致，主要引數如下： > **distance**：用於指定向外緩衝的距離，單位與向量資料自帶單位保持一致，在常見的投影座標系如*Web Mercator*（EPSG:3857）下就是以米為單位，因此需要注意一定要先將向量資料轉換為合適的投影座標系之後，再進行緩衝區分析才是合理有效的 > > **resolution**：因為在建立緩衝區時，對於構成向量物件的每一個點，都會以對應點為中心向外建立半徑=緩衝區距離的圓，而Polygon型別始終是由有限個點所構成的，因此需要近似拼接出圓形的輪廓，*resolution*引數就用於決定每個四分之一圓弧上使用**多少段**連續的線段來近似拼接以表示圓的形狀，預設引數值為16，足以近似模擬圓面積的99.8% 　　下面我們分別對點、線以及面繪製不同**resolution**引數取值下緩衝前後的對比圖：

圖2 　　可以看出，**resolution**引數對最終形成的緩衝區形態影響較大，但預設16的引數下已經可以較準確地逼近圓形，且緩衝距離還可以設定為負數，即幾何物件向內收縮： ```Python # 分別繪製多邊形、多邊形正向緩衝區、多邊形負向緩衝區 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陣列： ```Python 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()`來簡化它： ```Python 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度之後的美國： ```Python 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倍放縮，紅色代表縮放之後： ```Python 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()`不同引數下結果的區別： ```Python 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`進行疊加分析： ```Python overlay_result = gpd.overlay(df1=polygon1, df2=polygon2, how='union') overlay_result ``` 　　得到的結果如圖15： 圖15 　　可以發現，有些行存在缺失值而有些行又是完整的，我們分別繪製出這兩類記錄行： ```Python # 存在缺失值的行 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`進行疊加分析： ```Python overlay_result = gpd.overlay(df1=polygon1, df2=polygon2, how='intersection') overlay_result ``` 圖17 　　這時返回的結果中不再帶有缺失值，因為`intersection`只保留**df1**和**df2**彼此相交的部分： ```Python 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`進行疊加分析： ```Python overlay_result = gpd.overlay(df1=polygon1, df2=polygon2, how='difference') overlay_result ``` 圖19 　　這時返回的結果中不再有`value2`欄位，結合圖13可以知曉在`how='difference'`下的返回結果與`Arcgis`中的**擦除**功能一樣，返回的是**df1**中不與**df2**相交的部分，且以`Multi`的形式保留被切割開來的碎片向量： ```Python 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中不相交的部分都會被返回： ```Python overlay_result = gpd.overlay(df1=polygon1, df2=polygon2, how='symmetric_difference') overlay_result ``` 圖21 圖22 - **how='identity'** 　　最後我們來看看`how='identity'`，對應`Arcgis`中的**標識**功能： ```Python overlay_result = gpd.overlay(df1=polygon1, df2=polygon2, how='identity') overlay_result ``` 圖23 　　為了方便理解，我們同樣分別繪製出存在缺失值的行以及不存在缺失值的行： ```Python # 存在缺失值的行 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是否小於世界中位數水平： 圖28 　　接著我們以國家對應大洲列`continent`為分組依據，並對人口和GDP列進行求和，如圖29所示，在非向量列得到對應的聚合計算之後，向量列也被融合為`Multi-Polygon`： 圖29 - **explode()** 　　`explode()`功能與`dissolve()`相反，用於將`Multi-xxx`或`Geometry-Collection`型別的資料從一行拆分到多行，如下面的例子，非向量欄位會自動填充到每一行： 圖30 　　以上就是本文的全部內容，關於更多`geopandas`中**空間計算**的內容，我們將在下一篇中繼續