Building a Site Suitability Scoring Pipeline with GeoPandas and pvlib #

You have a few hundred candidate solar parcels and four decision layers — a modeled resource surface, a slope raster, the transmission network, and a stack of regulatory exclusion polygons — and you need one defensible number per site: a 0–100 suitability score that a development committee can rank on. This page walks the full pipeline end to end and sits under Solar Irradiance Raster Processing, which produces the analysis-ready Global Horizontal Irradiance (GHI) grid this workflow consumes. The scoring itself is a weighted sum; almost every real failure happens in the plumbing around the sum — layers arriving in different projections, raster values sampled in the wrong frame, one criterion silently dominating because it was never normalized, and banned sites still scoring because the exclusion layer was treated as a soft penalty instead of a hard mask.

Root-cause analysis #

Four compounding causes turn a one-line weighted average into a misleading ranking, and each maps to a specific fix stage below:

1. Layers in mismatched CRS. The sites layer might be in EPSG:4326, the GHI raster in EPSG:32610, and the grid lines in a state-plane CRS. GeoPandas will happily compute a distance() between geometries in different coordinate systems and return a meaningless number in degrees. Every distance, area, and buffer in the score is wrong before it starts, so strict coordinate reference system alignment into one projected metric CRS is the precondition for the whole pipeline.
2. Mixing raster sampling with vector overlay. Rasterizing a continuous GHI surface to polygons, or reprojecting site points into the raster’s frame ad hoc, introduces registration error and doubles the CRS surface area. The correct pattern is to bring every vector layer into one metric frame and then sample the rasters at the site points — never overlay a raster against a vector as if they were the same data model.
3. Unnormalized criteria dominating. Annual GHI is ~1600–2000 (kWh/m²), slope is ~0–30 (degrees), and distance-to-grid is ~0–40 (km). Feed those raw numbers into a weighted sum and GHI’s magnitude swamps everything — the weights become decorative. Each criterion must be min-max normalized to 0–1 before weighting, with cost criteria (slope, distance) inverted.
4. Exclusion zones not hard-masked. A wetland, a setback buffer, or a protected-area polygon is not a penalty to be outweighed by a great resource — it is a disqualifier. If the exclusion enters the score as one more weighted term, a high-GHI site inside a national park can still rank in the top decile. Exclusions and hard constraints must be a multiplicative 0/1 mask applied after the weighted sum.

The composite score for site $i$ is a masked, weighted sum of normalized criteria:

where each raw criterion $c_{i,k}$ is min-max normalized — directly for benefits (yield) and inverted for costs (slope, distance):

and $M_i$ is the product of all hard masks (exclusion membership, maximum slope, maximum grid distance). Because $M_i$ multiplies the whole sum, a single failed constraint drives the score to exactly 0 regardless of how good the resource is.

Pre-flight validation #

Before any sampling or scoring, surface the two structural failures — mismatched CRS and layers that do not spatially overlap — with a compact validator. It refuses to run rather than silently sampling nodata or differencing degrees against metres. This is the guard that makes a CI/CD run fail fast with a precise message.

import numpy as np
import rasterio
import geopandas as gpd


def preflight_suitability_layers(sites_gdf, ghi_raster_path, slope_raster_path,
                                 grid_gdf, exclusion_gdf, target_epsg=32610):
    """Raise on CRS drift or missing coverage before the score is computed."""
    problems = []
    for name, gdf in [("sites", sites_gdf), ("grid", grid_gdf),
                      ("exclusion", exclusion_gdf)]:
        if gdf.crs is None:
            problems.append(f"{name}: CRS is undefined")
        elif gdf.crs.to_epsg() != target_epsg:
            problems.append(f"{name}: EPSG:{gdf.crs.to_epsg()} != target EPSG:{target_epsg}")

    if not (sites_gdf.geom_type == "Point").all():
        problems.append("sites layer must be Point geometries (candidate centroids)")

    sb = sites_gdf.total_bounds  # (minx, miny, maxx, maxy) in the metric CRS
    for path in (ghi_raster_path, slope_raster_path):
        with rasterio.open(path) as src:
            code = src.crs.to_epsg() if src.crs else None
            if code != target_epsg:
                problems.append(f"{path}: raster EPSG:{code} != EPSG:{target_epsg}")
            b = src.bounds
            if sb[0] < b.left or sb[1] < b.bottom or sb[2] > b.right or sb[3] > b.top:
                problems.append(f"{path}: does not fully cover site extent — "
                                "points outside it will sample nodata")

    if problems:
        raise ValueError("Pre-flight failed:\n  - " + "\n  - ".join(problems))

