How to Build a Solar Suitability Map in Python
On this page
Problem statement
The exclusion mask says where a solar plant may go. It does not say where it should. That needs a suitability map: each buildable cell scored on the factors that make a site cheaper, more productive or easier to permit, and the scores combined into a ranking.
The standard method — a weighted linear combination of normalised factors — is simple, and its two weak points are well known. The weights are opinions, so the ranking has to be tested against other plausible weights. And a factor that barely varies across the region contributes noise with the authority of a map layer. Irradiance is the usual example: it decides which country to build in, and often not which field.
Quick answer
Normalise each factor to 0–1 (1 = best), weight, sum, and apply the mask:
import numpy as np, pandas as pd, rasterio, rasterio.merge
from rasterio.warp import reproject, Resampling
with rasterio.open("siting_mask.tif") as src: # exclusion-mask guide
mask, transform, crs, RES = src.read(1), src.transform, src.crs, src.res[0]
with rasterio.open("dist_substation_km.tif") as src: # distance-to-grid guide
d_sub = src.read(1)
buildable = mask == 1
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, mt = rasterio.merge.merge(srcs, bounds=(-4.32, 39.68, -3.88, 40.02))
dem = np.full(mask.shape, np.nan, dtype="float32")
reproject(mosaic[0], dem, src_transform=mt, src_crs=srcs[0].crs,
dst_transform=transform, dst_crs=crs, resampling=Resampling.bilinear)
dz_dy, dz_dx = np.gradient(dem, RES)
slope = np.hypot(dz_dx, dz_dy) * 100 # percent
aspect = (np.degrees(np.arctan2(dz_dx, dz_dy)) + 180) % 360 # 0 = north, 180 = south
def linear(x, best, worst):
"""1 at `best`, 0 at `worst`, clipped, in either direction."""
return np.clip((x - worst) / (best - worst), 0, 1)
factors = {
"grid": linear(d_sub, best=0, worst=10), # km to a ≥45 kV substation
"slope": linear(slope, best=0, worst=10), # % — mask already cut >10%
"aspect": np.where(slope < 2, 1, linear(np.abs(aspect - 180), best=0, worst=180)),
}
weights = {"grid": 0.6, "slope": 0.25, "aspect": 0.15}
score = sum(weights[k] * factors[k] for k in factors)
score = np.where(buildable, score, np.nan)
print(f"buildable cells: {buildable.sum():,} ({buildable.sum() * RES**2 / 1e6:,.0f} km²)")
print(f"score: median {np.nanmedian(score):.2f}, top decile above "
f"{np.nanpercentile(score, 90):.2f}")
for k, f in factors.items():
print(f" {k:<7} mean {np.nanmean(np.where(buildable, f, np.nan)):.2f} "
f"std {np.nanstd(np.where(buildable, f, np.nan)):.2f}")buildable cells: 948,464 (593 km²)
score: median 0.48, top decile above 0.72
grid mean 0.40 std 0.26
slope mean 0.57 std 0.24
aspect mean 0.64 std 0.32
The spread of each factor over the buildable land matters as much as its weight: a factor with a tiny standard deviation cannot move the ranking whatever weight it is given.
Start from the mask, choose factors that vary, weigh openly
1. Start from the exclusion mask
Score only the cells that survived the hard constraints. Scoring everything and masking afterwards is the same result with more room for mistakes.
2. Choose factors that vary across the region
Before adding a layer, look at its range over the buildable land. Irradiance varies a few percent across a whole country and very little within a 35 km region; distance to a substation runs from zero to 13.5 km over this region's buildable land.
3. Normalise each factor with an explicit rule
A linear ramp between a best and a worst value is transparent and easy to argue about. Write the thresholds down; they are as much a part of the model as the weights.
4. Weight, sum, and keep the parts
Store each normalised factor as well as the total. A score of 0.7 is uninformative until you can see whether it came from a flat field far from the grid or a sloping one next to a substation.
5. Test the ranking against other weights
Recompute with two or three alternative weightings and compare the top-ranked land. If the top sites change completely, the weights are doing the deciding and the map should say so.
6. Aggregate to candidate sites
Plants are built on contiguous parcels. Summarise the score over connected patches of buildable land above a minimum size, and rank sites rather than cells.
Whether irradiance belongs, and how weights decide
The data used on this page
The two inputs are the outputs of the exclusion-mask guide and the distance-to-grid guide: the buildable land of the region 4.30–3.90° W, 39.70–40.00° N south-west of Madrid on a 25 m grid in ETRS89 / UTM 30N, and the straight-line distance from every cell to an OpenStreetMap substation of 45 kV or more, measured with the region padded so that substations beyond its edge count. If either file is missing, the setup rebuilds it with exactly those guides' rules, from OpenStreetMap as the Overpass mirror answered on 9 October 2026. The elevation is the Copernicus DEM, read by the quick answer; Example 1 asks the EU's PVGIS service for PV yields.
import io
import os
import zipfile
import geopandas as gpd
import numpy as np
import osmnx as ox
import rasterio
import rasterio.merge
import 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
# 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-solar-suitability-map-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
ox.settings.max_query_area_size = 50e9 # one Overpass query for the padded box, not several
if not os.path.exists("siting_mask.tif"):
W, S, E, N, CRS, RES = -4.30, 39.70, -3.90, 40.00, "EPSG:25830", 25.0
region = gpd.GeoSeries([box(W, S, E, N)], crs=4326).to_crs(CRS)
x0, y0, x1, y1 = (np.floor(b / RES) * RES for b in region.total_bounds)
shape = (int((y1 - y0) / RES), int((x1 - x0) / RES))
tr = from_origin(x0, y1, RES, RES)
inside = features.rasterize(region, out_shape=shape, transform=tr) == 1
burn = lambda gs, touched=True: features.rasterize(((g, 1) for g in gs), out_shape=shape, transform=tr,
fill=0, all_touched=touched).astype(bool)
near = lambda m, d: ndimage.distance_transform_edt(~m) * RES < d
osm = lambda tags: ox.features_from_bbox((W - 0.03, S - 0.03, E + 0.03, N + 0.03), tags).to_crs(CRS).geometry
excluded = (near(burn(osm({"landuse": ["residential", "industrial", "commercial", "retail"]})), 500)
| near(burn(osm({"highway": ["motorway", "trunk", "primary"]})), 100)
| near(burn(osm({"natural": "water"})), 50))
srcs = [rasterio.open(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")]
mosaic, mt = rasterio.merge.merge(srcs, bounds=(W - 0.01, S - 0.01, E + 0.01, N + 0.01))
dem = np.full(shape, np.nan, dtype="float32")
reproject(mosaic[0], dem, src_transform=mt, src_crs=srcs[0].crs, dst_transform=tr, dst_crs=CRS,
resampling=Resampling.bilinear)
dz_dy, dz_dx = np.gradient(dem, RES)
natura = []
for layer in (0, 1):
gj = requests.get("https://bio.discomap.eea.europa.eu/arcgis/rest/services/ProtectedSites/"
f"Natura2000Sites/MapServer/{layer}/query", params={
"geometry": f"{W},{S},{E},{N}", "geometryType": "esriGeometryEnvelope", "inSR": 4326,
"outSR": 4326, "spatialRel": "esriSpatialRelIntersects", "outFields": "SITECODE",
"f": "geojson"}, timeout=120).json()
natura += list(gpd.GeoDataFrame.from_features(gj["features"], crs=4326).to_crs(CRS).geometry)
candidate = ~(excluded | burn(natura, touched=False) | (np.hypot(dz_dx, dz_dy) * 100 > 10) | ~inside)
labels, _ = ndimage.label(candidate)
keep = np.bincount(labels.ravel()) * RES ** 2 / 1e4 >= 10
keep[0] = False
with rasterio.open("siting_mask.tif", "w", driver="GTiff", dtype="uint8", count=1, width=shape[1],
height=shape[0], crs=CRS, transform=tr, nodata=255, compress="deflate") as dst:
dst.write(np.where(inside, keep[labels], 255).astype("uint8"), 1)
if not os.path.exists("dist_substation_km.tif"):
with rasterio.open("siting_mask.tif") as src:
mask, transform, crs, RES = src.read(1), src.transform, src.crs, src.res[0]
power = ox.features_from_bbox((-4.45, 39.55, -3.75, 40.15), {"power": ["line", "substation"]}).to_crs(crs)
def max_kv(v):
try:
return max(float(x) for x in str(v).split(";")) / 1000
except ValueError:
return np.nan
subs = power[(power["power"] == "substation") & (power["voltage"].map(max_kv) >= 45)].geometry.centroid
MARGIN = int(12_000 / RES)
wide = (mask.shape[0] + 2 * MARGIN, mask.shape[1] + 2 * MARGIN)
burnt = features.rasterize(((g, 1) for g in subs), out_shape=wide, fill=0, all_touched=True,
transform=transform * transform.translation(-MARGIN, -MARGIN))
d_sub = (ndimage.distance_transform_edt(burnt == 0) * 25.0 / 1000)[MARGIN:-MARGIN, MARGIN:-MARGIN]
with rasterio.open("dist_substation_km.tif", "w", driver="GTiff", dtype="float32", count=1,
width=mask.shape[1], height=mask.shape[0], crs=crs, transform=transform,
nodata=-1, compress="deflate") as dst:
dst.write(np.where(mask != 255, d_sub, -1).astype("float32"), 1)
with rasterio.open("siting_mask.tif") as a, rasterio.open("dist_substation_km.tif") as b:
m, d = a.read(1), b.read(1)
print(f"siting_mask.tif: {m.shape[1]} x {m.shape[0]}, buildable {(m == 1).sum() * 625 / 1e6:,.0f} km²; "
f"distance from buildable land to a ≥45 kV substation: median {np.median(d[m == 1]):.1f} km, "
f"max {d[m == 1].max():.1f} km")siting_mask.tif: 1386 x 1349, buildable 593 km²; distance from buildable land to a ≥45 kV substation: median 5.9 km, max 13.5 km
Example 1 — does irradiance belong in the score?
import requests
yields, spread = [], []
for lat in np.linspace(39.72, 39.98, 4):
for lon in np.linspace(-4.28, -3.92, 4):
j = requests.get("https://re.jrc.ec.europa.eu/api/v5_3/PVcalc", params={
"lat": round(lat, 3), "lon": round(lon, 3), "peakpower": 1, "loss": 14,
"angle": 35, "aspect": 0, "outputformat": "json"}, timeout=120).json()
yields.append(j["outputs"]["totals"]["fixed"]["E_y"])
spread.append(j["outputs"]["totals"]["fixed"]["SD_y"]) # year-to-year standard deviation
y = np.array(yields)
print(f"PV yield over 16 points: {y.min():,.0f}–{y.max():,.0f} kWh/kWp, "
f"range {(y.max() - y.min()) / y.mean():.1%}, CV {y.std() / y.mean():.1%}")
print(f"year-to-year standard deviation at a point: {np.mean(spread):.0f} kWh/kWp "
f"({np.mean(spread) / y.mean():.1%}), against {y.max() - y.min():.0f} across the region")PV yield over 16 points: 1,633–1,668 kWh/kWp, range 2.1%, CV 0.6%
year-to-year standard deviation at a point: 49 kWh/kWp (3.0%), against 35 across the region
Two percent between the best and worst corner of the region is less than the year-to-year variation at any one point: PVGIS puts that at 49 kWh/kWp a year, against 35 across the whole region. As a suitability layer it would add a smooth gradient that means nothing; the terrain effect on irradiance is already in the slope and aspect factors.
Example 2 — how much the weights decide
from scipy import stats
schemes = {"grid-led": {"grid": 0.6, "slope": 0.25, "aspect": 0.15},
"equal": {"grid": 1 / 3, "slope": 1 / 3, "aspect": 1 / 3},
"terrain-led": {"grid": 0.3, "slope": 0.45, "aspect": 0.25}}
scores = {n: sum(w[k] * factors[k] for k in factors)[buildable] for n, w in schemes.items()}
top = {n: s >= np.percentile(s, 90) for n, s in scores.items()}
for a, b in (("grid-led", "equal"), ("grid-led", "terrain-led")):
rho = stats.spearmanr(scores[a], scores[b]).statistic
overlap = (top[a] & top[b]).sum() / top[a].sum()
print(f"{a} vs {b}: Spearman {rho:.2f}, top-10% overlap {overlap:.0%}")grid-led vs equal: Spearman 0.85, top-10% overlap 74%
grid-led vs terrain-led: Spearman 0.84, top-10% overlap 72%
A rank correlation of 0.85 sounds like agreement, yet 26–28% of the top decile changes with the weights — and the top decile is the part anybody acts on. Report the overlap, not just the correlation.
Example 3 — rank candidate sites, not cells
from scipy import ndimage
labels, n = ndimage.label(buildable)
ids = np.arange(1, n + 1)
area_ha = ndimage.sum(buildable, labels, ids) * RES**2 / 1e4
mean_score = ndimage.mean(np.nan_to_num(score), labels, ids)
mean_dist = ndimage.mean(d_sub, labels, ids)
sites = pd.DataFrame({"area_ha": area_ha, "score": mean_score, "km_to_substation": mean_dist},
index=ids)
sites = sites[sites.area_ha >= 50].sort_values("score", ascending=False)
sites["capacity_MW"] = sites.area_ha / 100 * 50 # ≈ 51 MW/km², the regional median
print(f"{len(sites)} sites of 50 ha or more, {sites.capacity_MW.sum():,.0f} MW in total")
print(sites.head(8).round({"area_ha": 0, "score": 2, "km_to_substation": 1,
"capacity_MW": 0}).to_string())51 sites of 50 ha or more, 29,010 MW in total
area_ha score km_to_substation capacity_MW
54 162.0 0.89 1.3 81.0
50 67.0 0.87 0.5 33.0
48 104.0 0.86 2.0 52.0
102 60.0 0.80 0.6 30.0
35 771.0 0.75 2.4 385.0
2 397.0 0.75 1.8 198.0
42 162.0 0.73 4.0 81.0
6 3488.0 0.71 2.9 1744.0
The site table is the deliverable. Each row can be checked on the ground, discussed with a landowner, and traced back through the factor layers to the reason it ranks where it does. Note what ranking by mean score does: the four best sites are 60–162 ha, while a 3,488 ha patch — room for 1.7 GW — sits eighth at 0.71 because it averages over good and poor ground.
Normalisation thresholds are model choices
Why normalisation thresholds are model choices
"0 km is best, 10 km or more is worst" is a claim about connection costs. A 20 km threshold would flatten the grid factor and hand influence to slope. The thresholds deserve the same scrutiny as the weights, and the same sensitivity test.
Why weighted sums hide trade-offs
A weighted sum lets a very good score on one factor compensate for a very bad one on another. Sometimes that is right — a slightly steeper field next to a substation — and sometimes it is not. Where a factor has a threshold beyond which no compensation is acceptable, it belongs in the mask.
Why patch-level scores differ from cell-level ones
A large patch averages over its best and worst parts; a small patch can be uniformly excellent. Ranking by mean score favours small homogeneous patches, ranking by total good area favours large ones. Decide which the developer wants, and consider both columns.
Aspect on flat ground and grid capacity
- Aspect on flat ground is noise: gate it on a minimum slope, as above, or it randomly penalises flat fields.
- Land cover and ownership are often the most important soft factors and are missing from open data; cadastral parcels, where available, make the site table far more useful.
- Multi-criteria methods such as AHP formalise the weights but do not remove their subjectivity.
- Resolution: cell scores at 25 m are finer than the data behind distance and land cover justify; aggregate before presenting.
- Grid capacity, where the operator publishes it, should replace or modulate the distance factor.
Internal links
- How to build a siting exclusion mask in Python — the mask this map is built on
- How to measure distance to the grid — the heaviest factor here
- Renewable siting constraints explained — hard constraints against soft factors
- Power density and land use explained — the 50 MW/km² used for capacity
- How to find the optimal tilt and orientation — why aspect is graded
- How to run a site suitability analysis — the general method
FAQ
How do I make a suitability map in Python?
Normalise each factor raster to 0–1, multiply by its weight, sum, and restrict to the cells that pass the exclusion mask.
Should irradiance be a factor?
Across a country, yes. Within a small region it often varies by only a few percent and adds nothing.
How do I choose the weights?
From the costs they represent if you can, from stakeholders if you must — and always test the ranking against alternatives.
How do I check that the result is robust?
Compare the top decile, or the top-ranked sites, under several weighting schemes. Report the overlap.
Why rank sites rather than cells?
Plants are built on contiguous parcels; a ranking of patches is what a developer can act on.
What is the difference between a factor and a constraint?
A constraint excludes land absolutely; a factor ranks the land that remains.