Checking whether a list of points lie within a multipolygon

Question

I am trying to check whether a list of points (which are stored in df.Coords) are found inside my shapefile (multipolygon). I wrote a code below, however, I realised that python only reads in the coordinates of my final layer of my multipolygon.

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) 
    

I tried to search for some code on Stack Overflow (reference: Check if a point falls within a multipolygon with Python) and edited it (as shown below). However, I think I have edited it incorrectly.

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)

How can I check if my list of points (df.Coords) lie within my multipolygon (with 26 layers?)

Answer 1

Create a Geodataframe of the points and use spatial join to find if they are intersecting any polygon:

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="cyan", 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)