How to Run a Weighted Site Suitability Analysis in Python
Problem statement
You need to find land that suits something — a depot, a solar array, a clinic, a housing allocation — and the criteria are spatial: not too steep, near a road, near the people it serves, not in a lake, a stream corridor or a nature reserve.
A weighted suitability analysis turns that list into a map. Each criterion becomes a raster on one shared grid; the hard rules become masks; the soft preferences become 0–1 factors; a weighted sum ranks every cell; and contiguous patches of high scores become candidate sites.
This guide builds the whole pipeline on a real 30 m grid of Chittenden County, Vermont — 2,220 × 1,611 cells, 1,606 km² inside the county, TIGER roads, OpenStreetMap water and protected areas, 2020 census blocks — and measures each step. Two findings shape the method more than any weight does. The whole analysis ran in 2.55 seconds, so there is no reason not to rerun it for every decision. And swapping the elevation source from a surface model to a bare-earth model, with nothing else changed, replaced 54% of the best cells.
Quick answer
import numpy as np
import rasterio
from rasterio import features
from scipy import ndimage
with rasterio.open("slope_dtm_deg.tif") as src: # the reference grid
slope, transform, shape, res = src.read(1), src.transform, src.shape, src.res[0]
def burn(geoms, all_touched=False):
return features.rasterize(((g, 1) for g in geoms), out_shape=shape,
transform=transform, fill=0, dtype="uint8",
all_touched=all_touched).astype(bool)
in_area = burn(county.geometry)
allowed = in_area & ~burn(water.geometry) & ~burn(protected.geometry) & (slope <= 15)
dist_road = ndimage.distance_transform_edt(~burn(roads.geometry, all_touched=True)) * res
dist_pop = ndimage.distance_transform_edt(~burn(dense_blocks.geometry)) * res
def decreasing(x, good, bad):
return np.clip((bad - x) / (bad - good), 0, 1).astype("float32")
score = (0.40 * decreasing(slope, 0, 15)
+ 0.35 * decreasing(dist_road, 0, 2000)
+ 0.25 * decreasing(dist_pop, 0, 5000))
score = np.where(allowed, score, np.nan).astype("float32")Measured on the county: constraints left 66.2% of the land allowed (1,064 km²), the score's median over allowed cells was 0.798, and the whole run — loading, rasterising seven layers, three distance transforms, the overlay, writing a GeoTIFF and extracting patches — took 2.55 s with a 650 MB peak.
Step-by-step solution
1. Choose one grid and put everything on it
Pick a projected CRS in metres and a cell size, and make one raster the reference. Every other layer is either reprojected onto that grid or rasterised with its transform. Distances, slopes and areas are only meaningful in a projected CRS.
with rasterio.open("slope_dtm_deg.tif") as ref:
profile, transform, shape = ref.profile, ref.transform, ref.shapeHere the reference is EPSG:32145 (NAD83 / Vermont) at 30 m, covering the county plus 2 km: 3,576,420 cells, of which 1,784,173 are inside the county. One float32 layer of that size is 14.3 MB, so a dozen layers sit comfortably in memory.
2. Derive slope from a bare-earth DTM, not a surface model
Slope is usually the first factor, and the choice of elevation data decides it. Copernicus GLO-30 is a surface model: it includes buildings and tree canopy. USGS 3DEP is a terrain model: bare earth. Both were resampled to the same grid and measured against each other:
slope source median p95 share > 15°
Copernicus DSM 6.0 20.6 14.3%
USGS 3DEP DTM 4.5 20.3 11.8%
The DSM called 5.3% of genuinely gentle ground steep, mostly forest edges and building blocks. After weighting, only 46% of the DTM's top 5% of cells were also in the DSM's top 5% — a Jaccard index of 0.305. Use a DTM wherever one exists; see the slope guide for computing it correctly.
3. Turn hard rules into boolean masks
A constraint is a yes/no rule that no score can compensate for. Rasterise each one and combine with &:
water_m = burn(water.geometry) | burn(lake_blocks.geometry)
stream = burn(waterways.geometry.buffer(30)) # 30 m riparian setback
prot = burn(protected.geometry)
steep = slope > 15
allowed = in_area & ~water_m & ~stream & ~prot & ~steepMeasure what each one costs, alone and after the others:
constraint alone marginal
water / lake 13.5% 12.9%
stream 30 m buffer 4.4% 3.4%
protected area 9.0% 4.7%
slope > 15° 11.8% 7.9%
allowed 66.2% of the county
The marginal column is the one to argue about. The protected-area rule looks like it removes 9% of the county, but half of that land was already excluded as steep or wet; on its own it adds only 4.7%.
4. Compute distance factors with a distance transform
scipy.ndimage.distance_transform_edt measures, for every cell, the straight-line distance to the nearest target cell, in one pass:
road = burn(roads.geometry, all_touched=True)
dist_road = ndimage.distance_transform_edt(~road) * res # cells -> metresThree distance transforms over 3.6 million cells took 0.40 s together.
Rasterise lines with all_touched=True. By default only cells whose centre a line crosses are burned, which left 85,018 road cells against 106,231 with all_touched, and changed the distance value of 25.8% of the county's cells — by a cell or so each, since the median barely moved (255 m against 256 m), but along every road.
5. Standardise every factor to 0–1 with explicit thresholds
Each factor needs the same scale and direction, 1 for best:
def decreasing(x, good, bad):
return np.clip((bad - x) / (bad - good), 0, 1).astype("float32")
f_slope = decreasing(slope, 0, 15) # flat is best, 15° is unusable
f_road = decreasing(dist_road, 0, 2000) # on the road is best, 2 km is too far
f_pop = decreasing(dist_pop, 0, 5000)The thresholds come from the problem, not the data. Min–max scaling ties the scale to whichever cell is most extreme, and one outlier or stray NoData value then compresses every other cell into a sliver of the range.
6. Weight, sum and mask
weights = {"slope": 0.40, "road": 0.35, "pop": 0.25}
score = sum(w * f for w, f in zip(weights.values(), (f_slope, f_road, f_pop)))
score = np.where(allowed, score, np.nan).astype("float32")The overlay itself took 0.04 s. Over allowed cells the score ran from 0.506 at the 5th percentile to 0.966 at the 95th, with a median of 0.798 — high, because most of the allowed land is within a couple of kilometres of a road and gently sloping. That skew is normal; rank cells by percentile rather than reading 0.8 as "good".
7. Turn high-scoring cells into candidate sites
Nobody builds on a cell. Threshold the score, label contiguous patches, and keep those big enough to use:
good = np.nan_to_num(score, nan=0) >= np.nanquantile(score, 0.90)
labels, n = ndimage.label(good, structure=np.ones((3, 3))) # queen contiguity
hectares = np.bincount(labels.ravel())[1:] * res * res / 1e4
keep = np.flatnonzero(hectares >= 5) + 1Measured with the top 10% of cells: 3,909 patches, of which 41.4% were a single 0.09 ha cell. Only 271 reached 5 ha, but they held 83.3% of the high-scoring area, and the largest was 531 ha. With rook contiguity (no diagonals) the same cells formed 5,303 patches, over half of them single cells.
8. Write the score and the sites
with rasterio.open("suitability.tif", "w", **dict(profile, dtype="float32", nodata=np.nan)) as dst:
dst.write(score, 1)Vectorise the kept patches with rasterio.features.shapes, dissolve by label, and rank by mean score. The 271 patches became 271 polygons totalling 8,858 ha, which is a table a planning committee can actually read.
Code examples
Example 1 — building constraint masks and factors on a reference grid
import numpy as np
import rasterio
from rasterio import features
from scipy import ndimage
class SuitabilityGrid:
"""Every layer rasterised onto one reference grid."""
def __init__(self, reference_path):
with rasterio.open(reference_path) as ref:
self.profile = ref.profile
self.transform, self.shape = ref.transform, ref.shape
self.res = ref.res[0]
self.crs = ref.crs
def mask(self, gdf, buffer=0, all_touched=False):
geoms = gdf.to_crs(self.crs).geometry
if buffer:
geoms = geoms.buffer(buffer)
return features.rasterize(((g, 1) for g in geoms if g is not None and not g.is_empty),
out_shape=self.shape, transform=self.transform,
fill=0, dtype="uint8", all_touched=all_touched).astype(bool)
def distance_to(self, gdf, all_touched=True):
target = self.mask(gdf, all_touched=all_touched)
if not target.any():
raise ValueError("no target cells: the layer does not overlap the grid")
return (ndimage.distance_transform_edt(~target) * self.res).astype("float32")
def constraint_report(self, study, constraints):
n = study.sum()
allowed = study.copy()
for name, m in constraints.items():
others = np.zeros(self.shape, bool)
for other, m2 in constraints.items():
if other != name:
others |= m2
print(f"{name:22s} alone {np.sum(m & study) / n:6.1%} "
f"marginal {np.sum(m & study & ~others) / n:6.1%}")
allowed &= ~m
print(f"allowed: {allowed.sum() / n:.1%}")
return allowedExample 2 — the weighted overlay with thresholds held in one place
import numpy as np
FACTORS = {
# name: (weight, direction, good, bad)
"slope": (0.40, "decreasing", 0.0, 15.0),
"dist_road": (0.35, "decreasing", 0.0, 2000.0),
"dist_pop": (0.25, "decreasing", 0.0, 5000.0),
}
def standardise(x, direction, good, bad):
if direction == "decreasing":
return np.clip((bad - x) / (bad - good), 0, 1).astype("float32")
return np.clip((x - bad) / (good - bad), 0, 1).astype("float32")
def weighted_overlay(layers, allowed, spec=FACTORS):
total = sum(w for w, *_ in spec.values())
if not np.isclose(total, 1.0):
raise ValueError(f"weights sum to {total}")
score = np.zeros(allowed.shape, dtype="float32")
for name, (weight, direction, good, bad) in spec.items():
factor = standardise(layers[name], direction, good, bad)
print(f"{name:10s} weight {weight:.2f} mean factor {np.nanmean(factor[allowed]):.3f}")
score += np.float32(weight) * factor
score[~allowed] = np.nan
return scoreKeeping weights and thresholds in one dictionary makes them the thing reviewers see and argue about, and the thing you vary in a sensitivity run.
Example 3 — candidate sites from the score surface
import geopandas as gpd
import numpy as np
from rasterio import features
from scipy import ndimage
def candidate_sites(score, transform, crs, top_share=0.10, min_ha=5.0, res=30.0):
"""Patches in the top share of cells, at least min_ha, ranked by mean score."""
threshold = np.nanquantile(score, 1 - top_share)
good = np.nan_to_num(score, nan=0) >= threshold
labels, n = ndimage.label(good, structure=np.ones((3, 3), dtype=int))
hectares = np.bincount(labels.ravel())[1:] * res * res / 1e4
keep = np.flatnonzero(hectares >= min_ha) + 1
print(f"{n:,} patches above {threshold:.3f}; {len(keep)} of at least {min_ha} ha")
kept = np.where(np.isin(labels, keep), labels, 0).astype("int32")
rows = [{"geometry": geom, "properties": {"patch": int(value)}}
for geom, value in features.shapes(kept, mask=kept > 0, transform=transform)]
sites = gpd.GeoDataFrame.from_features(rows, crs=crs).dissolve("patch").reset_index()
sites["hectares"] = sites.area / 1e4
sites["mean_score"] = ndimage.mean(np.nan_to_num(score), labels=labels, index=sites["patch"])
return sites.sort_values("mean_score", ascending=False).reset_index(drop=True)Run end to end with Examples 1 and 2 on the county, it reported 3,951 patches above a score of 0.945 and 270 of at least 5 ha, in 1.52 s including the distance transforms. (The step-by-step run found 3,909 and 271: there the populated blocks were rasterised by cell centre, while distance_to burns them with all_touched=True — the same choice as in step 4, moving the result by about 1%.) The best five had mean scores between 0.977 and 0.979 and ranged from 7.7 to 435 ha — near-identical scores, very different sites, which is why the size belongs in the ranking table.
Explanation
Why the grid comes first
A weighted overlay is element-wise arithmetic: cell [i, j] of every layer must describe the same patch of ground. Rasterising every vector input with the reference transform, and reprojecting every raster onto it, is what makes that true. Skip it and the arithmetic either fails loudly on mismatched shapes or, worse, succeeds on layers with matching shapes but different origins.
Once everything is on one grid, the analysis is fast. Measured: rasterising seven layers took 0.86 s and the overlay 0.04 s. The expensive part of suitability work is deciding the thresholds and weights, not computing them — so rerun freely.
Why the elevation model moved the best sites
A digital surface model records the first thing a radar or laser hits: the canopy top or the roof. Where a forest meets a field, the surface jumps by the height of the trees across one or two 30 m cells, which reads as a slope of 20–30°. The DSM's median elevation here was 4.9 m above the DTM, and 18 m at the 95th percentile.
Those false slopes sit exactly where suitability analyses look — the flat, cleared land at the edge of woodland, near roads. So the DSM both excluded more land (14.3% steep against 11.8%) and scored the remaining edges lower, and the top 5% of cells under the two models shared only 30.5% of their union. No weight choice in this analysis moved the result that much.
Why marginal constraint costs matter
Constraints overlap: stream buffers sit in valleys, protected land is often steep, lakes have wet margins. The standalone share of each rule adds up to 38.7% of the county, while the combined exclusion is 33.8%. When someone asks what dropping the protected-area rule would buy, the honest answer is the marginal figure — 4.7% — not the 9.0% it appears to remove.
Why a cell is not a site
A score surface has texture at the scale of the grid. The top 10% of cells formed 3,909 separate patches, and over 40% of them were a single cell: a gap in the trees, a flat bench on a hillside, a pixel of rounding. Those cannot host anything. The minimum-area rule is a constraint that operates on the result rather than the inputs, and it removed 93% of the patches while keeping 83% of the high-scoring area.
Edge cases or notes
- Distance transforms are straight-line. A site 500 m from a road across a river is not 500 m from access; use network distance for anything that must be driven to.
- Rasterise lines with
all_touched=True, or thin diagonal roads leave gaps that inflate distances along their length. - Stream buffers belong in vector space (
buffer(30)) before rasterising, so the setback is 30 m rather than one cell of whatever size the grid is. - Cells at the grid edge see no roads or towns beyond it. Pad the grid past the study area, as here with 2 km, so distances near the boundary are honest.
- Weights are judgements. Record them with the thresholds and vary them: the result is a ranking under stated assumptions, not a fact.
- The score distribution is skewed by the constraints and thresholds; select by percentile, not by an absolute value such as 0.8.
- Queen and rook contiguity give different site counts — 3,909 against 5,303 patches here — so state which one you used.
- Float32 is enough. A score needs three significant figures, and it halves the memory against float64.
Internal links
- Site suitability explained: constraints, factors and weights — the concepts behind every step here
- Fixing a suitability map that comes out all one value or all NoData — when the overlay goes wrong
- How to calculate slope and aspect from a DEM in Python — the slope factor
- Digital elevation models explained: DEM, DSM and DTM — why the surface model moved the result
- How to rasterize a vector layer in Python — building constraint masks
- How to reproject a raster in Python with Rasterio — putting rasters on the reference grid
- How to convert a raster to a vector in Python — turning patches into sites
- Location analysis explained: catchments, accessibility and site selection — where suitability fits
FAQ
What is a weighted site suitability analysis?
A method that scores every cell of a grid by combining standardised criteria with weights, after removing land that fails hard constraints. The highest-scoring contiguous areas become candidate sites.
How long does a suitability analysis take in Python?
For a county at 30 m, seconds. Measured on 3.6 million cells with seven input layers and three distance transforms, the whole pipeline took 2.55 s with a 650 MB peak.
Should I use a DSM or a DTM for slope?
A DTM. Measured, the Copernicus surface model called 5.3% of gentle ground steep, and only 46% of the best 5% of cells under the DTM were still in the best 5% under the DSM.
How do I turn the suitability raster into sites?
Threshold the score, label contiguous patches with scipy.ndimage.label, and keep patches above a minimum area. Here 271 of 3,909 patches reached 5 ha but held 83.3% of the high-scoring land.
How should I standardise factors?
With explicit thresholds from the problem, such as flat below 0 degrees and unusable above 15, clipped to 0–1. Min–max scaling lets one extreme cell compress the scale for all the others.
Why do my constraint percentages add up to more than the excluded area?
Because constraints overlap. The four rules here removed 38.7% of the county when counted separately and 33.8% combined; report each rule's marginal cost.