Site Suitability Explained: Constraints, Factors and Weights

Problem statement

"Where is the best place for it?" is the oldest question in applied GIS, and a suitability analysis is how a map answers it. Every version has the same three parts:

• Constraints — yes/no rules that make a place ineligible: open water, protected land, slopes too steep to build on.
• Factors — graded qualities where more or less is better: gentler slope, a shorter distance to a main road, more people nearby.
• Weights — how much each factor counts relative to the others.

The method is simple arithmetic on a grid. Most of its mistakes come from treating those three parts as interchangeable. That was measured on a real 30 m grid of Chittenden County, Vermont: 1,704,534 land cells, two constraints (water, slope over 15°), three factors (slope, distance to a primary or secondary road, population density within 2 km), and one question — which places come out on top.

what went wrong                                          effect on the answer
constraint left to the factors (water scored, not masked)   7.5% of the top 5% is open water
factors summed in raw units (degrees + metres + people/km²)  8 of the 10 best sites change
0.1 of weight moved between any pair of factors              0 or 1 of the 10 best sites change

The weights — the part people argue about in meetings — mattered least. Treating a constraint as a factor and forgetting to standardise mattered most.

Quick answer

Mask with the constraints, rescale each factor to 0–1 with thresholds you can defend, then take a weighted sum:

import numpy as np

def decreasing(x, good, bad):
    """1 at or below `good`, 0 at or beyond `bad`, linear between."""
    return np.clip((bad - x) / (bad - good), 0, 1)

def increasing(x, bad, good):
    return np.clip((x - bad) / (good - bad), 0, 1)

factors = {
    "slope":  decreasing(slope_deg, good=0, bad=15),
    "road":   decreasing(dist_road_m, good=0, bad=3000),
    "demand": increasing(people_per_km2, bad=0, good=500),
}
weights = {"slope": 1/3, "road": 1/3, "demand": 1/3}

allowed = in_county & ~water & (slope_deg <= 15)          # constraints: a mask
score = sum(weights[k] * factors[k] for k in factors)       # factors: a weighted sum
score = np.where(allowed, score, np.nan)
print(f"constraints removed {1 - allowed.sum() / (in_county & ~water).sum():.1%} of the land; "
      f"{allowed.sum():,} cells ranked")

constraints removed 14.9% of the land; 1,450,846 cells ranked

The two constraints removed 14.9% of the county's land; the weighted sum then ranks the remaining 1,450,846 cells. The arithmetic is cheap — the judgement in the thresholds and weights is not.

Constraints as masks, factors rescaled and weighted

1. Sort every criterion into a constraint or a factor

Ask of each criterion: is there a value beyond which the place is simply not a candidate? If so, it is a constraint. If a worse value can be made up for by a better value elsewhere, it is a factor.

Some criteria are both. Slope was used twice here: slopes over 15° were excluded outright, and below 15° gentler was better.

2. Express constraints as a mask, never as a low score

A constraint that only lowers a score can still be outscored. Leaving Lake Champlain and the county's other water bodies unmasked, and relying on the factors to rank them low, measured:

no mask: top 5% = 76,525 cells; in water blocks 5,716 (7.5%)
mean score: water cells 0.509   land cells 0.394

Water scored higher than land. On the slope factor it is nearly perfect by definition (0.944 against 0.621 for land), and the lake shore is more densely settled than the county average, so its water scores higher on demand too (0.244 against 0.176). On road access land was only slightly ahead (0.385 against 0.338). Two of three factors said the lake was an excellent site (Example 4 prints all of these).

3. Rescale each factor to a common 0–1 range

Slope is measured in degrees, distance in metres and density in people per km². Adding them as they are is adding numbers in different units, and whichever has the widest spread decides the answer (see Explanation).

Use thresholds that mean something: 0° is ideal and 15° is as bad as it gets; 0 m from a main road is ideal and 3 km is as bad as it gets. Thresholds make a factor score readable — 0.5 means "halfway between ideal and unacceptable" — and they do not move when the study area is extended.

