通过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. 完成~

posted @ 2019-04-01 17:52  keeptg  阅读(5363)  评论(0编辑  收藏  举报