⚙️ How-To Renewable energy siting Raster & rasterio Renewable Energy Siting · step 12 Notebook By · Updated

How to Build a Siting Exclusion Mask in Python

On this page

Problem statement

Vector buffers are the natural first way to apply siting rules, and they stop scaling as soon as a rule depends on terrain. Slope comes from a DEM, and a DEM is a raster; so does aspect, irradiance and every other suitability factor that follows. The practical representation of an exclusion mask is therefore a raster on a fixed metric grid: one cell per 25 m, true where a plant may go.

A raster mask also makes the next steps trivial — area is a cell count, overlaps are boolean algebra, and patches too small to build on can be removed with connected-component labelling.

Quick answer

Put every constraint on one grid, as a distance or a value, and combine them with boolean operators:

import numpy as np, geopandas as gpd, osmnx as ox, rasterio, requests
from rasterio import features
from rasterio.transform import from_origin
from rasterio.warp import reproject, Resampling
from scipy import ndimage
from shapely.geometry import box

ox.settings.http_user_agent = "spatialworkflow-examples/1.0"
WEST, SOUTH, EAST, NORTH = -4.30, 39.70, -3.90, 40.00
CRS, RES = "EPSG:25830", 25.0

region = gpd.GeoSeries([box(WEST, SOUTH, EAST, NORTH)], crs=4326).to_crs(CRS)
x0, y0, x1, y1 = (np.floor(b / RES) * RES for b in region.total_bounds)
width, height = int((x1 - x0) / RES), int((y1 - y0) / RES)
transform = from_origin(x0, y1, RES, RES)
inside = features.rasterize(region, out_shape=(height, width), transform=transform) == 1

def burn(geoms):
    """1 where a feature is, 0 elsewhere, on the common grid."""
    return features.rasterize(((g, 1) for g in geoms), out_shape=(height, width),
                              transform=transform, fill=0, all_touched=True).astype(bool)

def distance_m(mask):
    """Distance in metres from every cell to the nearest True cell."""
    return ndimage.distance_transform_edt(~mask) * RES

def osm(tags):
    return ox.features_from_bbox((WEST - 0.03, SOUTH - 0.03, EAST + 0.03, NORTH + 0.03),
                                 tags).to_crs(CRS).geometry

built_g = osm({"landuse": ["residential", "industrial", "commercial", "retail"]})
roads_g = osm({"highway": ["motorway", "trunk", "primary"]})
water_g = osm({"natural": "water"})
built, roads, water = burn(built_g), burn(roads_g), burn(water_g)
print(f"grid: {width} × {height} cells of {RES:.0f} m, {inside.sum() * RES**2 / 1e6:,.0f} km² inside")

excluded = (distance_m(built) < 500) | (distance_m(roads) < 100) | (distance_m(water) < 50)
print(f"excluded by settlements, roads and water: "
      f"{(excluded & inside).sum() * RES**2 / 1e6:,.0f} km²")
grid: 1386 × 1349 cells of 25 m, 1,140 km² inside
excluded by settlements, roads and water: 253 km²

The distance transform replaces thousands of buffer operations with one pass over the grid. It sees features just outside the region only where the grid extends beyond the region: the OSM query is padded by more than the largest setback, but burn draws only inside the grid, which here is the bounding box of the projected region. Where that matters, burn onto a grid widened by the largest setback and crop, as the distance-to-grid guide does for substations that can be kilometres outside.

One grid, padded fetches, distance transforms, slope

1. Define one grid and snap to it

Choose the CRS, the resolution and the origin once, and put every layer on that grid. Snapping the origin to a multiple of the resolution keeps later layers — irradiance, land cover — aligned without resampling artefacts.

2. Pad the data fetch by the largest setback

A town 400 m outside the region still excludes land inside it. Download features for the region plus a margin, and clip at the end.

3. Rasterise features, then take distances

rasterize with all_touched=True burns every cell a line crosses; the Euclidean distance transform then gives the distance from each cell to the nearest burned one. A setback becomes a threshold on distance.

4. Derive slope from a DEM on the same grid

Reproject the DEM onto the grid (bilinear), compute the gradient in metres per metre, and express it in percent or degrees — whichever the rule uses.

5. Combine, then remove small patches

OR the exclusions together, invert to get the candidates, and drop connected patches smaller than a buildable plant — a few hectares for solar.

6. Save the mask as a GeoTIFF with its grid

Every suitability layer that follows reads the same transform. Store the mask as uint8 with a nodata value for outside the region.

Raster layers on one 25 m grid combined with OR into a siting exclusion mask. one grid, one boolean per ruledistance to built-up < 500 mdistance transformORdistance to major road < 100 mdistance transformORdistance to water < 50 mdistance transformORNatura 2000rasterised polygonsORslope > 10%Copernicus GLO-30ORpatches < 10 handimage.label= excluded1,386 × 1,349 cells; a new setback is a new threshold, not a new geometry spatialworkflow.io
Each rule becomes one boolean layer on the same grid; the mask is their OR.

A full exclusion mask from fetched data

The data used on this page

Everything is fetched by the page's own code: OpenStreetMap settlements, main roads, water, power lines through OSMnx; the Copernicus GLO-30 DEM from its public bucket; and Natura 2000 sites from the EEA's map service, for a region south-west of Madrid, 4.30–3.90° W, 39.70–40.00° N. The setup points OSMnx at an Overpass mirror, because the main Overpass server was refusing connections when this page was run, and unpacks the mirror's answers as saved on 9 October 2026: the water layer alone gained a feature between two runs.

import io
import os
import osmnx as ox
import requests
import zipfile

# OpenStreetMap as Overpass returned it on 9 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-siting-exclusion-mask-python-2026-10-09.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"
ox.settings.requests_timeout = 600

Example 1 — slope from the Copernicus DEM

import rasterio.merge

tiles = [f"https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N39_00_{t}_00_DEM/"
         f"Copernicus_DSM_COG_10_N39_00_{t}_00_DEM.tif" for t in ("W005", "W004")]
srcs = [rasterio.open(t) for t in tiles]
mosaic, mosaic_transform = rasterio.merge.merge(
    srcs, bounds=(WEST - 0.01, SOUTH - 0.01, EAST + 0.01, NORTH + 0.01))

dem = np.full((height, width), np.nan, dtype="float32")
reproject(mosaic[0], dem, src_transform=mosaic_transform, src_crs=srcs[0].crs,
          dst_transform=transform, dst_crs=CRS, resampling=Resampling.bilinear)

dz_dy, dz_dx = np.gradient(dem, RES)
slope_pct = np.hypot(dz_dx, dz_dy) * 100
print(f"elevation {np.nanmin(dem[inside]):.0f}–{np.nanmax(dem[inside]):.0f} m")
for limit in (5, 10, 15):
    share = (slope_pct[inside] > limit).mean()
    print(f"slope over {limit:>2}%: {share:.1%} of the region")
elevation 407–1078 m
slope over  5%: 48.1% of the region
slope over 10%: 19.4% of the region
slope over 15%: 8.8% of the region

Four fifths of the region is under 10%, and half of it under 5%. A 10% limit is typical for fixed-tilt solar; single-axis trackers usually need less. Copernicus GLO-30 is a surface model, so a little of the steep ground is the edge of buildings and tree rows rather than terrain.

Example 2 — the full mask and its patches

def natura2000():
    parts = []
    for layer in (0, 1):
        gj = requests.get("https://bio.discomap.eea.europa.eu/arcgis/rest/services/"
                          f"ProtectedSites/Natura2000Sites/MapServer/{layer}/query", params={
            "geometry": f"{WEST},{SOUTH},{EAST},{NORTH}", "geometryType": "esriGeometryEnvelope",
            "inSR": 4326, "outSR": 4326, "spatialRel": "esriSpatialRelIntersects",
            "outFields": "SITECODE", "f": "geojson"}, timeout=120).json()
        parts.append(gpd.GeoDataFrame.from_features(gj["features"], crs=4326))
    return gpd.pd.concat(parts).to_crs(CRS).geometry

natura = natura2000()
protected = features.rasterize(((g, 1) for g in natura), out_shape=(height, width),
                               transform=transform, fill=0).astype(bool)
steep = slope_pct > 10
exclusion = excluded | protected | steep | ~inside
candidate = ~exclusion

labels, n = ndimage.label(candidate)
sizes = np.bincount(labels.ravel()) * RES**2 / 1e4          # hectares
keep = sizes >= 10
keep[0] = False
buildable = keep[labels]

km2 = lambda m: m.sum() * RES**2 / 1e6
print(f"candidate land: {km2(candidate):,.0f} km² in {n:,} patches")
print(f"patches of 10 ha or more: {keep.sum():,}, {km2(buildable):,.0f} km²")
print(f"lost to fragmentation: {km2(candidate) - km2(buildable):,.0f} km² "
      f"in {n - keep.sum():,} slivers")
candidate land: 612 km² in 5,655 patches
patches of 10 ha or more: 108, 593 km²
lost to fragmentation: 19 km² in 5,547 slivers

5,547 of the 5,655 patches are too small to hold a plant: slivers between road buffers, corners cut off by a water setback. Together they are 19 km² that a raw eligible-land figure would count as usable, while 108 patches hold almost all the rest.

Example 3 — write it out, and check against the vector result

profile = {"driver": "GTiff", "dtype": "uint8", "count": 1, "width": width,
           "height": height, "crs": CRS, "transform": transform, "nodata": 255,
           "compress": "deflate"}
out = np.where(inside, buildable.astype("uint8"), 255).astype("uint8")
with rasterio.open("siting_mask.tif", "w", **profile) as dst:
    dst.write(out, 1)

no_slope = (excluded | protected) & inside
print(f"without slope or patch filters: {km2(inside) - km2(no_slope):,.0f} km² available")
setbacks = gpd.GeoSeries(list(built_g.buffer(500)) + list(roads_g.buffer(100)) +
                         list(water_g.buffer(50)) + list(natura), crs=CRS)
area = region.iloc[0]
vector_km2 = (area.area - setbacks.union_all().intersection(area).area) / 1e6
print(f"vector buffers on the same four layers:  {vector_km2:,.0f} km² available")
print(f"from {len(roads_g):,} road features, {len(built_g):,} land-use polygons, "
      f"{len(water_g):,} water features and {len(natura):,} Natura 2000 sites")
without slope or patch filters: 767 km² available
vector buffers on the same four layers:  768 km² available
from 1,383 road features, 352 land-use polygons, 239 water features and 5 Natura 2000 sites

The raster and vector routes agree to within a square kilometre: 767 against 768 km² on identical inputs. The raster side is the smaller because all_touched rasterisation errs towards exclusion by up to a cell along every feature. Renewable siting constraints explained builds the vector version with more layers — railways and rivers — so its total is not directly comparable.

From the 1,140 km² region through each exclusion to 593 km² of buildable patches. Toledo region, km²region1,140−253after settlements, roads, water887−120after Natura 2000767−155after slope > 10%612−19 in 5,547 sliverspatches ≥ 10 ha593108 patches hold almost all of the buildable land spatialworkflow.io
Every rule takes its share; fragmentation takes the last one.

Distance transforms rather than buffers

Why a distance transform instead of buffers

Buffering 1,383 road features and 352 land-use polygons and unioning the results is slow and produces heavy geometries. Rasterising once and computing the Euclidean distance transform is linear in the number of cells, and a change of setback is a new threshold, not a new geometry operation — which is what makes sensitivity analysis cheap.

Why all_touched=True

The default rasterisation burns only cells whose centre is inside a polygon or on a line's path, so a thin road can disappear between cell centres. Burning every touched cell errs on the side of exclusion by up to one cell width, which is the right direction for a constraint.

Why patch size is a constraint

A plant needs a contiguous site for its fence, internal roads and connection. A 2 ha field wedged between two setbacks is eligible land that nobody will build on, and leaving it in the total inflates the technical potential.

Share of the Toledo region over three slope limits, from the Copernicus 30 m DEM. how much land each slope limit removesslope over 5%48.1% of the regionslope over 10%19.4% — the usual fixed-tilt limitslope over 15%8.8%030%60%Trackers need gentler ground: the limit choice is a technology choice spatialworkflow.io
A 5% limit for trackers removes two and a half times the land a 10% limit does.

Resolution, DEM type and memory

  • Resolution: 25 m is fine for regional siting; setbacks under 50 m or narrow features need finer cells or a vector pass.
  • DEM type: Copernicus GLO-30 is a surface model — forests and buildings raise it — so slopes at forest edges can be exaggerated.
  • Slope rules are sometimes in degrees: 10% is 5.7°.
  • Memory: the distance transform needs a float array the size of the grid; for a country, tile it with an overlap equal to the largest setback.
  • Aspect belongs in the suitability score, not in the mask, unless a rule explicitly forbids north-facing ground.

FAQ

How do I create an exclusion mask for solar or wind siting?

Rasterise each constraint layer onto one metric grid, threshold distances and slopes by the rules, OR the results, and invert.

Should I use vector buffers or rasters?

Rasters once terrain or any other gridded factor is involved; vectors for a few simple layers or very precise setbacks.

What resolution should the grid have?

About a quarter of the smallest setback you care about. 25 m suits regional siting.

How do I remove patches that are too small?

Label connected components with scipy.ndimage.label, measure each one's area, and drop those below the minimum site size.

Why does my raster mask exclude more than my vector buffers?

all_touched=True rasterisation and cell-size rounding err towards exclusion by up to one cell.

Where do I get slope for the mask?

From a DEM such as Copernicus GLO-30, reprojected to the mask's grid before computing the gradient.