4. Choose weights, and write down why

Direct weights are fine for three or four factors. For more, a pairwise comparison (the Analytic Hierarchy Process) forces one judgement at a time: is slope more important than road access, and by how much?

The judgements used here — slope moderately more important than roads (3), slightly more than demand (2), demand slightly more than roads (2) — give:

AHP weights slope/road/demand = [0.54  0.163 0.297]  lambda_max 3.0092  CR 0.0079

A consistency ratio (CR) below 0.1 means the judgements agree with each other. A deliberately circular set — slope beats roads, roads beat demand, demand beats slope — gave CR 1.149, which is the method telling you the weights are incoherent.

5. Combine, then apply the mask

The weighted linear combination is $\text{score} = \sum \text{weight} \times \text{factor}$, and the mask sets excluded cells to NoData. Keep NoData as NaN in a float array so it cannot be mistaken for a score of zero.

6. Select sites at the scale of the decision

A single 30 m cell is not a site. Summarise the score over the footprint you would actually build on — here, a moving 1 km window with at least 90% eligible land — and pick the best non-overlapping windows.

That also makes the result more stable than a cell-level top 5%, because one unusual cell cannot promote a whole site.

7. Test how much the answer depends on the weights

Move 0.1 of weight between each pair of factors, starting from equal weights, and count how many of the ten best sites survive:

+0.1 slope  -0.1 road  : 9/10 kept
+0.1 slope  -0.1 demand: 9/10 kept
+0.1 road   -0.1 slope : 10/10 kept
+0.1 road   -0.1 demand: 10/10 kept
+0.1 demand -0.1 slope : 10/10 kept
+0.1 demand -0.1 road  : 9/10 kept
AHP weights: 8/10 kept
raw units:   2/10 kept

If a site only appears under one weighting, report it as sensitive rather than as a recommendation.

Thresholds, AHP weights and the consistency ratio

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, its water and its people come from the TIGER/Line 2020 blocks: the county is the union of its blocks, water is every block with no land area (ALAND20 == 0) — mostly Lake Champlain — and demand is people per km² within 2 km, from each block's POP20 placed at its representative point. Main roads are the TIGER 2020 primary and secondary roads. Slope comes from the Copernicus GLO-30 surface model, two 1° tiles resampled bilinearly to the grid.

import os
import urllib.request

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.features import rasterize
from rasterio.merge import merge
from rasterio.transform import from_origin
from rasterio.warp import reproject
from scipy import ndimage
from scipy.signal import fftconvolve

CRS, CELL = 32145, 30.0                                   # NAD83 / Vermont, 30 m cells
downloads = {
    "tl_2020_50_tabblock20.zip": "https://www2.census.gov/geo/tiger/TIGER2020/TABBLOCK20/tl_2020_50_tabblock20.zip",
    "tl_2020_50_prisecroads.zip": "https://www2.census.gov/geo/tiger/TIGER2020/PRISECROADS/tl_2020_50_prisecroads.zip",
}
for tile in ("N44_00_W074_00", "N44_00_W073_00"):       # Copernicus GLO-30 DSM, 1 degree tiles
    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"]                           # Chittenden County
west, south, east, north = blocks.total_bounds + np.array([-2000, -2000, 2000, 2000])   # 2 km pad
width, height = int(np.ceil((east - west) / CELL)), int(np.ceil((north - south) / CELL))
transform = from_origin(west, north, CELL, CELL)
grid = dict(out_shape=(height, width), transform=transform)

in_county = rasterize(blocks.geometry, fill=0, default_value=1, dtype="uint8", **grid).astype(bool)
water = rasterize(blocks[blocks.ALAND20 == 0].geometry, fill=0, default_value=1, dtype="uint8", **grid).astype(bool)

