Find points within polygons, plot something about the points' attributes

points-in-polys.py

# %%
#### This code is only required to generate some dummy data

import numpy as np
import geopandas as gpd
from shapely.geometry import Point, Polygon
import matplotlib.pyplot as plt

# Make coords for corners of a square
square = np.array([[0, 0], [1, 0], [1, 1], [0, 1]])

# Copy n times to different locations at random sizes
n = 4
squares = np.array([(square + np.random.rand() * 5) * np.random.rand() for i in range(n)])

# Convert to Polygons
df_squares = gpd.GeoDataFrame(geometry=[Polygon(square) for square in squares], crs=4326)

# Create m random points inside
m = 1000
maxx, maxy = squares[:,:,0].max(), squares[:,:,1].max()
df_points = gpd.GeoDataFrame(
    data={"attribute":np.random.rand(m)},
    geometry=[Point(maxx * np.random.rand(), maxy * np.random.rand()) for i in range(m)],
    crs=4326)

# Save
df_squares.to_file("polygons.shp")
df_points.to_file("points.shp")

df_squares.plot(color="none")
df_points.plot(markersize=2, ax=plt.gca())

# %%
#### This is the code you need

import geopandas as gpd
import matplotlib.pyplot as plt

## Find which points are within each polygon

# Load polygons and points
df_polys = gpd.read_file("polygons.shp")
df_points = gpd.read_file("points.shp")

# Spatially join the polygons with the points
# This will match each point to each polygon it resides inside
# This assumes they're in the same CRS. If not, do df.to_crs(target_crs) first
df_polys["poly_geom"] = df_polys.geometry
df_inside = gpd.sjoin(df_points, df_polys)

# Clean up the resultant dataframe
df_inside = (df_inside
    .drop(columns="FID")
    .reset_index()
    .rename(columns={"index_right": "index_poly", "index": "index_point"}))

## Plots

# For each polygon, plot a histogram of the point attributes inside
poly_idxs = df_inside.index_poly.unique()
fig, axes = plt.subplots(ncols=len(poly_idxs), figsize=(len(poly_idxs)*3, 3))

for ax, poly_idx in zip(axes, sorted(poly_idxs)):
    df_poly_points = df_inside[df_inside.index_poly==poly_idx]
    df_poly_points["attribute"].plot(kind="hist", ax=ax)
    ax.set_title(f"Polygon number {poly_idx}")

plt.tight_layout()

# Or maybe you want to plot a summary statistic about the point attributes on a map
# Create the summary statistic, e.g. the mean attribute inside a polygon
df_summary = (df_inside
    .groupby("index_poly")[["attribute"]]
    .mean()
    .rename(columns={"attribute":"attribute_mean"}))

# Merge with the polygon geometry so we can plot on a map
df_summary = df_summary.merge(df_polys[["geometry"]], left_index=True, right_index=True)
df_summary = gpd.GeoDataFrame(df_summary, geometry="geometry", crs=df_polys.crs)

# Plot
df_summary.plot(column="attribute_mean", legend=True)

# Annotate
for idx, row in df_summary.iterrows():
    plt.annotate(
        text=f"{row.attribute_mean:.2f}", 
        xy=row.geometry.centroid.coords[0], 
        horizontalalignment='center')