How to Run a Weighted Site Suitability Analysis in Python

On this page

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,219 × 1,610 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 pipeline ran in 1.7 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 44% 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")
print(f"allowed {allowed[in_area].mean():.1%} of the county ({allowed.sum() * res * res / 1e6:,.0f} km²); "
      f"median score {np.nanmedian(score):.3f}")

allowed 69.8% of the county (1,121 km²); median score 0.706

Measured on the county: these three constraints left 69.8% of it allowed (1,121 km²) and the score's median over allowed cells was 0.706; with the stream setback and the lake blocks as well (step 3) it was 66.2%. The full pipeline in Example 3 — rasterising the layers, the constraint report, two distance transforms, the overlay and extracting patches — took 1.7 s.

Everything after the first box is array arithmetic on one grid, which is why it takes seconds.

One grid, a bare-earth slope, masks, distances, weights

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.shape

Here the reference is EPSG:32145 (NAD83 / Vermont) at 30 m, covering the county plus 2 km: 3,572,590 cells, of which 1,784,179 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 4.7% of the county steep where the DTM found gentle ground, mostly forest edges and building blocks. After weighting, only 56% of the DTM's top 5% of cells were also in the DSM's top 5% — a Jaccard index of 0.396. 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 & ~steep

Measure what each one costs, alone and after the others:

constraint            alone   marginal
water / lake          13.5%    12.8%
stream 30 m buffer     4.4%     3.4%
protected area         9.1%     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 -> metres

Three distance transforms over 3.6 million cells took 0.60 s together.

Rasterise lines with all_touched=True. By default only cells whose centre a line crosses are burned, which left 93,187 road cells against 116,262 with all_touched, and changed the distance value of 24.9% of the county's cells — by a cell or so each, since the median barely moved (234 m against 228 m), but along every road.

5. Standardise every factor to 0–1 with explicit thresholds

The weights in the next step only mean something because the factors already share a scale.

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 is a few array operations. Over allowed cells the score ran from 0.388 at the 5th percentile to 0.950 at the 95th, with a median of 0.705 — 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) + 1

Measured with the top 10% of cells: 2,363 patches, of which 39.4% were a single 0.09 ha cell. Only 172 reached 5 ha, but they held 89.6% of the high-scoring area, and the largest was 614 ha. With rook contiguity (no diagonals) the same cells formed 3,439 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. A couple of hundred polygons, each with an area and a mean score, is a table a planning committee can actually read.

Report the marginal cost of a rule; the standalone figure double-counts overlapping land.

Constraint masks and a weighted overlay on a 30 m grid

The data used on this page

Chittenden County, Vermont, on a 30 m grid in NAD83 / Vermont (EPSG:32145) covering the county plus 2 km. The county and its people are the TIGER/Line 2020 blocks; blocks with no land (ALAND20 == 0) are lake, and "dense" blocks — the population the site should be near — are those with at least 1,000 people per km² of land. Roads are the county's TIGER 2020 road layer. Water (natural=water), streams and rivers (waterway=stream|river) and protected land (boundary=protected_area) come from OpenStreetMap. Slope comes from the USGS 3DEP 1 arc-second bare-earth DTM, with the Copernicus GLO-30 surface model for comparison; the DTM slope is written to slope_dtm_deg.tif, the reference grid. Overpass via the maps.mail.ru mirror; the setup reads its three answers as saved on 9 October 2026, because OpenStreetMap's stream count alone moved by two between runs.

import io
import os
import requests
import urllib.request
import zipfile

import geopandas as gpd
import numpy as np
import osmnx as ox
import rasterio
from rasterio.enums import Resampling
from rasterio.merge import merge
from rasterio.transform import from_origin
from rasterio.warp import reproject
from scipy import ndimage

# 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-site-suitability-analysis-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.overpass_rate_limit = False
ox.settings.http_user_agent = "spatialworkflow.io guide examples (https://www.spatialworkflow.io/)"
CRS, CELL = 32145, 30.0
downloads = {
    "tl_2020_50_tabblock20.zip": "https://www2.census.gov/geo/tiger/TIGER2020/TABBLOCK20/tl_2020_50_tabblock20.zip",
    "tl_2020_50007_roads.zip": "https://www2.census.gov/geo/tiger/TIGER2020/ROADS/tl_2020_50007_roads.zip",
}
for tile in ("n45w074", "n45w073"):
    downloads[f"USGS_1_{tile}.tif"] = ("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/"
                                       f"current/{tile}/USGS_1_{tile}.tif")
for tile in ("N44_00_W074_00", "N44_00_W073_00"):
    name = f"Copernicus_DSM_COG_10_{tile}_DEM"
    downloads[name + ".tif"] = f"https://copernicus-dem-30m.s3.amazonaws.com/{name}/{name}.tif"
for name, url in downloads.items():
    if not os.path.exists(name):
        urllib.request.urlretrieve(url, name)

