Creating polygon grid using GeoPandas

Question

I'm working with a GeoDataFrame of points and I need to create a grid of rectangular polygons.

I know how to have bounds of the GeoDataFrame with gdf.geometry.bounds. Now using this information and with two variables length and width (L, W) of rectangular polygon, how can I create a grid?

Answer 1

Previous answer gives correct result, but I allowed myself to improve the code, to avoid many unnecessary lines, as well as iterating on indexes instead of values on the list.

import geopandas as gpd
from shapely.geometry import Polygon
import numpy as np
points = gpd.read_file('points.shp')

xmin, ymin, xmax, ymax = points.total_bounds

length = 1000
wide = 1200

cols = list(np.arange(xmin, xmax + wide, wide))
rows = list(np.arange(ymin, ymax + length, length))

polygons = []
for x in cols[:-1]:
    for y in rows[:-1]:
        polygons.append(Polygon([(x,y), (x+wide, y), (x+wide, y+length), (x, y+length)]))

grid = gpd.GeoDataFrame({'geometry':polygons})
grid.to_file("grid.shp")

The main idea might be the same, but we are creating now much less useless variables and whole code is clearer to understand

Answer 2

There are many solutions.

One of them

 import geopandas as gpd
 from shapely.geometry import Polygon
 import numpy as np
 points = gpd.read_file('points.shp')
 xmin,ymin,xmax,ymax =  points.total_bounds
 width = 2000
 height = 1000
 rows = int(np.ceil((ymax-ymin) /  height))
 cols = int(np.ceil((xmax-xmin) / width))
 XleftOrigin = xmin
 XrightOrigin = xmin + width
 YtopOrigin = ymax
 YbottomOrigin = ymax- height
 polygons = []
 for i in range(cols):
    Ytop = YtopOrigin
    Ybottom =YbottomOrigin
    for j in range(rows):
        polygons.append(Polygon([(XleftOrigin, Ytop), (XrightOrigin, Ytop), (XrightOrigin, Ybottom), (XleftOrigin, Ybottom)])) 
        Ytop = Ytop - height
        Ybottom = Ybottom - height
    XleftOrigin = XleftOrigin + width
    XrightOrigin = XrightOrigin + width

grid = gpd.GeoDataFrame({'geometry':polygons})
grid.to_file("grid.shp")

You can also truncate the grid (convex hull):

But one of the most interesting is to use the module gpd_lite_toolboox)

Answer 3

One adjustment to the answer from @gene. otherwise, it worked great for my purposes.

I ran this against some points I had to do some binning with, and noticed that the two points that defined the upper and left limits of the extent (used as the origin point for drawing the first and subsequent polygons) do not actually intersect the resulting grid.

To adjust this I just "buffered" out the x and y min and max values by 1 (1 meter), and that did the trick. I understand this is kind of a "Kentucky Windage" solution though, so if there is a better way, let me know!

xmin, ymin, xmax, ymax = gdf.geometry.total_bounds
xmin = xmin - 1
ymin = ymin - 1
xmax = xmax + 1
ymax = ymax + 1