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

1 簡介

　　在基於geopandas的空間資料分析系列文章第8篇中，我們對geopandas開展空間計算的部分內容進行了介紹，涉及到緩衝區分析、向量資料簡化、仿射變換、疊加分析與空間融合等常見空間計算操作，而本文就將針對geopandas中剩餘的其他常用空間計算操作進行介紹。

　　本文是基於geopandas的空間資料分析系列文章的第9篇，也是整個系列文章主線部分內容的最後一篇，通過本文，你將學習到geopandas中的更多常用空間計算方法。

2 基於geopandas的空間計算

　　承接上文內容，geopandas中封裝的空間計算方法除了系列上一篇文章中介紹的那幾種外，還有其他的幾類，下面我們繼續來學習：

2.1 空間連線

　　類比常規表格資料的連線操作，在空間資料分析中也存在類似表連線的操作，譬如我們手頭有一張包含設施點資料的矢量表，以及另一張包含行政區劃面資料的矢量表，當我們想要通過某些操作來統計出每個行政區劃面內部的設施點資訊時，空間連線就可以非常方便快捷地實現這類需求。

　　我們都清楚常規表格資料的連線，是按照設定的連線方式，將每張表中指定的某列或某些列數值相等的記錄行合併為同一行，最後彙整成連線結果表返回：

圖1

　　而空間連線不同於常規表連線，其合併同一行的依據不是檢查指定的列數值是否相等，而是基於不同矢量表其向量列

圖2

　　在geopandas中我們利用sjoin函式來實現空間連線，其使用方式類似pandas中的merge接近，主要引數如下：

left_df：GeoDataFrame，傳入空間連線對應的左表

right_df：GeoDataFrame，傳入空間連線對應的右表

how：字元型，用於決定連線方式，'inner'表示內連線，且連線結果表中的向量列來自左表；'left'表示左連線，且結果表中的向量列來自左表；'right'表示右連線，最終結果表中的向量列來自右表

op：字元型，用於設定拓撲判斷的規則，'intersects'代表相交，即幾何物件之間存在共有的邊或內部點；'contains'

代表包含，即一個幾何物件至少有一個點位於另一個幾何物件內部，且其本身沒有任何點落在另一個結幾何物件的外部；'within'表示在內部，是'contains'的相反情況，即左表被右表向量'contains'

lsuffix：字元型，代表當左右表連線之後存在重名列時，為左表重名的列新增的字尾，預設為'left'

rsuffix：字元型，意義類似lsuffix，預設為'right'

　　瞭解過sjoin()中的核心引數後，我們來通過實際例子理解它們的具體作用，how的作用與pandas中效果的一致，這裡不多解讀，我們來重點學習op各引數的不同效果：

• 引數op

　　intersects是空間連線中最常使用的模式，即相比較的兩個幾何物件有至少1個公共點就會被匹配上，下面我們以柏林公交站點資料為例，首先我們先讀入柏林行政區劃面資料，其中欄位Gemeinde_n是每個行政區劃的名稱：

# 讀入柏林行政區劃面檔案

Berlin = gpd.read_file('Berlin/Bezirke__Berlin.shp')

Berlin.head() # Gemeinde_n代表鎮，即Berlin中每個面檔案對應的行政區劃名稱

圖3

　　接著再讀入柏林全部交通車站資料，其中fclass列代表對應車站的類別：

Berlin_transport = gpd.read_file('Berlin/gis_osm_transport_free_1.shp')

Berlin_transport.head()

圖4

　　對站點的空間分佈進行視覺化：

圖5

　　接著我們就利用sjoin()將區劃面作為左表，站點作為右表，在op='intersects'引數設定下進行空間連線，再銜接groupby，以統計出各區劃面內部的公交站點數量：

gpd.sjoin(left_df=Berlin,

          right_df=Berlin_transport.query("fclass=='bus_stop'"),

          op='intersects') \

   .groupby('Gemeinde_n') \

   .size()

圖6

　　再設定op='contains'，因為進行連線的物件是左表面要素，右表點要素，所以這裡的效果等價於op='intersects'：

圖7

　　但當op='within'時，按照拓撲規則，如果依舊是左表面要素，右表點要素，得到的結果就會為空，反過來則正常：

圖8

　　類似的，其他型別幾何物件之間的空間連線你也可以根據自己的需要進行操作，值得一提的是，利用sjoin()進行空間左、右、內連線時，因為結果表依舊是GeoDataFrame，所以只會保留一列向量列，按照上文中引數介紹部分的描述，只有右連線時結果表中的向量列才來自右表，但無論採取什麼連線方式，結果表中未被保留的向量列對應的index會被作為單獨的一列儲存下來，幫助我們可以按圖索驥利用loc方式索引出需要的資料：

圖9

2.2 拓撲關係判斷

　　geopandas中除了在上一篇文章中介紹的疊加分析以及上文介紹的空間連線中基於拓撲關係判斷實現多表資料聯動之外，還針對GeoSeries與GeoDataFrame設計了一系列方法，可以直接進行向量資料之間的拓撲關係判斷並返回對應的bool型判斷結果，以contains()為例，在比較向量資料之間拓撲關係時，向量資料與待比較向量資料之間主要有以下幾種格式：

• 長度n與長度1進行比較

　　當主體向量列長度為n，而輸入待比較的向量列長度為1時，返回的bool值是待比較向量列與主題向量列一一進行比較後的結果：

圖10