blocks = gpd.read_file("tl_2020_50_tabblock20.zip").to_crs(CRS)
blocks = blocks[blocks.COUNTYFP20 == "007"]
county = gpd.GeoDataFrame(geometry=[blocks.union_all()], crs=CRS)
lake_blocks = blocks[blocks.ALAND20 == 0]
dense_blocks = blocks[(blocks.ALAND20 > 0) & (blocks.POP20 / (blocks.ALAND20 / 1e6) >= 1000)]
roads = gpd.read_file("tl_2020_50007_roads.zip").to_crs(CRS)

area_ll = county.to_crs(4326).union_all()
def osm(tags, kinds):
    f = ox.features_from_polygon(area_ll, tags=tags).to_crs(CRS)
    return f[f.geom_type.isin(kinds)][["geometry"]].reset_index(drop=True)
water = osm({"natural": "water"}, ["Polygon", "MultiPolygon"])
waterways = osm({"waterway": ["stream", "river"]}, ["LineString", "MultiLineString"])
protected = osm({"boundary": "protected_area"}, ["Polygon", "MultiPolygon"])

west, south, east, north = county.total_bounds + np.array([-2000, -2000, 2000, 2000])
shape = (int(np.ceil((north - south) / CELL)), int(np.ceil((east - west) / CELL)))
transform = from_origin(west, north, CELL, CELL)


def slope_from(names):
    sources = [rasterio.open(n) for n in names]
    mosaic, mt = merge(sources)
    dem = np.full(shape, np.nan, dtype="float32")
    reproject(mosaic[0], dem, src_transform=mt, src_crs=sources[0].crs, src_nodata=sources[0].nodata,
              dst_transform=transform, dst_crs=CRS, resampling=Resampling.bilinear, dst_nodata=np.nan)
    dz_dy, dz_dx = np.gradient(dem, CELL)
    return np.degrees(np.arctan(np.hypot(dz_dx, dz_dy))).astype("float32")


slope_dtm = slope_from([n for n in downloads if n.startswith("USGS_1_")])
slope_dsm = slope_from([n for n in downloads if n.endswith("_DEM.tif")])
with rasterio.open("slope_dtm_deg.tif", "w", driver="GTiff", height=shape[0], width=shape[1], count=1,
                   dtype="float32", crs=CRS, transform=transform, nodata=np.nan) as dst:
    dst.write(slope_dtm, 1)
print(f"{len(water)} water polygons, {len(waterways)} streams and rivers, {len(protected)} protected areas, "
      f"{len(dense_blocks)} dense blocks, {len(roads):,} road segments")

833 water polygons, 2737 streams and rivers, 66 protected areas, 845 dense blocks, 6,546 road segments

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 allowed

Example 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 score

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


import time

start = time.perf_counter()
grid = SuitabilityGrid("slope_dtm_deg.tif")
study = grid.mask(county)
allowed = grid.constraint_report(study, {
    "water / lake": grid.mask(water) | grid.mask(lake_blocks),
    "stream 30 m buffer": grid.mask(waterways, buffer=30),
    "protected area": grid.mask(protected),
    "slope > 15°": slope_dtm > 15,
})
layers = {"slope": slope_dtm, "dist_road": grid.distance_to(roads), "dist_pop": grid.distance_to(dense_blocks)}
score = weighted_overlay(layers, allowed)
sites = candidate_sites(score, grid.transform, grid.crs)
print(f"{time.perf_counter() - start:.2f} s end to end; best five: mean score "
      f"{sites.mean_score.head(5).min():.3f}-{sites.mean_score.head(5).max():.3f}, "
      f"{sites.hectares.head(5).min():.1f}-{sites.hectares.head(5).max():.0f} ha")

water / lake           alone  13.5%   marginal  12.8%
stream 30 m buffer     alone   4.4%   marginal   3.4%
protected area         alone   9.1%   marginal   4.7%
slope > 15°            alone  11.8%   marginal   7.9%
allowed: 66.2%
slope      weight 0.40  mean factor 0.644
dist_road  weight 0.35  mean factor 0.885
dist_pop   weight 0.25  mean factor 0.506
2,365 patches above 0.918; 174 of at least 5.0 ha
1.07 s end to end; best five: mean score 0.965-0.971, 32.9-583 ha

Run end to end with Examples 1 and 2 on the county, it reported 2,365 patches above a score of 0.918 and 174 of at least 5 ha, in about a second including the distance transforms. (The step-by-step run in Example 4 found 2,363 and 172: there the dense blocks were rasterised by cell centre, while distance_to burns them with all_touched=True — the same choice as in step 4.) The best five had mean scores between 0.965 and 0.971 and ranged from 32.9 to 583 ha — near-identical scores, very different sites, which is why the size belongs in the ranking table.

Diagonal neighbours cut the patch count by a third without changing the sites worth naming.

Example 4 — the step-by-step measurements