Building the suitability score #

The main function assembles every layer in EPSG:32610, samples GHI and slope at each site point, models a relative PV yield from the sampled GHI with pvlib, adds distance-to-grid and the hard exclusion mask, then normalizes and combines into the 0–100 score. The pvlib step is deliberately a screening-grade relative index — a temperature-corrected pvwatts estimate per unit capacity — not a bankable hourly run; for that, feed the shortlist into Solar PV Yield Simulation with a full ModelChain.

import numpy as np
import geopandas as gpd
import rasterio
import pvlib


def _minmax(values, invert=False):
    """Scale to 0–1; invert for cost criteria. Zero when a criterion has no spread."""
    v = np.asarray(values, dtype="float64")
    lo, hi = np.nanmin(v), np.nanmax(v)
    if not np.isfinite(lo) or (hi - lo) < 1e-9:
        return np.zeros_like(v)            # no discrimination — do not let it drive the score
    scaled = (v - lo) / (hi - lo)
    return 1.0 - scaled if invert else scaled


def relative_pv_yield(ghi_kwh_m2_yr, temp_air_c, wind_ms=2.0, gamma_pdc=-0.0035):
    """Screening-grade relative annual DC yield per kWp from annual GHI."""
    g_eff = np.asarray(ghi_kwh_m2_yr, "float64") * 1000.0 / 8760.0   # kWh/m²/yr → mean W/m²
    temp_cell = pvlib.temperature.faiman(g_eff, temp_air_c, wind_ms)
    dc_w = pvlib.pvsystem.pvwatts_dc(g_eff, temp_cell, pdc0=1000.0, gamma_pdc=gamma_pdc)
    return dc_w * 8760.0 / 1000.0          # relative kWh/kWp/yr index


def score_site_suitability(sites_gdf, ghi_raster_path, slope_raster_path,
                           grid_gdf, exclusion_gdf, target_epsg=32610,
                           weights=None, temp_air_c=15.0,
                           max_slope_deg=15.0, max_grid_km=20.0):
    weights = weights or {"yield": 0.50, "slope": 0.20, "grid": 0.30}
    assert abs(sum(weights.values()) - 1.0) < 1e-6, "criterion weights must sum to 1.0"

    # 1. Assemble one metric frame — reproject every vector layer to EPSG:32610
    sites = sites_gdf.to_crs(target_epsg).copy()
    grid = grid_gdf.to_crs(target_epsg)
    exclusion = exclusion_gdf.to_crs(target_epsg)
    preflight_suitability_layers(sites, ghi_raster_path, slope_raster_path,
                                 grid, exclusion, target_epsg)

    # 2. Sample rasters AT the site points — never overlay raster against vector
    coords = [(geom.x, geom.y) for geom in sites.geometry]
    with rasterio.open(ghi_raster_path) as ghi_src:
        ghi_nodata = ghi_src.nodata
        sites["ghi_kwh_m2_yr"] = [rec[0] for rec in ghi_src.sample(coords)]
    with rasterio.open(slope_raster_path) as slope_src:
        sites["slope_deg"] = [rec[0] for rec in slope_src.sample(coords)]
    if ghi_nodata is not None:
        sites.loc[sites["ghi_kwh_m2_yr"] == ghi_nodata, "ghi_kwh_m2_yr"] = np.nan

    # 3. Resource criterion: model relative yield from sampled GHI with pvlib
    sites["rel_yield"] = relative_pv_yield(sites["ghi_kwh_m2_yr"], temp_air_c)

    # 4a. Grid criterion: distance to nearest transmission line (tie-safe)
    near = gpd.sjoin_nearest(sites[["geometry"]], grid[["geometry"]],
                             how="left", distance_col="grid_dist_m")
    sites["grid_dist_km"] = (near.groupby(near.index)["grid_dist_m"].min()
                             .reindex(sites.index) / 1000.0)

    # 4b. HARD mask: exclusion membership + slope/distance constraints (0 or 1)
    hit = gpd.sjoin(sites[["geometry"]], exclusion[["geometry"]],
                    how="left", predicate="intersects")
    in_excl = sites.index.isin(hit.index[hit["index_right"].notna()])
    hard_mask = (~in_excl
                 & (sites["slope_deg"] <= max_slope_deg)
                 & (sites["grid_dist_km"] <= max_grid_km)).astype("float64")

    # 5. Normalize each criterion to 0–1 BEFORE weighting (cost criteria inverted)
    n_yield = _minmax(sites["rel_yield"])                  # benefit
    n_slope = _minmax(sites["slope_deg"], invert=True)     # cost
    n_grid = _minmax(sites["grid_dist_km"], invert=True)   # cost
    composite = (weights["yield"] * n_yield
                 + weights["slope"] * n_slope
                 + weights["grid"] * n_grid)

    sites["excluded"] = in_excl
    sites["suitability"] = np.round(100.0 * composite * hard_mask, 1)
    ranked = sites.sort_values("suitability", ascending=False, kind="stable")
    ranked["rank"] = range(1, len(ranked) + 1)
    return ranked

