How to Grid a Shapefile in Python

Imagine you want to analyze your data over a country, but you want your analyses to be done on a grid e.g. every 5 square kilometers. Here I’ll show you how to get it done.

Take a shapefile and import it as a geodataframe:

import geopandas
import numpy as np
from shapely.geometry import Polygon

switzerland = geopandas.read_file(DATA_DIRECTORY / "swissBOUNDARIES3D_1_5_LV95_LN02.gpkg", layer="tlm_landesgebiet")

In this case we take a shapefile of Switzerland and specify the layer to be the landesgebeit, or the country area. In this gpkg there are multiple layers, each denoting a different feature e.g. city boundaries, cantonal boundaries etc; therefore, it’s important to choose the right layer.

It’s also important that we only select Switzerland (icc = CH) and none of the surrounding countries, which are included in the layer.

switzerland['icc'].value_counts()
icc
LI    1
CH    1
DE    1
IT    1

Name: count, dtype: int64


switzerland = switzerland.loc[switzerland['icc']=='CH']

Next, we check our coordinate reference system (CRS), get our total bounds of the entire layer, and choose our grid size:

switzerland.crs
<Projected CRS: EPSG:2056>
Name: CH1903+ / LV95
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)


minx,miny,maxx,maxy = switzerland.total_bounds
length = 5000
width = 5000

Since our CRS uses meters we specify 5000 meters for the length and width of our grid cell. Then we prepare our rows and columns, and iterate through them to create a Shapely polygon and append it to a list.

cols = list(np.arange(minx, maxx+ wide, wide))
rows = list(np.arange(miny, maxy+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)]))

With that done, we make a GeoDataFrame from our polygons list, setting the polygons as the geometry column, and set it’s CRS to the original CRS.

grid = geopandas.GeoDataFrame({'geometry': polygons})
grid.set_crs(epsg="2056", inplace=True)

We overlay the grid onto our shapefile with the intersection keyword and then plot the first hundred squares to check the result.

gridded_switzerland = geopandas.overlay(grid, switzerland, how="intersection")

To check if it worked, we can plot the first hundred squares.

There we have it, 5 x 5 km grid cells clipped to the shape of the outline of Switzerland. Finally to save the file as a .gpkg:

test.to_file(DATA_DIRECTORY / "griddedSwitzerland.gpkg")

A gridding function

Here we have all of that as a single function you can use:


def grid_shapefile(cell_len, cell_width, target_gdf):
  """Returns a gridded version of the input GeoDataFrame as a GeoDataFrame. 
  Requires:
  import geopandas
  import numpy as np
  from shapely.geometry import Polygon
  Requires the target_gdf to have a defined crs.
  """
  minx,miny,maxx,maxy = target_gdf.total_bounds
  cols = list(np.arange(minx, maxx+ cell_width, cell_width))
  rows = list(np.arange(miny, maxy+cell_len, cell_len))
  polygons = []
  
  for x in cols[:-1]:
    for y in rows[:-1]:
        polygons.append(Polygon([(x,y), (x+cell_width,y),(x+cell_width, y+cell_len),(x,y+cell_len)]))
  grid = geopandas.GeoDataFrame({'geometry': polygons}).set_crs(target_gdf.crs)
  return geopandas.overlay(grid, target_gdf, how="intersection")

A gridding function

Leave a Reply Cancel reply