roads = gpd.read_file("tl_2020_50_prisecroads.zip").to_crs(CRS)          # primary and secondary roads
on_major_road = rasterize(roads.geometry, fill=0, default_value=1, dtype="uint8",
                          all_touched=True, **grid).astype(bool)

sources = [rasterio.open(n) for n in downloads if n.endswith("_DEM.tif")]
mosaic, mosaic_transform = merge(sources)
dsm = np.full((height, width), np.nan, dtype="float32")
reproject(mosaic[0], dsm, src_transform=mosaic_transform, src_crs=sources[0].crs,
          dst_transform=transform, dst_crs=CRS, resampling=Resampling.bilinear, dst_nodata=np.nan)
dz_dy, dz_dx = np.gradient(dsm, CELL)
slope_deg = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))

# people per km² within 2 km: block populations at their representative points, summed over a 2 km disc
people = np.zeros((height, width), dtype="float64")
pts = blocks[blocks.POP20 > 0]
rows, cols = rasterio.transform.rowcol(transform, pts.representative_point().x, pts.representative_point().y)
np.add.at(people, (np.asarray(rows), np.asarray(cols)), pts.POP20.to_numpy())
r = int(2000 / CELL)
yy, xx = np.mgrid[-r:r + 1, -r:r + 1]
disc = (xx ** 2 + yy ** 2 <= r ** 2).astype("float64")
people_per_km2 = np.clip(fftconvolve(people, disc, mode="same"), 0, None) / (disc.sum() * CELL ** 2 / 1e6)

dist_road_m = ndimage.distance_transform_edt(~on_major_road) * CELL
land = in_county & ~water
print(f"grid {height:,} x {width:,}; {land.sum():,} land cells, {water[in_county].sum():,} water cells; "
      f"{blocks.POP20.sum():,} people")

grid 2,219 x 1,610; 1,704,534 land cells, 79,645 water cells; 168,323 people

Example 1 — factors with explicit thresholds and a constraint mask

import numpy as np
from scipy import ndimage


def distance_to(mask, cell_size):
    """Euclidean distance in map units from every cell to the nearest True cell."""
    return ndimage.distance_transform_edt(~mask) * cell_size


def standardise(x, bad, good):
    """Linear rescale to 0..1; works for increasing (good > bad) and decreasing factors."""
    out = (np.asarray(x, dtype="float32") - bad) / (good - bad)
    return np.clip(out, 0, 1)


def weighted_overlay(factors, weights, allowed):
    """Weighted linear combination; excluded cells become NaN, not 0."""
    total = sum(weights.values())
    if not np.isclose(total, 1.0):
        raise ValueError(f"weights sum to {total:.3f}, not 1")
    missing = set(weights) ^ set(factors)
    if missing:
        raise KeyError(f"factor/weight names do not match: {sorted(missing)}")
    score = sum(weights[name] * factors[name] for name in weights)
    return np.where(allowed, score, np.nan).astype("float32")

factors = {
    "slope": standardise(slope_deg, bad=15, good=0),
    "road": standardise(distance_to(on_major_road, 30), bad=3000, good=0),
    "demand": standardise(people_per_km2, bad=0, good=500),
}
allowed = in_county & ~water & (slope_deg <= 15)
score = weighted_overlay(factors, {"slope": 1/3, "road": 1/3, "demand": 1/3}, allowed)
try:
    weighted_overlay(factors, {"slope": 0.4, "road": 0.4, "demand": 0.3}, allowed)
except ValueError as error:
    print(f"ValueError: {error}")

ValueError: weights sum to 1.100, not 1

On the 2,219 × 1,610 grid — 3.6 million cells — the three factors, the distance transform and the overlay took about 0.2 s on a shared machine (Example 4). The two checks in weighted_overlay catch the silent errors that otherwise produce a plausible-looking map: weights that do not sum to one, as above, and factor names that do not match the weights.