Three parameter choices are load-bearing. sjoin_nearest can emit duplicate rows when a site is equidistant from two lines, so the groupby(...).min() collapses ties back to one distance per site — a detail that this vectorized nearest-neighbour search shares with proximity distance calculations across the grid layer. The default weights put half the signal on resource because it is the criterion with the widest bankability spread, but they are an argument precisely so a portfolio can re-run under alternate priorities. And max_slope_deg=15.0 is treated as a hard cutoff rather than a penalty, because construction cost and racking constraints on steep ground are non-negotiable — the same reasoning that drives hillshade and slope analysis for turbine siting. Exclusion polygons should originate from the authoritative regulatory boundary mapping layer so setbacks and protected areas are current.

Fallback routing & performance tuning #

• Sites off the raster footprint. A point outside the GHI grid samples nodata and yields NaN, which _minmax ignores but which then scores 0 — indistinguishable from a genuinely poor site. Route missing-GHI sites to a logged regional mean, or drop them with the fraction recorded, rather than letting a silent NaN masquerade as a real result.
• Degenerate normalization. When a criterion has no spread — every candidate is the same distance from grid — _minmax returns all zeros and that term contributes nothing. Log the min/max range of each criterion; if one collapses, redistribute its weight rather than shipping a score that is secretly built on two criteria.
• Scale the joins with the spatial index. GeoPandas builds an STRtree automatically for sjoin_nearest and sjoin; keep grid_gdf and exclusion_gdf pre-clipped to the study-area bounding box so the tree is small and the nearest query stays fast on tens of thousands of sites.
• Reproject rasters once, up front. If GHI or slope is not already in EPSG:32610, warp the raster once (see Solar Irradiance Raster Processing) instead of reprojecting site points into the raster CRS per run — moving points into a different frame is exactly the mixed-frame bug this pipeline exists to prevent.
• Test weight sensitivity. Run the score under two or three weight sets and export the rank spread. A parcel that lands top-decile under only one weighting is not a robust pick; a parcel that stays top-decile across all of them is.

Downstream validation #

Before the ranked layer is exported to GeoPackage or handed to a mapping portal, gate it with an assertion function suitable for a CI/CD step. It re-checks the invariants that the four failure modes attack — score bounds, the hard mask actually holding, CRS retention, and unique ranks.

def assert_suitability_output(ranked, target_epsg=32610):
    """CI/CD gate: fail the build if the ranked scores are not defensible."""
    assert ranked.crs is not None and ranked.crs.to_epsg() == target_epsg, \
        "output CRS drifted from EPSG:32610"
    score = ranked["suitability"]
    assert score.between(0, 100).all(), "suitability outside 0–100 — normalization broke"
    # Every excluded site MUST score exactly 0 — proves the hard mask held.
    assert (ranked.loc[ranked["excluded"], "suitability"] == 0).all(), \
        "excluded site scored > 0 — hard mask leaked into the weighted sum"
    assert ranked["ghi_kwh_m2_yr"].notna().any(), \
        "all GHI samples are nodata — sites miss the raster footprint"
    assert ranked["rank"].is_unique, "duplicate ranks — tie handling failed"

The exclusion assertion is the one that matters most for permitting review: it is a machine-checkable proof that no disqualified parcel can surface in the shortlist, which is exactly the kind of guarantee a regulator or lender will ask you to demonstrate. Pin geopandas, shapely, rasterio, and pvlib versions in pyproject.toml so a default-behaviour change in a spatial join or a pvlib model cannot silently shift the ranking between runs, and persist the weights and constraint thresholds alongside the output as provenance.

Related #

• Solar Irradiance Raster Processing — produces the CRS-aligned GHI grid this scorer samples.
• Solar PV Yield Simulation — the full pvlib ModelChain run for shortlisted sites once screening narrows the field.
• Automating Hillshade and Slope Analysis for Wind Turbine Siting — derives the slope raster the terrain criterion consumes.
• Proximity Distance Calculations — vectorized nearest-line distance underpinning the grid-proximity criterion.