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 2 variables length and width (L, W) of rectangular polygon, how can I create a grid ?


There are many solutions.

One of them

enter image description here

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

enter image description here

You can also truncate the grid (convex hull):

enter image description here

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

  • For calculating number of columns and rows I add abs() because of possibility to have negative coordinates with wgs84 : rows = abs(int(np.ceil((ymax-ymin) / height))) and cols = abs(int(np.ceil((xmax-xmin) / width))) – Tim C. Jan 26 '18 at 14:01
  • 1
    @gene I disagree with you since number will always be positive since it is max minus min, not needing abs. – gcamargo May 18 '18 at 18:45
  • Where is abs in the answer ? – gene May 19 '18 at 19:15

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

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

  • 2
    When you want to work with non-integer coordinates or grid sizes use np.arange instead of range: cols = list(np.arange(xmin, xmax + wide, wide)) ... – quassy Feb 27 at 0:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.