python plot a regular grid of points in a polygon

Question

I'd like to plot a regular grid of points within a polygon.

Take the following as a polygon:

import geopandas as gpd
from shapely.geometry import Polygon

x = [0, 100, 150, 100, 50]
y = [0, 115, 200, 250, 50]

polygon_geom = Polygon(zip(x, y))
crs = 'epsg:27700'
polygon = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[polygon_geom])

Using that extent I'd like to plot points on a regular grid, the points are separated by a distance of 10 in N-S and E-W.

Then it should only be a matter of using a contains query to keep those inside the polygon.

Michael Delgado · Answer

I modified your points so the polygon is valid by dropping the last set of points:

import geopandas as gpd
from shapely.geometry import Polygon

x = [0, 100, 150, 100]
y = [0, 115, 200, 250]

polygon_geom = Polygon(zip(x, y))
crs = 'epsg:27700'
polygon = gpd.GeoDataFrame(index=[0], crs=crs, geometry=[polygon_geom])

Not sure exactly the grid spec you're looking for, but you could do something along these lines:

import numpy as np, pandas as pd, shapely.vectorized
import matplotlib.pyplot as plt

# define grid spec from min to max value with spacing of grid_res=10
grid_res = 10

# this isn't necessary in your case, but if your polygon had
# boundaries not falling on grid edges, you might want to ensure
# the grid started on specific cell boundaries or midpoints
minx = np.floor(polygon_geom.bounds[0] / grid_res) * grid_res
maxx = np.ceil(polygon_geom.bounds[2] / grid_res) * grid_res
miny = np.floor(polygon_geom.bounds[1] / grid_res) * grid_res
maxy = np.ceil(polygon_geom.bounds[3] / grid_res) * grid_res

minx, maxx, miny, maxy  # (0.0, 150.0, 0.0, 250.0)

# you can form the grid by setting up the coordinates and then filling
# the grid with numpy.meshgrid
x = np.arange(minx, maxx + grid_res / 2, grid_res)
y = np.arange(miny, maxy + grid_res / 2, grid_res)
XX, YY = np.meshgrid(x, y)

# use shapely.vectorized to find all points within the polygon
in_polygon = shapely.vectorized.contains(polygon_geom, XX, YY)

# filter the points and flatten them to 1D vectors
x_in_polygon = XX[in_polygon].ravel()
y_in_polygon = YY[in_polygon].ravel()

This can then be plotted like this:

In [41]: fig, ax = plt.subplots()
    ...: polygon.plot(color='none', edgecolor='k', ax=ax)
    ...: ax.scatter(x_in_polygon, y_in_polygon, zorder=-1, s=0.4)

python plot a regular grid of points in a polygon

1 Answers