Example 2 — AHP weights and the consistency ratio

import numpy as np

RANDOM_INDEX = {3: 0.58, 4: 0.90, 5: 1.12, 6: 1.24, 7: 1.32, 8: 1.41}


def ahp_weights(matrix, names):
    """Weights from a reciprocal pairwise comparison matrix (principal eigenvector)."""
    A = np.asarray(matrix, dtype=float)
    n = A.shape[0]
    if not np.allclose(A * A.T, 1.0):
        raise ValueError("matrix is not reciprocal: A[j, i] must equal 1 / A[i, j]")
    values, vectors = np.linalg.eig(A)
    i = int(np.argmax(values.real))
    w = np.abs(vectors[:, i].real)
    w /= w.sum()
    lam = values[i].real
    cr = ((lam - n) / (n - 1)) / RANDOM_INDEX[n]
    print(f"lambda_max {lam:.4f}  CR {cr:.4f}  "
          f"{'consistent' if cr < 0.1 else 'INCONSISTENT - revisit the judgements'}")
    return {name: round(float(v), 3) for name, v in zip(names, w)}

print(ahp_weights([[1, 3, 2], [1/3, 1, 1/2], [1/2, 2, 1]], ["slope", "road", "demand"]))
print(ahp_weights([[1, 3, 1/3], [1/3, 1, 3], [3, 1/3, 1]], ["slope", "road", "demand"]))

lambda_max 3.0092  CR 0.0079  consistent
{'slope': 0.54, 'road': 0.163, 'demand': 0.297}
lambda_max 4.3333  CR 1.1494  INCONSISTENT - revisit the judgements
{'slope': 0.333, 'road': 0.333, 'demand': 0.333}

The circular matrix returns perfectly equal weights. Without the CR check, incoherent judgements look like a deliberate choice of equal weighting.

Example 3 — how many of the best sites survive a change of weights

import numpy as np
from scipy import ndimage


def best_sites(score, allowed, k=10, window=33, min_eligible=0.9):
    """The k best non-overlapping windows (window x window cells) by mean score."""
    filled = np.nan_to_num(score, nan=0.0)
    share = ndimage.uniform_filter(allowed.astype("float32"), size=window)
    mean = ndimage.uniform_filter(filled, size=window) / np.maximum(share, 1e-6)
    mean = np.where(share >= min_eligible, mean, -1.0)
    picks = []
    for _ in range(k):
        r, c = np.unravel_index(np.argmax(mean), mean.shape)
        picks.append((int(r), int(c)))
        mean[max(0, r - window):r + window, max(0, c - window):c + window] = -1.0
    return picks


def sites_kept(base, other, tolerance=17):
    return sum(any(abs(a[0] - b[0]) <= tolerance and abs(a[1] - b[1]) <= tolerance
                   for b in other) for a in base)


def weight_sensitivity(factors, weights, allowed, step=0.1):
    base = best_sites(weighted_overlay(factors, weights, allowed), allowed)
    for up in weights:
        for down in weights:
            if up == down or weights[down] < step:
                continue
            w = dict(weights, **{up: weights[up] + step, down: weights[down] - step})
            kept = sites_kept(base, best_sites(weighted_overlay(factors, w, allowed), allowed))
            print(f"+{step} {up:7s} -{step} {down:7s}: {kept}/{len(base)} kept")


equal = {"slope": 1/3, "road": 1/3, "demand": 1/3}
weight_sensitivity(factors, equal, allowed)

+0.1 slope   -0.1 road   : 9/10 kept
+0.1 slope   -0.1 demand : 9/10 kept
+0.1 road    -0.1 slope  : 10/10 kept
+0.1 road    -0.1 demand : 10/10 kept
+0.1 demand  -0.1 slope  : 10/10 kept
+0.1 demand  -0.1 road   : 9/10 kept

