通過geopandas.sjoin()函數按多邊形范圍分割點

最近有一批點和多變型的數據，需要將點按照多邊形的區域進行分割。

經過若干嘗試，終於通過geopandas的sjoin函數得以實現。

這里首先感謝博主“張da統帥”的分享，使得本人獲得該實現方法的靈感，原帖見：https://beiyuan.me/geospatial-analysis-with-python-4/

以下是實現步驟及代碼：

 

1. 如果點文件為帶有坐標的文本文件，則先將其寫入GeoDataFrame類對象，若為shpfile文件，則直接用geopandas進行讀取。

# 載入geopandas庫
import geopandas as gpd
# 用geopandas讀取點shpfile文件
point = gpd.read_file(point_path)

 

2. 用geopandas讀取多邊形shp文件

# 讀取多邊形shpfile文件
polygon = gpd.read_file(polygon_path)

 

3. 通過sjoin()函數獲取點與多邊形相交的系列，組成新的GeoDataFrame類對象

new_gdf = gpd.sjoin(point, polygon, op='intersects')

這里new_gdf包含了所有point和polygon相交元素的信息。

new_gpd的geometry列值為點的坐標，兩個對象的其他列都合並在了new_gdf中（可以通過how關鍵字參數選擇只包含point/polygon的列值，具體見：http://geopandas.org/reference/geopandas.sjoin.html）

op關鍵字參數有 {‘intersects’, ‘contains’, ‘within’}三種可選，由於geopandas給出的文本鏈接不可用（可能是因為某種神秘力量吧），所以三種參數對應方法暫不明了。

 

4. 按照polygon的范圍對point進行分割

# new——gdf中的每一行代表一個點，每個點中都包含了其所在多邊形的所有列（geometry列除外）
# 由於index屬性是每個DataFrame對象都有的，因此利用index值判斷點屬於哪個多邊形
# polygon的index值在new_gdf中默認為index_rigth, 具體見：http://geopandas.org/reference/geopandas.sjoin.html

# 獲取每個點對應的polygon的index值
indexs = new_gdf.index_right.values
# 去掉s中重復的index值（由於有多個點在一個多邊形中的情況）
s = list(set(indexs))
# 從new_gdf中揀出在index值對應多邊形中的點，存入point_split列表（也可用point_split.to_file()方法直接寫入shpfile文件）
point_split = []
for index in indexs:
　　point_split.append(new_gdf.loc[gdf['index_right'] == index, ['想要保存的列名0'， ‘想要保存的列名1’， ...]]

 

5. 完成~