Polygon Gridding

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()