On a 30 m grid a 33-cell window is about 1 km, and a 17-cell tolerance treats two windows whose centres are within about 500 m as the same site. The baseline and six shifted runs took about a second together, so there is no reason to skip the test.

Example 4 — the measurements quoted on this page

import time

equal = {"slope": 1/3, "road": 1/3, "demand": 1/3}
ahp = ahp_weights([[1, 3, 2], [1/3, 1, 1/2], [1/2, 2, 1]], ["slope", "road", "demand"])
land = in_county & ~water


def top5(score):
    return score >= np.nanpercentile(score, 95)


def jaccard(a, b):
    return (a & b).sum() / (a | b).sum()


times = []
for _ in range(3):
    start = time.perf_counter()
    f = {"slope": standardise(slope_deg, bad=15, good=0),
         "road": standardise(distance_to(on_major_road, 30), bad=3000, good=0),
         "demand": standardise(people_per_km2, bad=0, good=500)}
    weighted_overlay(f, equal, allowed)
    times.append(time.perf_counter() - start)
print(f"factors + distance transform + overlay: {min(times):.2f} s to {max(times):.2f} s over 3 runs")

# water left to the factors: only the slope constraint applied
nomask = weighted_overlay(factors, equal, in_county & (slope_deg <= 15) | in_county & water)
top = top5(nomask)
print(f"no mask: top 5% = {top.sum():,} cells; in water blocks {(top & water).sum():,} "
      f"({(top & water).sum() / top.sum():.1%})")
print(f"mean score: water cells {np.nanmean(nomask[in_county & water]):.3f}   "
      f"land cells {np.nanmean(nomask[land & (slope_deg <= 15)]):.3f}")
for name in factors:
    print(f"  {name:6}: water {factors[name][in_county & water].mean():.3f}  "
          f"land {factors[name][land & (slope_deg <= 15)].mean():.3f}")

# raw units: each input oriented so that larger is better, then summed with equal weights
raw_inputs = {"slope": -slope_deg, "road": -distance_to(on_major_road, 30), "demand": people_per_km2}
raw = weighted_overlay(raw_inputs, equal, allowed)
base_sites = best_sites(score, allowed)
print(f"sites kept with AHP weights: {sites_kept(base_sites, best_sites(weighted_overlay(factors, ahp, allowed), allowed))}/10; "
      f"raw units: {sites_kept(base_sites, best_sites(raw, allowed))}/10")
top_std, top_raw, top_ahp = top5(score), top5(raw), top5(weighted_overlay(factors, ahp, allowed))
print(f"top 5% Jaccard with the standardised map: raw units {jaccard(top_std, top_raw):.3f}, "
      f"AHP {jaccard(top_std, top_ahp):.3f}")
ok = allowed & np.isfinite(raw)
for name, x in (("slope", slope_deg), ("road", distance_to(on_major_road, 30)), ("demand", people_per_km2)):
    print(f"  raw {name:6}: 0 .. {x[ok].max():,.1f}  sd {x[ok].std():,.1f}  corr with raw score: "
          f"all {np.corrcoef(raw_inputs[name][ok], raw[ok])[0, 1]:.3f}, "
          f"top 5% {np.corrcoef(raw_inputs[name][top_raw], raw[top_raw])[0, 1]:.3f}; "
          f"sd in raw top 5% {x[top_raw].std():,.1f}")

# a Boolean overlay against the weighted one
dist = distance_to(on_major_road, 30)
boolean = allowed & (slope_deg <= 8) & (dist <= 1000) & (people_per_km2 >= 100)
print(f"Boolean overlay keeps {boolean.sum() / allowed.sum():.1%} of the eligible land "
      f"({boolean.sum() * 900 / 1e6:.1f} km²)")
fails = top_std & ~boolean
print(f"weighted top 5% ({top_std.sum():,} cells) failing a Boolean threshold: {fails.sum() / top_std.sum():.1%}, "
      f"{(fails & (dist > 1000)).sum() / top_std.sum():.1%} of them on road distance alone or with others")