start = time.perf_counter()
burn = lambda gdf, all_touched=False: grid.mask(gdf, all_touched=all_touched)
in_area = burn(county)
print(f"grid {shape[0]:,} x {shape[1]:,} = {shape[0] * shape[1]:,} cells, {in_area.sum():,} in the county; "
      f"one float32 layer {slope_dtm.astype('float32').nbytes / 1e6:.1f} MB")

# surface model against terrain model
for label, sl in (("Copernicus DSM", slope_dsm), ("USGS 3DEP DTM", slope_dtm)):
    v = sl[in_area]
    print(f"{label:15s} median {np.nanmedian(v):.1f}  p95 {np.nanpercentile(v, 95):.1f}  "
          f"share > 15° {np.mean(v > 15):.1%}")
print(f"DSM steep where the DTM is not: {np.mean((slope_dsm > 15) & (slope_dtm <= 15), where=in_area):.1%}")


def overlay(sl):
    ok = allowed_base & (sl <= 15)
    s = (0.40 * decreasing(sl, 0, 15) + 0.35 * decreasing(dist_road, 0, 2000) + 0.25 * decreasing(dist_pop, 0, 5000))
    return np.where(ok, s, np.nan)


water_m = burn(water) | burn(lake_blocks)
allowed_base = in_area & ~water_m & ~burn(waterways.assign(geometry=waterways.buffer(30))) & ~burn(protected)
road_centre, road_touch = burn(roads), burn(roads, all_touched=True)
dist_road = ndimage.distance_transform_edt(~road_touch) * 30.0
d_centre = ndimage.distance_transform_edt(~road_centre) * 30.0
dist_pop = ndimage.distance_transform_edt(~burn(dense_blocks)) * 30.0
a, b = overlay(slope_dtm), overlay(slope_dsm)
ta, tb = a >= np.nanquantile(a, .95), b >= np.nanquantile(b, .95)
print(f"top 5%: {(ta & tb).sum() / ta.sum():.0%} of the DTM's also in the DSM's; Jaccard {(ta & tb).sum() / (ta | tb).sum():.3f}")
print(f"road cells by centre {road_centre.sum():,}, all_touched {road_touch.sum():,}; distance changed in "
      f"{np.mean(d_centre[in_area] != dist_road[in_area]):.1%} of county cells; medians "
      f"{np.median(d_centre[in_area]):.0f} m and {np.median(dist_road[in_area]):.0f} m")

t = time.perf_counter()
for _ in range(3):
    ndimage.distance_transform_edt(~road_touch)
print(f"three distance transforms: {time.perf_counter() - t:.2f} s")
q = np.nanquantile(a, [.05, .5, .95])
print(f"score over allowed cells: p5 {q[0]:.3f}, median {q[1]:.3f}, p95 {q[2]:.3f}")

good = np.nan_to_num(a, nan=0) >= np.nanquantile(a, 0.90)
for label, structure in (("queen", np.ones((3, 3))), ("rook", None)):
    labels, n = ndimage.label(good, structure=structure)
    ha = np.bincount(labels.ravel())[1:] * 0.09
    print(f"{label}: {n:,} patches, {np.mean(ha == 0.09):.1%} single cells, {np.sum(ha >= 5)} of 5 ha or more holding "
          f"{ha[ha >= 5].sum() / ha.sum():.1%} of the area, largest {ha.max():.0f} ha")

grid 2,219 x 1,610 = 3,572,590 cells, 1,784,179 in the county; one float32 layer 14.3 MB
Copernicus DSM  median 6.0  p95 20.6  share > 15° 14.3%
USGS 3DEP DTM   median 4.5  p95 20.3  share > 15° 11.8%
DSM steep where the DTM is not: 4.7%
top 5%: 56% of the DTM's also in the DSM's; Jaccard 0.396
road cells by centre 93,187, all_touched 116,262; distance changed in 24.9% of county cells; medians 234 m and 228 m
three distance transforms: 0.30 s
score over allowed cells: p5 0.388, median 0.705, p95 0.950
queen: 2,363 patches, 39.4% single cells, 172 of 5 ha or more holding 89.6% of the area, largest 614 ha
rook: 3,439 patches, 51.1% single cells, 182 of 5 ha or more holding 88.0% of the area, largest 563 ha

Grid first, and a cell is not a site

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, the whole pipeline in Example 3 took 1.7 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°.

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 39.6% 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.8% 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.1% 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 2,363 separate patches, and nearly 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 90% of the high-scoring area.

Surface models add false slopes along forest edges, exactly where flat developable land is found.

Straight-line distances, rasterised lines and edge cells

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
How to build a landslide susceptibility index in Python — the same weighted overlay, applied to hazard
How to build a solar suitability map — the same method for solar sites, with a weight-sensitivity test

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 1.7 s.

Should I use a DSM or a DTM for slope?

A DTM. Measured, the Copernicus surface model called 4.7% of the county steep where the DTM found gentle ground, and only 56% 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 172 of 2,363 patches reached 5 ha but held 89.6% 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.8% of the county when counted separately and 33.8% combined; report each rule's marginal cost.