参考以下网址提出的脚本
Python按shp文件提取格点和插值图_python shp-CSDN博客
原脚本在判断大量点是否在相同多边形内部时,速度较慢,原因是shapely.geometry.Point(pt).within(geometry.shape(poly)) 方法会为每个点创建一个 Point
对象,并使用 Shapely 提供的 within()
方法来判断点是否在多边形内部。这种方法适用于需要判断单个点是否在多边形内部的情况。
根据ChatGPT,shapely.geometry.Point(pt).within(geometry.shape(poly))
和 shapely.ops.cascaded_union.contains(point)
分别用于判断点是否在多边形内部。
-
shapely.geometry.Point(pt).within(geometry.shape(poly))
:这个方法会为每个点创建一个
Point
对象,并使用 Shapely 提供的within()
方法来判断点是否在多边形内部。这种方法适用于需要判断单个点是否在多边形内部的情况。 -
shapely.ops.cascaded_union.contains(point)
:这个方法首先对一组多边形进行联合操作,然后再判断点是否在结果多边形内部。
cascaded_union
是对整个多边形集合进行了优化,因此对于大量的点或者多边形,这种方法可能会更有效率。这种方法适用于需要反复判断大量点是否在相同多边形内部的情况。
在实际使用中,如果只是针对少量点进行判断,两种方法的性能差别可能不太明显。但如果需要对大量点进行判断,建议使用 cascaded_union.contains(point)
,因为它只需要进行一次联合操作,然后反复使用结果多边形来判断点是否在内部,相比于每次都创建 Point
对象来说,效率更高。
方案2在判断大量点时,计算效率显著提升(方案1提取一个国家的格点可能要一个多小时,方案2仅需1-2分钟)
现附上一个示例脚本留档,方便后续处理类似的工作
import numpy as np
import shapefile
from shapely.ops import cascaded_union
from shapely.geometry import Polygon
from shapely.geometry import Point
# 用函数表示
def mask_grid(shp, ilon, ilat): # 本函数实现shp文件范围内的格点选取
shp = shapefile.Reader(shp) #, encoding='ISO-8859-1'
rec = shp.shapeRecords()
polygon = []
for r in rec:
polygon.append(Polygon(r.shape.points)) # 获取各个地市点边界
poly = cascaded_union(polygon) # 并集
ext = list(poly.exterior.coords) # 外部点,即最外面轮廓(省界)
# x = [i[0] for i in ext]
# y = [i[1] for i in ext]
# plt.plot(x, y, 'r')
# 提取边界内格点
grid_lon, grid_lat = np.meshgrid(ilon, ilat)
flat_lon = grid_lon.flatten()
flat_lat = grid_lat.flatten()
in_shape_points = []
for lon, lat in zip(flat_lon, flat_lat):
print(lon,lat)
point = Point(lon, lat)
if poly.contains(point): # 使用shapely.geometry.Point.within()方法检查点是否在多边形内
in_shape_points.append((lon, lat))
sel_lon, sel_lat = list(zip(*in_shape_points)) # 将结果重新组合成两个列表
return sel_lon, sel_lat
shp='路径+BR_Pais_2022.shp'
lon = np.arange(-90, -30+1)
lat = np.arange(-45, 15+1)
sel_lon, sel_lat=mask_grid(shp, lon, lat)
print(len(sel_lon))
print(len(sel_lat))
大致能满足按不规则边界提取格点/站点的功能
但是有些格点有遗漏(暂时不确定原因)