# what the best cells look like, and how the factors relate
for name in factors:
    print(f"  {name:6}: mean factor in the top 5% {factors[name][top_std].mean():.2f}, "
          f"across eligible land {factors[name][allowed].mean():.2f}")
F = np.vstack([factors[n][allowed] for n in factors])
c = np.corrcoef(F)
print(f"factor correlations: road-demand {c[1, 2]:.3f}, slope-road {c[0, 1]:.3f}, slope-demand {c[0, 2]:.3f}")
print(f"masked top 5% against the unmasked one: Jaccard {jaccard(top_std, top):.3f}")

lambda_max 3.0092  CR 0.0079  consistent
factors + distance transform + overlay: 0.17 s to 0.19 s over 3 runs
no mask: top 5% = 76,525 cells; in water blocks 5,716 (7.5%)
mean score: water cells 0.509   land cells 0.394
  slope : water 0.944  land 0.621
  road  : water 0.338  land 0.385
  demand: water 0.244  land 0.176
sites kept with AHP weights: 8/10; raw units: 2/10
top 5% Jaccard with the standardised map: raw units 0.559, AHP 0.762
  raw slope : 0 .. 15.0  sd 4.3  corr with raw score: all 0.071, top 5% -0.008; sd in raw top 5% 3.2
  raw road  : 0 .. 9,562.9  sd 2,099.3  corr with raw score: all 0.994, top 5% -0.102; sd in raw top 5% 196.0
  raw demand: 0 .. 2,875.0  sd 251.6  corr with raw score: all 0.440, top 5% 0.951; sd in raw top 5% 630.4
Boolean overlay keeps 13.3% of the eligible land (173.8 km²)
weighted top 5% (72,543 cells) failing a Boolean threshold: 4.2%, 4.0% of them on road distance alone or with others
  slope : mean factor in the top 5% 0.83, across eligible land 0.62
  road  : mean factor in the top 5% 0.88, across eligible land 0.39
  demand: mean factor in the top 5% 0.95, across eligible land 0.18
factor correlations: road-demand 0.534, slope-road 0.069, slope-demand 0.131
masked top 5% against the unmasked one: Jaccard 0.905

Boolean and weighted overlays disagree

Why a Boolean overlay and a weighted overlay disagree

A Boolean overlay turns every criterion into a threshold and keeps the cells that pass all of them. With slope ≤ 8°, within 1 km of a main road and at least 100 people per km², it kept 13.3% of the eligible land (173.8 km²) — and it treats a cell that barely passes all three exactly like one that passes them easily.

A weighted overlay lets a strength in one factor compensate for a weakness in another. Of the weighted top 5% (72,543 cells), 4.2% would have failed a Boolean threshold, nearly all of them (4.0%) because they were more than 1 km from a main road while being flat and densely settled.

Neither is wrong. Compensation is the point of a weighted overlay; if a shortfall should never be compensated, it belongs in the constraints.

Why unstandardised factors are decided by their units

Adding degrees, metres and people per km² gives a number dominated by whichever input has the largest spread:

factor   range              sd        corr with score:  all cells   top 5%
slope    0 .. 15.0 deg        4.3                          0.071     -0.008
road     0 .. 9,562.9 m   2,099.3                          0.994     -0.102
demand   0 .. 2,875.0/km²   251.6                          0.440      0.951

Across the whole county the raw score is almost exactly distance to the road, because metres vary by thousands. Among the cells that reach the top 5%, which all lie near a road, distance varies with a standard deviation of only 196 m — and density, varying by hundreds, takes over.

So the raw-unit map is not merely weighted wrongly: which factor decides the answer changes depending on where you look, and slope decides nothing at all. The top 5% overlapped the standardised top 5% by a Jaccard index of only 0.559, and 2 of the 10 best sites matched.

