(數據科學學習手札88)基於geopandas的空間數據分析——空間計算篇(下)


本文示例代碼及數據已上傳至我的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中除了在上一篇文章中介紹的疊加分析以及上文介紹的空間連接中基於拓撲關系判斷實現多表數據聯動之外,還針對GeoSeriesGeoDataFrame設計了一系列方法,可以直接進行矢量數據之間的拓撲關系判斷並返回對應的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,均為GeoSeriesGeoDataFrame的方法:

intersects():檢查相交關系

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

within():檢查主體矢量是否在待檢查矢量的內部

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

crosses():檢查交叉關系,常見如線與線之間的交叉

disjoint():檢查不相交關系,即兩個矢量之間沒有任何接觸

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

overlaps():檢查重疊關系

2.3 空間裁切

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

  在geopandas中我們可以使用clip()函數來基於蒙版矢量對目標矢量進行裁切,其主要參數如下:

gdfGeoDataFrameGeoSeries,代表將要被裁切的矢量數據集

maskGeoDataFrameGeoSeriesshapely中的PolygonMulti-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相關的新特性或額外知識,依舊會不定期作為系列文章的補充,總結發布出來與大家分享。


  與熱愛的技術一起成長🧐。


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM