In many applications, we have a potentially large polygon, and we need to split it into equal-sized squares, for subsequent parallel processing for example.
In this post, we’ll see how to create such squares. We will simply construct a square grid surrounding our polygon, and we will keep the intersecting squares.
Let’s begin writing a function that creates an enclosing grid. The code is pretty straightforward, relying on the
bounds attribute of a polygon, and the
linspace function in Numpy.
from shapely.geometry import Polygon import numpy as np def grid_bounds(geom, delta): minx, miny, maxx, maxy = geom.bounds nx = int((maxx - minx)/delta) ny = int((maxy - miny)/delta) gx, gy = np.linspace(minx,maxx,nx), np.linspace(miny,maxy,ny) grid =  for i in range(len(gx)-1): for j in range(len(gy)-1): poly_ij = Polygon([[gx[i],gy[j]],[gx[i],gy[j+1]],[gx[i+1],gy[j+1]],[gx[i+1],gy[j]]]) grid.append( poly_ij ) return grid
Now, we would like to keep only the squares that intersect our geometry. This can be implemented with a batch of intersection queries.
In order to do this efficiently, we will use a feature of Shapely, called prepared geometry.
To apply the
intersects function to a list of squares, we then rely on Python’s
filter built-in function (see official Python docs).
from shapely.prepared import prep def partition(geom, delta): prepared_geom = prep(geom) grid = list(filter(prepared_geom.intersects, grid_bounds(geom, delta))) return grid
Now, it’s time to run a test and plot the results. In the example below we are using Geopandas just for it’s plotting function.
import geopandas as gpd import matplotlib.pyplot as plt geom = Polygon([[0,0],[0,2],[1.5,1],[0.5,-0.5],[0,0]]) grid = partition(geom, 0.1) fig, ax = plt.subplots(figsize=(15, 15)) gpd.GeoSeries(grid).boundary.plot(ax=ax) gpd.GeoSeries([geom]).boundary.plot(ax=ax,color="red") plt.show()