python按照shp文件提取格点

 参考以下网址提出的脚本

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) 分别用于判断点是否在多边形内部。

1. shapely.geometry.Point(pt).within(geometry.shape(poly)):

   这个方法会为每个点创建一个 Point 对象，并使用 Shapely 提供的 within() 方法来判断点是否在多边形内部。这种方法适用于需要判断单个点是否在多边形内部的情况。

2. 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))

大致能满足按不规则边界提取格点/站点的功能

但是有些格点有遗漏（暂时不确定原因）