• 長度1與長度n進行比較

　　與前面一種情況類似，只不過這裡是將主體向量列與待比較向量列一一比較之後的結果：

圖11

• 長度m與長度m-n（n>0）進行比較

　　這裡所說的情況指主體向量與待比較向量長度都不為1，且主體向量列的長度大於待比較向量，這時返回的結果只會對主體向量列前m-n個要素與待比較向量對應位置一一比較，主體向量被截斷未能進行比較的部分預設返回False：

圖12

• 長度m-n（n>0）與長度n進行比較

　　這時的情況就與前面一種類似，即從頭開始兩兩位置匹配上的要素才會進行比較及結果的輸出，多出的得不到匹配的要素會自動返回False：

圖13

　　geopandas中進行拓撲關係判斷的基本原則瞭解完了，下面羅列出常用的一些拓撲關係判斷API，均為GeoSeries或GeoDataFrame的方法：

intersects()：檢查相交關係

contains()：檢查包含關係，即主體向量完全包裹住待比較的向量且它們的邊界互不接觸，譬如面對點的包含

within()：檢查主體向量是否在待檢查向量的內部

touches()：檢查觸碰關係，即兩個向量之間至少有一個1個公共點，但它們的內部無任何相交區域

crosses()：檢查交叉關係，常見如線與線之間的交叉

disjoint()：檢查不相交關係，即兩個向量之間沒有任何接觸

geom_equals()：檢查是否完全相同

overlaps()：檢查重疊關係

2.3 空間裁切

　　在空間資料分析中，裁切也是非常常用的操作，譬如我們想要獲取某個公交站周圍500米半徑內部的路網向量，就可以使用到裁切。

　　在geopandas中我們可以使用clip()函式來基於蒙版向量對目標向量進行裁切，其主要引數如下：

gdf：GeoDataFrame或GeoSeries，代表將要被裁切的向量資料集

mask：GeoDataFrame、GeoSeries或shapely中的Polygon、Multi-Polygon物件，代表蒙版向量

keep_geom_type：同疊加分析overlay中的同名引數

　　基於實際例子進行演示，我們讀入資料berlin_footway_WGS84.shp，包含了柏林全部的步道路網線資料，並轉換到適合柏林地區的投影EPSG:32633：

圖14

　　接下來我們從上文中使用到的柏林車站點資料中篩選出租車站點，與步道路網資料統一座標參考系，生成500米緩衝區，並利用上一篇文章中介紹過的unary_union來得到MultiPolygon物件：

圖15

　　萬事俱備，接下來我們使用clip()來裁切所有計程車站點500米緩衝區內部的步行道路網：

# 裁切所有計程車站點500米緩衝區內部的路網線資料

taxi_station_500buffer_roads = gpd.clip(gdf=Berlin_footway,

                                        mask=taxi_station_500_buffer)

　　在互動模式下同時繪製出緩衝區以及裁切出的路網：

圖16

　　可以看出我們需要的道路網都被正確裁切出來。

• 與疊加分析進行對比

　　需要注意的是，clip()中的mask引數，即蒙版向量，無論是GeoDataFrame還是GeoSeries亦或是純粹的shapely向量，在執行裁切時，都會被整合為一個向量物件整體，因此與之前文章介紹過的overlay()疊加分析有著本質上的不同。

　　舉個實際的例子，當我們想算出整個柏林被計程車站點500米緩衝區所覆蓋的步道路網總長度時，可以在上文裁切計算結果的基礎上直接求得：

圖17

　　但當我們想要針對每個站點求出各自500米緩衝區內部的步道路網長度時，就需要疊加分析，因為疊加分析的向量疊置操作是在df1與df2各自行元素兩兩之間建立起的：

圖18

　　檢視裁切與疊加分析分別結果表路網向量總長度也可以看出疊加分析中的結果是針對每個站點分別計算的，因此對於彼此重疊的站點500米緩衝區就會出現重複重疊的路段：

圖19

3 寫在最後

　　從2020年2月8日釋出了geopandas空間資料分析系列第一篇文章，到今天這篇為止，geopandas中全部實用的主線內容（截至0.7.0版本），都在這斷斷續續撰寫完成的9篇文章中介紹完畢，不敢說是geopandas中文資料裡最好的，但穿插了眾多例子和舉一反三的內容，絕對是幫助大家理解學習geopandas非常實在的參考資料。

　　撰寫本系列文章的初衷，一是因為我對pandas的高度熟悉，二是由於喜歡程式設計，對ArcGIS之類主要靠點選相應按鈕完成任務且容易出錯的空間分析軟體不太喜歡，所以在瞭解到有這麼一個與pandas有著莫大淵源且可以做很多實用的空間計算操作的Python庫時，萌發出濃鬱的學習興趣，便將整個對geopandas相關內容學習精進的過程記錄下來，通過部落格與微信公眾號與廣大的讀者朋友共同交流學習，期間認識了很多業內大牛和朋友，收穫了很多很多。

　　geopandas是一個非常優秀的工具，它給了我們進行空間計算的多一種選擇，我目前所有工作中涉及到的可以用geopandas解決的問題，都會在jupyter中建立順滑的工作流。geopandas也是一個不斷髮展不斷迭代優化的開源專案，本系列主線內容雖已完結，但之後關於geopandas相關的新特性或額外知識，依舊會不定期作為系列文章的補充，總結髮布出來與大家分享。

　　與熱愛的技術一起成長