Why the weights mattered less than expected

Shifting 0.1 of weight between any two factors kept 9 or 10 of the 10 best sites. Moving from equal weights to the AHP weights — slope from 0.33 to 0.54 — kept 8, and the two top-5% sets overlapped by a Jaccard index of 0.762.

The best cells were good on everything. Their mean factor scores were 0.83 for slope, 0.88 for roads and 0.95 for demand, against 0.62, 0.39 and 0.18 across all eligible land. When a site wins on every factor, reweighting reorders the winners without replacing them.

Part of that is geography: road access and demand correlate at 0.534 here, because people live along main roads. Slope was nearly independent of both (0.069 and 0.131). Where factors pull in opposite directions — cheap land far from customers — the same test shows much more movement, which is why it is worth running rather than assuming.

Why water scored better than land

A constraint encodes knowledge that none of the factors contain. Nothing in slope, road distance or population density says "you cannot build on a lake", and on two of the three the lake is attractive: flat by definition, and bordered by the county's densest settlement.

That is the general failure. Any exclusion that coincides with good factor values — floodplains in flat valleys, protected land beside towns, airfields next to highways — rises to the top unless it is masked.

Extent-dependent rescaling and NoData

• Min–max rescaling depends on the extent. Stretching each factor between its own minimum and maximum changes every score when the study area grows; fixed thresholds do not.
• NoData must survive the sum. A NaN factor makes that cell NaN, which is correct; filling NaN with 0 first silently makes missing data look unsuitable.
• A surface model is not bare ground. The slope here came from the Copernicus GLO-30 DSM, which includes tree canopy and buildings; check which elevation model a slope factor was derived from.
• Distance transforms assume square cells in a projected CRS. In degrees they measure nothing in particular.
• "Top 5%" is relative. It always returns 5% of the eligible cells, however poor they are; add an absolute floor if a site must meet a standard.
• Weights are value judgements. AHP makes them consistent, not correct; record who made them.
• Cell-level rankings are noisy. Aggregate to the footprint of the decision before naming sites.
• Suitability is not demand. A suitable site says nothing about whether customers will come; that is a gravity model's question.

Internal links

• How to run a weighted site suitability analysis in Python — the full workflow on rasters
• Fixing a suitability map that comes out all one value or all NoData — when the overlay collapses
• Location analysis explained: catchments, accessibility and site selection — where suitability fits among the other questions
• Slope and aspect explained — how the slope factor is derived
• How to rasterize a vector layer in Python — turning roads and water into grids
• The raster data model explained — NoData, dtype and the transform
• Location-allocation explained: p-median, coverage and what each optimises — choosing several sites together
• Gravity and Huff models explained — whether people will actually come
• Renewable siting constraints explained — exclusions against factors for wind and solar

FAQ

What is the difference between a constraint and a factor?

A constraint excludes a place outright and is applied as a mask. A factor grades places that remain eligible and can be traded off against other factors through the weights.

Do I have to standardise factors before combining them?

Yes. Summed in raw units on the Chittenden grid, the county-wide score correlated 0.994 with distance to road alone, and only 2 of the 10 best sites matched the standardised result.

How should I choose the weights?

State them directly for a handful of factors, or use pairwise comparisons (AHP) and check the consistency ratio is below 0.1. Then test how many of the best sites survive a shift of 0.1 between factors.

Should I use a Boolean overlay or a weighted overlay?

Use Boolean rules for true exclusions and a weighted overlay for everything that can be traded off. A pure Boolean overlay kept 13.3% of the eligible land and could not rank any of it.

Why did my best sites end up in a lake?

Because an exclusion was left to the factors instead of being masked. Unmasked, water filled 7.5% of the top 5% here, since flat and central is what the factors reward.

What cell size should a suitability analysis use?

Fine enough to resolve the narrowest factor that matters, usually slope, and then aggregate the score to the footprint of the real decision before ranking sites.