🔧 Fix Renewable energy siting CRS & projections Renewable Energy Siting · step 19 Notebook By · Updated

Siting Mask Removes Everything

On this page

Problem statement

The exclusion mask comes back with no candidate land at all — or with all of it. Every rule looked reasonable, each layer drew correctly on the map, and yet the combined mask is uniformly one value. A mask that removes everything is not a finding about the region; it is almost always one of four mechanical errors: a setback applied in degrees, a layer rasterised in the wrong CRS, a slope computed on a geographic grid, or a boolean inverted somewhere between "excluded" and "allowed".

Each error has a distinctive signature, and a per-layer report of excluded area finds the culprit in one run.

Quick answer

Print the share of the region that each layer excludes on its own. The broken layer is the one at 0% or 100%:

import numpy as np, geopandas as gpd, osmnx as ox
ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"   # overpass-api.de mirror
ox.settings.overpass_rate_limit = False
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 = "EPSG:25830"
region = gpd.GeoSeries([box(WEST, SOUTH, EAST, NORTH)], crs=4326)

built = ox.features_from_bbox((WEST, SOUTH, EAST, NORTH),
                              {"landuse": ["residential", "industrial", "commercial", "retail"]})
built = built[built.geom_type.isin(["Polygon", "MultiPolygon"])]
print(f"{len(built)} built-up polygons, CRS {built.crs}")

def share(zone, crs):
    """Fraction of the region covered by `zone`, as a ratio of areas in the zone's own CRS.

    A 500-degree polygon cannot be reprojected meaningfully, so compare it where it lives.
    """
    r = region.to_crs(crs).iloc[0]
    return zone.intersection(r).area / r.area

wrong = built.geometry.buffer(500).union_all()                 # 500 *degrees*
right = built.to_crs(CRS).geometry.buffer(500).union_all()     # 500 metres
print(f"500 buffered in EPSG:4326: excludes {share(wrong, 4326):.0%} of the region")
print(f"500 buffered in EPSG:25830: excludes {share(right, CRS):.0%} of the region")
287 built-up polygons, CRS epsg:4326
500 buffered in EPSG:4326: excludes 100% of the region
500 buffered in EPSG:25830: excludes 18% of the region

A 500-unit buffer in a geographic CRS is 500 degrees: every town becomes a disc larger than the Earth. GeoPandas warns about it, and the warning is easy to lose in a notebook full of them.

Report each layer, metric CRS, rasterise on the grid

1. Report each layer's exclusion separately

Before combining anything, compute the excluded share of each layer on its own. A layer at 100% or 0% is the bug; a plausible set of percentages narrows the problem to the combination step.

2. Put every vector layer in a metric CRS before buffering

to_crs first, buffer second. Check gdf.crs.is_projected in code, not by eye.

3. Rasterise in the grid's CRS

rasterio.features.rasterize does not reproject. Geometries in degrees burned onto a grid in metres land in a few cells near the origin — or nowhere — and the layer silently excludes nothing.

4. Compute slope on a metric grid

np.gradient(dem, spacing) needs the spacing in the same unit as the elevations. On a raw geographic DEM the spacing is in degrees, and dividing metres by thousandths of a degree gives slopes in the thousands of percent.

5. Name the polarity in the variable

excluded and allowed rather than mask. Rasterise with an explicit fill, combine exclusions with OR, and invert exactly once.

Four causes of an all-excluded or all-allowed siting mask with the signature and fix of each. report each layer’s share firstsignaturecause → fixone layer excludes 100%buffer in degrees → to_crs beforebufferone layer excludes 0.0%rasterised in the wrong CRS →reproject firstslope layer 99% over limitgradient spacing not in metres →reproject DEMplausible map of the townspolarity inverted → name arrays,invert onceEach bug is off by orders of magnitude, so the mask saturates spatialworkflow.io
Four mechanical errors, four signatures — check each layer's share before combining.

Rasterising in the wrong CRS

The data used on this page

The same region south-west of Madrid, 4.30–3.90° W, 39.70–40.00° N, with OpenStreetMap's residential, industrial, commercial and retail land use as the built-up layer, fetched by the quick answer through OSMnx and the maps.mail.ru Overpass mirror. The setup unpacks that answer as saved on 10 October 2026, so the polygon count stays put.

import io
import os
import requests
import zipfile

# OpenStreetMap as Overpass returned it on 10 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-mask-removes-everything-2026-10-10.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

Example 1 — rasterising in the wrong CRS

from rasterio import features
from rasterio.transform import from_origin

x0, y0, x1, y1 = region.to_crs(CRS).total_bounds
RES = 25.0
shape = (int((y1 - y0) / RES), int((x1 - x0) / RES))
transform = from_origin(x0, y1, RES, RES)

def burned(geoms):
    return features.rasterize(((g, 1) for g in geoms), out_shape=shape,
                              transform=transform, fill=0).astype(bool)

print(f"geometries in EPSG:4326 on a metric grid: {burned(built.geometry).mean():.2%} of cells")
print(f"geometries in EPSG:25830:                 {burned(built.to_crs(CRS).geometry).mean():.2%} of cells")
geometries in EPSG:4326 on a metric grid: 0.00% of cells
geometries in EPSG:25830:                 3.53% of cells

No error, no warning: coordinates like (−4.0, 39.9) fall hundreds of kilometres outside a grid whose origin is at (404,000, 4,429,000), and nothing is burned.

Example 2 — slope on a geographic DEM

import rasterio

url = ("https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N39_00_W005_00_DEM/"
       "Copernicus_DSM_COG_10_N39_00_W005_00_DEM.tif")
with rasterio.open(url) as src:
    win = src.window(-4.30, 39.70, -4.00, 40.00)
    dem = src.read(1, window=win).astype("float32")
    step_deg = src.res[0]

for label, spacing in (("spacing in degrees", step_deg),
                       ("spacing = 1 (pixels)", 1.0),
                       ("spacing in metres (≈30 m)", step_deg * 111_320)):
    gy, gx = np.gradient(dem, spacing)
    slope = np.hypot(gx, gy) * 100
    print(f"{label:<27} median slope {np.median(slope):>12,.1f}%, "
          f"cells over 10%: {(slope > 10).mean():.0%}")
spacing in degrees          median slope    531,173.8%, cells over 10%: 99%
spacing = 1 (pixels)        median slope        147.5%, cells over 10%: 99%
spacing in metres (≈30 m)   median slope          4.8%, cells over 10%: 19%

Only the metric spacing gives a believable median: 4.8%, with 19% of cells over the limit. The pixel spacing is a factor of 30 too steep and excludes 99% of the land; the degree spacing makes every cell a cliff of over 500,000%. (The metric line above uses the north–south spacing for both axes; at 40°N the east–west spacing is 23% smaller — reproject the DEM rather than approximating, as in the exclusion-mask guide.)

Example 3 — the per-layer report

gy, gx = np.gradient(dem, 1.0)                                  # the pixel-spacing mistake
steep_wrong = np.hypot(gx, gy) * 100 > 10
checks = {
    "built-up, 500 m buffer": burned(right.geoms if hasattr(right, "geoms") else [right]),
    "built-up, wrong CRS": burned(built.geometry),
    "slope > 10%, pixel spacing": steep_wrong,
}
for name, excluded in checks.items():
    frac = excluded.mean()
    flag = "  <-- suspicious" if frac < 0.001 or frac > 0.95 else ""
    print(f"{name:<28} {frac:>6.1%}{flag}")
built-up, 500 m buffer        18.3%
built-up, wrong CRS            0.0%  <-- suspicious
slope > 10%, pixel spacing    98.6%  <-- suspicious

A report like this, printed every time the mask is built, turns a silent failure into a line that says which layer to look at.

A 500 m buffer around towns compared with a 500-degree buffer that covers the whole map. buffer(500): metres or degreesEPSG:25830 — 18%towns + 500 mEPSG:4326 — 100%towns + 500°Same code, same number; GeoPandas warns once, in a notebook full of warnings spatialworkflow.io
The same number, two units: 500 metres or 500 degrees.

Why these bugs give all or nothing

Why these bugs produce 0% or 100% rather than 20% or 40%

Each one is off by orders of magnitude, not by a few percent. Degrees against metres is a factor of about 100,000; a geographic slope a factor of about 100,000 as well; a CRS mismatch moves geometries hundreds of kilometres. At that scale the mask saturates — which is also what makes the per-layer share such an effective test.

Why rasterize does not reproject

rasterize burns coordinates onto a transform. It has no idea what CRS either is in, and GeoJSON-like geometries carry none. Rasterio leaves reprojection to the caller, which is correct and unforgiving.

Why the inversion bug is the hardest to see

A mask with the polarity reversed looks like a plausible map — it just shows the towns instead of the fields. It survives visual review because it is not uniform. Name the arrays by meaning and assert that the excluded share of the settlements layer is small.

Per-layer excluded share for a correct layer, a layer rasterised in the wrong CRS and a slope computed with pixel spacing. the per-layer report finds the broken layerbuilt-up, 500 m buffer18.3% — plausiblebuilt-up, wrong CRS0.0% — suspiciousslope > 10%, pixel spacing98.6% — suspicious050%100%Flag anything under 0.1% or over 95% before combining the layers spatialworkflow.io
A saturated layer names itself — if each layer is reported on its own.

Degree buffers, nodata and mixed CRSs

  • Buffers of 0.005 "because that's about 500 m" produce ellipses, not a blanket, and are wrong by up to 23% at 40°N — see the constraints guide.
  • Nodata in the DEM propagates NaN into the slope; NaN > 10 is False, but ~(slope <= 10) is True. Write comparisons so NaN falls on the side you intend.
  • Mixed CRSs in one list of layers are common after concatenating sources; assert a single CRS before combining.
  • Empty layers (a query that returned nothing) exclude 0% legitimately — report the feature count alongside the share.
  • Clip to the region last; clipping before buffering loses setbacks from features just outside.

FAQ

Why does my buffer cover everything?

It was computed in a geographic CRS, so the distance was in degrees. Reproject to a metric CRS first.

Why does my rasterised layer have no cells?

The geometries are in a different CRS from the grid's transform. Reproject them before rasterising.

Why are all my slopes over 100%?

The DEM is in degrees and the gradient spacing was not in metres. Reproject the DEM to a metric CRS.

How do I find which layer is wrong?

Report each layer's excluded share on its own. The one at 0% or 100% is the bug.

Why does my mask show the towns instead of the fields?

The polarity is inverted. Name arrays excluded and allowed and invert exactly once.