检查点列表是否位于多边形内（Python）

Question

我正在尝试检查在我的形状文件（多重多边形）中是否找到了点列表（存储在

df.Coords

中）。我在下面编写了一段代码，但是，我意识到 python 只读取多边形最后一层的坐标。

import shapefile
from shapely.geometry import shape, Point

# read your shapefile
r = shapefile.Reader("C:\\filename.shp")

# get the shapes
shapes = r.shapes()

# build a shapely polygon from your shape
hold = []
for k in range(len(shapes)-1): 
    polygon = shape(shapes[k])
    hold.append(polygon)


    for x in df.Coords: 

        if polygon.contains(Point(x)):

            hold.append(x) 

我尝试在stackoverflow上搜索一些代码（参考：https://gis.stackexchange.com/questions/208546/check-if-a-point-falls-within-a-multipolygon-with-python/208574 ）并进行编辑（如下所示）。但是，我认为我编辑错误。

from shapely.geometry import Polygon, Point, MultiPolygon
import shapefile 

polygon = shapefile.Reader("C:\\filename.shp") 

polygon = polygon.shapes() 

shpfilePoints = [ shape.points for shape in polygon ]

#print(shpfilePoints)


polygons = shpfilePoints
hold = []

for polygon in polygons:
    poly = Polygon(polygon)
    hold.append(poly)
    for x in hold: 

        for y in df.Coords:
            if x.contains(Point(y)):
                hold.append(y) 
print('hold',hold)

如何检查我的点列表 (

df.Coords

) 是否位于我的多边形内（有 26 层？）

Answer 1

创建点的地理数据框并使用空间连接来查找它们是否与任何多边形相交：

import numpy as np
import geopandas as gpd
poly_df = gpd.read_file(r"C:\Users\bera\Desktop\gistest\world.geojson")

#Create a point dataframe with 1000 rows
xMin, yMin, xMax, yMax = poly_df.total_bounds
x_coords = np.random.uniform(low=xMin, high=xMax, size=1000)
y_coords = np.random.uniform(low=yMin, high=yMax, size=1000)
points = gpd.points_from_xy(x=x_coords, y=y_coords)
point_df = gpd.GeoDataFrame(geometry=points, crs=poly_df.crs)

#Plot the polygons and points
ax = poly_df.plot(figsize=(10,10), color="lightgreen", zorder=1)
point_df.plot(ax=ax, color="gray", zorder=2, markersize=10)

#Spatial join the polygon df to the points. 
  #All attributes from the polygon df will be appended to the output point df
sj = point_df.sjoin(poly_df, how="left", predicate="intersects")

#Select points with a not nan polygon index = points within polygons and plot them
sj.loc[~gpd.pd.isna(sj.index_right)].plot(ax=ax, markersize=2, color="red", zorder=2)