Source code: Rooftop solar potential from LiDAR in Python: pvlib, PVGIS and shading for every roof

"""Rooftop solar potential in Python: Environment Agency 1 m LiDAR DSM + Overture Maps building
footprints + PVGIS weather + pvlib sun positions, for every roof in a 1 km x 0.8 km block of
Levenshulme, Manchester.

Pipeline: roof cells from footprints and the normalised DSM -> Horn slope and aspect -> horizon
angles in 72 directions by ray marching over the DSM (shading from neighbours, trees and the roof
itself) -> a "sun map" of the PVGIS typical meteorological year -> annual beam + diffuse +
reflected irradiation per 1 m roof cell -> kWh per year per building with stated assumptions.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/rooftop-solar-potential-lidar-python/
Data: Environment Agency LIDAR Composite DSM and DTM 1 m (Open Government Licence v3.0);
Overture Maps buildings (ODbL, includes OpenStreetMap and Microsoft ML footprints); PVGIS
(European Commission JRC, free to reuse with attribution).
"""
import sys
import time
import warnings
from importlib.metadata import version
from pathlib import Path

import duckdb
import geopandas as gpd
import numpy as np
import pandas as pd
import pvlib
import rasterio
import requests
import shapely
from pyproj import Transformer
from rasterio.features import rasterize
from rasterio.transform import from_origin
from scipy import ndimage
import matplotlib
matplotlib.use("Agg")   # render to files only; no GUI backend
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm, ListedColormap

warnings.filterwarnings("ignore")
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})
DATA = Path("data"); DATA.mkdir(exist_ok=True)
OUT = Path("output"); OUT.mkdir(exist_ok=True)

# Study area in British National Grid (EPSG:27700): terraces, semis, schools and sheds in Levenshulme
AOI = (386800, 393700, 387800, 394500)        # xmin, ymin, xmax, ymax: 1 km x 0.8 km of roofs to assess
FETCH = (386000, 393000, 388000, 395000)      # 2 km x 2 km of LiDAR, so shading from up to 300 m away is seen
EA_WCS = {
    "dtm": ("lidar-composite-digital-terrain-model-dtm-1m",
            "13787b9a-26a4-4775-8523-806d13af58fc__Lidar_Composite_Elevation_DTM_1m"),
    "dsm": ("lidar-composite-digital-surface-model-first-return-dsm-1m",
            "df4e3ec3-315e-48aa-aaaf-b5ae74d7b2bb__Lidar_Composite_Elevation_FZ_DSM_1m"),
}

# Assumptions: change these and re-run
PANEL_EFF = 0.20        # module efficiency at standard test conditions (typical 2026 mono-Si: 20-22%)
PERF_RATIO = 0.80       # performance ratio: inverter, cabling, temperature, soiling and mismatch losses
COVERAGE = 0.70         # share of a usable roof plane that panels can actually cover (edges, vents, access)
MAX_SLOPE = 60          # degrees; steeper cells are walls, dormers or LiDAR edge artefacts
MIN_IRR = 800           # kWh/m2/yr; cells below this are not worth covering in Manchester
ALBEDO = 0.20           # ground reflectance for the reflected component
MIN_HEIGHT = 2.5        # m above the DTM for a cell to count as roof
HOUSEHOLD_KWH = 2700    # Ofgem typical domestic electricity consumption (medium, 2024 onwards)

AZ = np.arange(2.5, 360, 5.0)                  # 72 horizon directions, clockwise from north
STEPS = np.r_[np.arange(2, 30, 1), np.arange(30, 100, 2), np.arange(100, 301, 5)]  # ray-march distances (m)
NOISE = 0.3             # m; LiDAR noise allowance so a rough roof does not shade itself


def section(title):
    print(f"\n=== {title} ===")


section("versions")
print(f"python {sys.version.split()[0]} | pvlib {pvlib.__version__} | rasterio {rasterio.__version__} | "
      f"duckdb {duckdb.__version__} | geopandas {gpd.__version__} | numpy {np.__version__} | "
      f"pandas {pd.__version__} | scipy {version('scipy')}")

# ---------------------------------------------------------------- Step 1: EA LiDAR DSM and DTM
section("Step 1: Environment Agency 1 m LiDAR DSM and DTM")


def ea_wcs(product, x0, y0, size=1000):
    """One 1 km x 1 km GeoTIFF from the Environment Agency WCS (no key needed)."""
    service, cov = EA_WCS[product]
    url = (f"https://environment.data.gov.uk/spatialdata/{service}/wcs?service=WCS&version=2.0.1"
           f"&request=GetCoverage&CoverageId={cov}&format=image/tiff"
           f"&subset=E({x0},{x0 + size})&subset=N({y0},{y0 + size})")
    for attempt in range(4):
        try:
            r = requests.get(url, timeout=180)
            r.raise_for_status()
            return r.content
        except requests.RequestException:
            if attempt == 3:
                raise
            time.sleep(5 * (attempt + 1))


def ea_mosaic(product):
    """Download 1 km tiles and paste them into one 1 m array on the FETCH grid."""
    tf = from_origin(FETCH[0], FETCH[3], 1, 1)
    a = np.full((FETCH[3] - FETCH[1], FETCH[2] - FETCH[0]), np.nan, "float32")
    for x0 in range(FETCH[0], FETCH[2], 1000):
        for y0 in range(FETCH[1], FETCH[3], 1000):
            tile = DATA / f"ea_{product}_{x0}_{y0}.tif"
            if not tile.exists():
                tile.write_bytes(ea_wcs(product, x0, y0))
            with rasterio.open(tile) as src:
                t = src.read(1).astype("float32")
                t[(t == src.nodata) | (t < -1000)] = np.nan
                c, r = (int(round(v)) for v in ~tf * (src.bounds.left, src.bounds.top))
                a[r:r + t.shape[0], c:c + t.shape[1]] = t
    return a, tf


dsm, TF = ea_mosaic("dsm")
dtm, _ = ea_mosaic("dtm")
print(f"grid {dsm.shape[1]} x {dsm.shape[0]} cells of 1 m | no data: DSM {np.isnan(dsm).mean() * 100:.2f}%, "
      f"DTM {np.isnan(dtm).mean() * 100:.2f}%")
dsm = np.where(np.isnan(dsm), dtm, dsm)
dtm = np.where(np.isnan(dtm), dsm, dtm)
dsm = ndimage.median_filter(dsm, size=3)        # knock down single-pixel LiDAR noise before taking slopes
ndsm = dsm - dtm
print(f"DTM {np.nanmin(dtm):.1f}..{np.nanmax(dtm):.1f} m | DSM {np.nanmin(dsm):.1f}..{np.nanmax(dsm):.1f} m | "
      f"cells over {MIN_HEIGHT} m above ground: {(ndsm > MIN_HEIGHT).mean() * 100:.1f}%")

# AOI window inside the FETCH grid
R0, C0 = FETCH[3] - AOI[3], AOI[0] - FETCH[0]
H, W = AOI[3] - AOI[1], AOI[2] - AOI[0]
A = (slice(R0, R0 + H), slice(C0, C0 + W))
ATF = from_origin(AOI[0], AOI[3], 1, 1)

# ---------------------------------------------------------------- Step 2: Overture building footprints
section("Step 2: Overture Maps building footprints with DuckDB")
to_ll = Transformer.from_crs(27700, 4326, always_xy=True)
lon0, lat0 = to_ll.transform(AOI[0], AOI[1])
lon1, lat1 = to_ll.transform(AOI[2], AOI[3])
BBOX = (min(lon0, lon1) - 0.001, min(lat0, lat1) - 0.001, max(lon0, lon1) + 0.001, max(lat0, lat1) + 0.001)
LON, LAT = to_ll.transform((AOI[0] + AOI[2]) / 2, (AOI[1] + AOI[3]) / 2)
print(f"AOI centre {LAT:.4f} N, {abs(LON):.4f} W | query bbox {', '.join(f'{v:.4f}' for v in BBOX)}")
bld_path = DATA / "overture_buildings.gpkg"
if not bld_path.exists():
    release = requests.get("https://stac.overturemaps.org/catalog.json", timeout=60).json()["latest"]
    con = duckdb.connect()
    con.sql("INSTALL spatial; LOAD spatial; INSTALL httpfs; LOAD httpfs; SET s3_region = 'us-west-2';")
    xmin, ymin, xmax, ymax = BBOX
    overlaps = f"bbox.xmin <= {xmax} AND bbox.xmax >= {xmin} AND bbox.ymin <= {ymax} AND bbox.ymax >= {ymin}"
    files = [r[0] for r in con.sql(f"""
        SELECT assets.aws.alternate.s3.href FROM read_parquet('https://stac.overturemaps.org/{release}/collections.parquet')
        WHERE collection = 'building' AND {overlaps}""").fetchall()]
    df = con.sql(f"""
        SELECT id, subtype, class, height, num_floors,
               list_filter(sources, lambda s: s.property = '')[1].dataset AS source,
               ST_AsWKB(geometry) AS wkb
        FROM read_parquet({files}) WHERE {overlaps}""").df()
    df["release"] = release
    gpd.GeoDataFrame(df.drop(columns="wkb"), geometry=shapely.from_wkb(df["wkb"].apply(bytes)),
                     crs=4326).to_file(bld_path)
bld = gpd.read_file(bld_path).to_crs(27700)
print(f"Overture release {bld['release'].iloc[0]}: {len(bld):,} buildings in the query box")
c = bld.geometry.centroid
bld = bld[(c.x > AOI[0]) & (c.x < AOI[2]) & (c.y > AOI[1]) & (c.y < AOI[3])].reset_index(drop=True)
bld["footprint_m2"] = bld.area
print(f"buildings with a centroid in the AOI: {len(bld):,} | sources: "
      + ", ".join(f"{k} {v:,}" for k, v in bld["source"].value_counts().items()))
print(f"footprint area: median {bld['footprint_m2'].median():.0f} m2, total {bld['footprint_m2'].sum() / 1e4:.1f} ha "
      f"({bld['footprint_m2'].sum() / (H * W) * 100:.0f}% of the AOI)")

bid = rasterize(((g, i + 1) for i, g in enumerate(bld.geometry)), out_shape=(H, W), transform=ATF,
                fill=0, dtype="int32")
inner = ndimage.binary_erosion(bid > 0, iterations=1)      # drop the edge cell: footprints and LiDAR never align perfectly
roof = inner & (ndsm[A] > MIN_HEIGHT)
print(f"footprint cells over {MIN_HEIGHT} m in the nDSM: {(ndsm[A] > MIN_HEIGHT)[bid > 0].mean() * 100:.1f}% "
      f"| roof cells after a 1-cell inward buffer: {roof.sum():,}")

# ---------------------------------------------------------------- Step 3: slope and aspect
section("Step 3: roof slope and aspect (Horn's method)")


def horn(z, res=1.0):
    """Slope (degrees) and aspect (degrees clockwise from north, the direction the slope faces)."""
    p = np.pad(z, 1, mode="edge")
    a, b, c = p[:-2, :-2], p[:-2, 1:-1], p[:-2, 2:]
    d, f = p[1:-1, :-2], p[1:-1, 2:]
    g, h, i = p[2:, :-2], p[2:, 1:-1], p[2:, 2:]
    dzdx = ((c + 2 * f + i) - (a + 2 * d + g)) / (8 * res)       # rise towards the east
    dzdy = ((a + 2 * b + c) - (g + 2 * h + i)) / (8 * res)       # rise towards the north
    slope = np.degrees(np.arctan(np.hypot(dzdx, dzdy)))
    aspect = (np.degrees(np.arctan2(-dzdx, -dzdy)) + 360) % 360  # downslope direction
    return slope, aspect


slope_f, aspect_f = horn(dsm)
slope, aspect = slope_f[A], aspect_f[A]
s, a = slope[roof], aspect[roof]
pitched = s >= 10
print(f"roof cells: flat (<10 deg) {(~pitched).mean() * 100:.1f}% | pitched 10-{MAX_SLOPE} deg "
      f"{(pitched & (s <= MAX_SLOPE)).mean() * 100:.1f}% | steeper {(s > MAX_SLOPE).mean() * 100:.1f}%")
print(f"pitched roof slope: median {np.median(s[pitched & (s <= MAX_SLOPE)]):.0f} deg, "
      f"interquartile {np.percentile(s[pitched & (s <= MAX_SLOPE)], 25):.0f}-"
      f"{np.percentile(s[pitched & (s <= MAX_SLOPE)], 75):.0f} deg")
facing = pd.Series(pd.cut((a[pitched] + 22.5) % 360, np.arange(0, 361, 45), labels=["N", "NE", "E", "SE", "S", "SW", "W", "NW"], right=False))
print("pitched roof cells by the way they face: "
      + ", ".join(f"{k} {v * 100:.0f}%" for k, v in facing.value_counts(normalize=True, sort=False).items()))

# ---------------------------------------------------------------- Step 4: horizon angles
section("Step 4: horizon angles in 72 directions by ray marching over the DSM")
rr, cc = np.nonzero(roof)
WF = dsm.shape[1]
flat_idx = (rr + R0) * WF + (cc + C0)
zf = dsm.ravel()
z0 = zf[flat_idx]
t0 = time.perf_counter()
hor = np.zeros((len(AZ), len(rr)), "float32")             # horizon elevation angle in degrees (>= 0)
for k, az in enumerate(AZ):
    best = np.zeros(len(rr), "float32")
    sa, ca = np.sin(np.radians(az)), np.cos(np.radians(az))
    for d in STEPS:
        off = -int(round(d * ca)) * WF + int(round(d * sa))
        np.maximum(best, (zf[flat_idx + off] - z0 - NOISE) / d, out=best)
    hor[k] = np.degrees(np.arctan(best))
print(f"{len(AZ)} directions x {len(STEPS)} steps (2-300 m) for {len(rr):,} roof cells in "
      f"{time.perf_counter() - t0:.1f} s")

# Sky view relative to an unobstructed plane of the same tilt (isotropic sky, cos^2 weighting)
beta, asp = np.radians(s), np.radians(a)
h_plane = np.arctan(np.maximum(-np.tan(beta)[None, :] * np.cos(np.radians(AZ)[:, None] - asp[None, :]), 0))
h_obs = np.maximum(np.radians(hor), h_plane)
svf_rel = (np.cos(h_obs) ** 2).mean(0) / (np.cos(h_plane) ** 2).mean(0)
del h_plane, h_obs
print(f"relative sky view: median {np.median(svf_rel):.2f}, 10th percentile {np.percentile(svf_rel, 10):.2f}")

# ---------------------------------------------------------------- Step 5: PVGIS weather and the sun map
section("Step 5: PVGIS typical meteorological year and a sun map")
tmy_path = DATA / "pvgis_tmy.csv"
if not tmy_path.exists():
    tmy, meta = pvlib.iotools.get_pvgis_tmy(LAT, LON, usehorizon=False, map_variables=True)
    tmy.to_csv(tmy_path)
tmy = pd.read_csv(tmy_path, index_col=0, parse_dates=True)
times = tmy.index + pd.Timedelta(minutes=30)             # hourly values: put the sun in the middle of the hour
sp = pvlib.solarposition.get_solarposition(times, LAT, LON)
tmy["zenith"], tmy["azimuth"], tmy["elevation"] = sp["apparent_zenith"].values, sp["azimuth"].values, \
    sp["apparent_elevation"].values
print(f"TMY {len(tmy):,} hours | annual GHI {tmy['ghi'].sum() / 1000:.0f} kWh/m2, DNI {tmy['dni'].sum() / 1000:.0f}, "
      f"DHI {tmy['dhi'].sum() / 1000:.0f} | diffuse share of GHI {tmy['dhi'].sum() / tmy['ghi'].sum() * 100:.0f}%")
sun = tmy[(tmy["elevation"] > 0) & (tmy["dni"] > 0)].copy()
z = np.radians(sun["zenith"]); azr = np.radians(sun["azimuth"])
sun["SA"] = sun["dni"] * np.cos(z)                       # Wh/m2 that a horizontal plane would get
sun["SB"] = sun["dni"] * np.sin(z) * np.cos(azr)
sun["SC"] = sun["dni"] * np.sin(z) * np.sin(azr)
sun["az_bin"] = (sun["azimuth"] // 5).astype(int)
sun["el_bin"] = (sun["elevation"] // 2).astype(int)
sun["w_el"] = sun["dni"] * sun["elevation"]
sunmap = sun.groupby(["az_bin", "el_bin"])[["SA", "SB", "SC", "dni", "w_el"]].sum().reset_index()
sunmap["el"] = sunmap["w_el"] / sunmap["dni"]            # DNI-weighted sun elevation in the bin
print(f"{len(sun):,} sunlit hours grouped into {len(sunmap)} sky cells (5 deg azimuth x 2 deg elevation) | "
      f"beam on a horizontal plane {sunmap['SA'].sum() / 1000:.0f} kWh/m2")


def annual_irradiance(slope_deg, aspect_deg, horizon_deg, svf):
    """Annual beam + diffuse + reflected irradiation (kWh/m2 of tilted surface)."""
    b, g = np.radians(slope_deg), np.radians(aspect_deg)
    cb, sb, cg, sg = np.cos(b), np.sin(b), np.cos(g), np.sin(g)
    beam = np.zeros(np.shape(slope_deg), "float64")
    for row in sunmap.itertuples():
        inc = row.SA * cb + sb * (row.SB * cg + row.SC * sg)
        lit = horizon_deg[row.az_bin] < row.el
        beam += np.where(lit, np.maximum(inc, 0), 0)
    diffuse = tmy["dhi"].sum() * (1 + cb) / 2 * svf
    reflected = tmy["ghi"].sum() * ALBEDO * (1 - cb) / 2
    return beam / 1000, diffuse / 1000, reflected / 1000


# ---------------------------------------------------------------- Step 6: check the model on an open plane
section("Step 6: check on an unshaded 35 deg south-facing plane")
b0, d0, r0 = annual_irradiance(np.array([35.0]), np.array([180.0]), np.zeros((72, 1)), np.array([1.0]))
ours = float(b0[0] + d0[0] + r0[0])
dni_h = tmy["dni"].where(tmy["zenith"] < 90, 0)
poa = pvlib.irradiance.get_total_irradiance(35, 180, tmy["zenith"], tmy["azimuth"], dni_h, tmy["ghi"], tmy["dhi"],
                                            albedo=ALBEDO, model="isotropic")["poa_global"]
dni_extra = pvlib.irradiance.get_extra_radiation(times)
poa_perez = pvlib.irradiance.get_total_irradiance(35, 180, tmy["zenith"], tmy["azimuth"], dni_h, tmy["ghi"],
                                                  tmy["dhi"], dni_extra=dni_extra.values, albedo=ALBEDO,
                                                  airmass=pvlib.atmosphere.get_relative_airmass(tmy["zenith"]).values,
                                                  model="perez")["poa_global"]
loc = pvlib.location.Location(LAT, LON)
cs = loc.get_clearsky(times, model="ineichen")
cs_poa = pvlib.irradiance.get_total_irradiance(35, 180, tmy["zenith"], tmy["azimuth"],
                                               cs["dni"].where(tmy["zenith"].values < 90, 0).values,
                                               cs["ghi"].values, cs["dhi"].values, albedo=ALBEDO)["poa_global"]
hourly_path = DATA / "pvgis_hourly_35S_2014_2023.csv"
if not hourly_path.exists():
    hourly, _ = pvlib.iotools.get_pvgis_hourly(LAT, LON, start=2014, end=2023, surface_tilt=35, surface_azimuth=180,
                                               components=True, usehorizon=False, pvcalculation=True,
                                               peakpower=1, loss=14)
    hourly.to_csv(hourly_path)
hourly = pd.read_csv(hourly_path, index_col=0, parse_dates=True)
hourly["poa"] = hourly[["poa_direct", "poa_sky_diffuse", "poa_ground_diffuse"]].sum(axis=1)
n_years = hourly.index.year.nunique()
pvgis_poa = hourly["poa"].sum() / 1000 / n_years
pvgis_kwh_kwp = hourly["P"].sum() / 1000 / n_years
print(f"sun-map model (this script)        : {ours:6.0f} kWh/m2/yr  (beam {b0[0]:.0f}, diffuse {d0[0]:.0f}, "
      f"reflected {r0[0]:.0f})")
print(f"pvlib hourly, isotropic sky        : {poa.sum() / 1000:6.0f} kWh/m2/yr")
print(f"pvlib hourly, Perez sky            : {poa_perez.sum() / 1000:6.0f} kWh/m2/yr")
print(f"PVGIS SARAH {hourly.index.year.min()}-{hourly.index.year.max()} mean       : {pvgis_poa:6.0f} kWh/m2/yr")
print(f"pvlib Ineichen clear sky (no cloud): {cs_poa.sum() / 1000:6.0f} kWh/m2/yr")
print(f"specific yield: ours {ours * PERF_RATIO:.0f} kWh/kWp (PR {PERF_RATIO}) | PVGIS PV model, 14% system loss "
      f"{pvgis_kwh_kwp:.0f} kWh/kWp")

# ---------------------------------------------------------------- Step 7: annual irradiation on every roof cell
section("Step 7: annual irradiation on every roof cell, with shading")
t0 = time.perf_counter()
beam, diffuse, refl = annual_irradiance(s, a, hor, svf_rel)
irr = beam + diffuse + refl
b_open, _, _ = annual_irradiance(s, a, np.zeros((72, 1)), np.ones(1))
print(f"{len(sunmap)} sky cells x {len(s):,} roof cells in {time.perf_counter() - t0:.1f} s")
shade_loss = 1 - beam.sum() / b_open.sum()
print(f"annual irradiation on roof cells: median {np.median(irr):.0f} kWh/m2, 90th percentile "
      f"{np.percentile(irr, 90):.0f}, max {irr.max():.0f}")
print(f"beam lost to shading from neighbours, trees and chimneys: {shade_loss * 100:.1f}% of the unshaded beam")
by_face = pd.DataFrame({"face": facing.astype(str), "irr": irr[pitched]})
print("median kWh/m2/yr on pitched roofs by aspect: "
      + ", ".join(f"{k} {v:.0f}" for k, v in by_face.groupby("face", sort=False)["irr"].median()
                  .reindex(["N", "NE", "E", "SE", "S", "SW", "W", "NW"]).items())
      + f" | flat roofs {np.median(irr[~pitched]):.0f}")

# ---------------------------------------------------------------- Step 8: kWh per building and ranking
section("Step 8: kWh per year per building")
usable = (s <= MAX_SLOPE) & (irr >= MIN_IRR)
tilted_m2 = 1 / np.cos(np.radians(s))                      # a 1 m x 1 m cell on a 40 deg roof is 1.31 m2 of roof
kwh = np.where(usable, irr * tilted_m2 * COVERAGE * PANEL_EFF * PERF_RATIO, 0)
kwp = np.where(usable, tilted_m2 * COVERAGE * PANEL_EFF, 0)  # 1 m2 at 20% = 0.2 kWp at 1000 W/m2
cell_b = bid[rr, cc]
agg = pd.DataFrame({"b": cell_b, "roof_m2": tilted_m2, "usable_m2": np.where(usable, tilted_m2, 0),
                    "kwp": kwp, "kwh": kwh, "irr_w": np.where(usable, irr * tilted_m2, 0)}).groupby("b").sum()
bld["roof_m2"] = agg["roof_m2"].reindex(bld.index + 1).fillna(0).values
bld["usable_m2"] = agg["usable_m2"].reindex(bld.index + 1).fillna(0).values
bld["kwp"] = agg["kwp"].reindex(bld.index + 1).fillna(0).values
bld["kwh_yr"] = agg["kwh"].reindex(bld.index + 1).fillna(0).values
bld["mean_irr"] = (agg["irr_w"] / agg["usable_m2"].replace(0, np.nan)).reindex(bld.index + 1).values
print(f"assumptions: efficiency {PANEL_EFF:.0%}, performance ratio {PERF_RATIO}, coverage {COVERAGE:.0%}, "
      f"usable if slope <= {MAX_SLOPE} deg and >= {MIN_IRR} kWh/m2/yr")
print(f"usable share of roof area: {bld['usable_m2'].sum() / bld['roof_m2'].sum() * 100:.0f}%")
print(f"AOI total: {bld['kwp'].sum() / 1000:.1f} MWp, {bld['kwh_yr'].sum() / 1e6:.1f} GWh/yr = "
      f"{bld['kwh_yr'].sum() / HOUSEHOLD_KWH:,.0f} typical households ({HOUSEHOLD_KWH:,} kWh each)")
houses = bld[(bld["footprint_m2"] > 30) & (bld["footprint_m2"] < 120)]
print(f"house-sized footprints (30-120 m2): {len(houses):,} | median {houses['kwp'].median():.1f} kWp, "
      f"{houses['kwh_yr'].median():,.0f} kWh/yr | {(houses['kwh_yr'] >= HOUSEHOLD_KWH).mean() * 100:.0f}% "
      f"could cover {HOUSEHOLD_KWH:,} kWh/yr | {(houses['kwh_yr'] == 0).mean() * 100:.0f}% have no usable roof")
bld["rank"] = bld["kwh_yr"].rank(ascending=False, method="first").astype(int)
top = bld.sort_values("rank").head(10)
print("\ntop 10 roofs by annual output:")
print(top.assign(id=top["id"].str[:8], kwh_kwp=top["kwh_yr"] / top["kwp"],
                 type=top["class"].fillna(top["subtype"]).fillna("-"))[
    ["rank", "id", "type", "footprint_m2", "usable_m2", "kwp", "kwh_yr", "kwh_kwp", "mean_irr"]].astype({"kwp": float})
    .round({"kwp": 1, "footprint_m2": 0, "usable_m2": 0, "kwh_yr": 0, "kwh_kwp": 0, "mean_irr": 0})
    .astype({c: int for c in ["footprint_m2", "usable_m2", "kwh_yr", "kwh_kwp", "mean_irr"]})
    .to_string(index=False))
bld.drop(columns="geometry").sort_values("rank").to_csv(OUT / "roof_ranking.csv", index=False)
print(f"\nwrote {OUT / 'roof_ranking.csv'} ({len(bld):,} buildings)")

# ---------------------------------------------------------------- Figures
section("figures")
EXT = [0, W, 0, H]                                            # metres from the AOI south-west corner
dz_y, dz_x = np.gradient(dsm[A])
hill = np.clip(0.5 + 0.5 * (-dz_x * 0.7 + dz_y * 0.7) / np.sqrt(1 + dz_x ** 2 + dz_y ** 2), 0, 1)
off = lambda geo: geo.translate(-AOI[0], -AOI[1])
ZOOM = (330, 330 + 260, 330, 330 + 200)                      # x0, x1, y0, y1 in metres from the AOI SW corner


def put(values, fill=np.nan):
    """Scatter per-roof-cell values back onto the AOI grid."""
    g = np.full((H, W), fill, "float32"); g[rr, cc] = values
    return g


# 1: study area
fig, axs = plt.subplots(1, 2, figsize=(12, 5.2), layout="constrained")
for ax, zoom in zip(axs, (False, True)):
    ax.imshow(hill, cmap="gray", extent=EXT, vmin=0, vmax=1)
    im = ax.imshow(np.ma.masked_less(ndsm[A], 1), cmap="Spectral_r", vmin=1, vmax=15, extent=EXT, alpha=0.75)
    off(bld.geometry).boundary.plot(ax=ax, color="k", lw=0.25 if not zoom else 0.6)
    if zoom:
        ax.set_xlim(ZOOM[0], ZOOM[1]); ax.set_ylim(ZOOM[2], ZOOM[3])
    else:
        ax.add_patch(plt.Rectangle((ZOOM[0], ZOOM[2]), ZOOM[1] - ZOOM[0], ZOOM[3] - ZOOM[2], fill=False,
                                   ec="#3288bd", lw=1.5))
    ax.set_xlabel("metres east")
axs[0].set_ylabel("metres north")
axs[0].set_title(f"LiDAR height above ground and {len(bld):,} Overture footprints", loc="left")
axs[0].set_xlim(0, W); axs[0].set_ylim(0, H); axs[1].set_ylabel("")
axs[1].set_title("Zoom: terraces, back extensions and street trees", loc="left")
fig.colorbar(im, ax=axs, shrink=0.8, label="DSM - DTM (m)")
fig.savefig(OUT / "01-study-area-ndsm-footprints.png"); plt.close(fig)

# 2: slope and aspect
fig, axs = plt.subplots(1, 3, figsize=(13, 4.4), layout="constrained", gridspec_kw={"width_ratios": [1, 1, 0.8]})
for ax in axs[:2]:
    ax.imshow(hill, cmap="gray", extent=EXT, vmin=0, vmax=1)
    ax.set_xlim(ZOOM[0], ZOOM[1]); ax.set_ylim(ZOOM[2], ZOOM[3]); ax.set_xlabel("metres east")
im = axs[0].imshow(put(s), cmap="Spectral_r", vmin=0, vmax=60, extent=EXT, interpolation="nearest")
fig.colorbar(im, ax=axs[0], shrink=0.8, label="slope (degrees)")
axs[0].set_title("Roof slope", loc="left"); axs[0].set_ylabel("metres north")
asp_cls = put(np.where(pitched, ((a + 22.5) % 360 // 45), 8))
cm_asp = ListedColormap(["#3288bd", "#66c2a5", "#abdda4", "#fee08b", "#d53e4f", "#f46d43", "#fdae61", "#5e4fa2",
                         "#bdbdbd"])
im = axs[1].imshow(asp_cls, cmap=cm_asp, norm=BoundaryNorm(np.arange(-0.5, 9.5), 9), extent=EXT,
                   interpolation="nearest")
cb = fig.colorbar(im, ax=axs[1], shrink=0.8, ticks=range(9))
cb.ax.set_yticklabels(["N", "NE", "E", "SE", "S", "SW", "W", "NW", "flat"])
axs[1].set_title("Roof aspect (pitched roofs)", loc="left")
axs[2].hist(s[s <= 70], bins=np.arange(0, 71, 2), color="#3288bd")
axs[2].axvline(MAX_SLOPE, color="#d53e4f", ls="--", lw=1)
axs[2].set_xlabel("slope (degrees)"); axs[2].set_ylabel("roof cells (1 m2)")
axs[2].set_title("All roof cells in the AOI", loc="left")
fig.savefig(OUT / "02-roof-slope-aspect.png"); plt.close(fig)

# 3: shadows at noon on three dates
fig, axs = plt.subplots(1, 3, figsize=(13, 4.0), sharey=True, layout="constrained")
zr = slice(R0 + H - ZOOM[3], R0 + H - ZOOM[2]); zc = slice(C0 + ZOOM[0], C0 + ZOOM[1])
base = dsm[zr, zc]
for ax, day in zip(axs, ["2026-12-21", "2026-03-20", "2026-06-21"]):
    t = pd.Timestamp(f"{day} 12:10", tz="UTC")
    p = pvlib.solarposition.get_solarposition(pd.DatetimeIndex([t]), LAT, LON).iloc[0]
    tan_el = np.tan(np.radians(p["apparent_elevation"]))
    sa_, ca_ = np.sin(np.radians(p["azimuth"])), np.cos(np.radians(p["azimuth"]))
    shade = np.zeros(base.shape, bool)
    for d in np.arange(1, 301):
        dr, dc = -int(round(d * ca_)), int(round(d * sa_))
        shade |= dsm[zr.start + dr:zr.stop + dr, zc.start + dc:zc.stop + dc] - base - NOISE > d * tan_el
    ax.imshow(hill[H - ZOOM[3]:H - ZOOM[2], ZOOM[0]:ZOOM[1]], cmap="gray", vmin=0, vmax=1,
              extent=[ZOOM[0], ZOOM[1], ZOOM[2], ZOOM[3]])
    ax.imshow(np.ma.masked_where(~shade, shade), cmap=ListedColormap(["#3f007d"]), alpha=0.55,
              extent=[ZOOM[0], ZOOM[1], ZOOM[2], ZOOM[3]])
    off(bld.geometry).boundary.plot(ax=ax, color="#fee08b", lw=0.5)
    ax.set_xlim(ZOOM[0], ZOOM[1]); ax.set_ylim(ZOOM[2], ZOOM[3])
    ax.set_title(f"{pd.Timestamp(day):%d %b} 12:10 UTC: sun {p['apparent_elevation']:.0f} deg up, "
                 f"{shade.mean() * 100:.0f}% in shade", loc="left", fontsize=9)
    ax.set_xlabel("metres east")
    print(f"{day} 12:10 UTC: sun elevation {p['apparent_elevation']:.1f} deg, azimuth {p['azimuth']:.0f} deg, "
          f"{shade.mean() * 100:.0f}% of the zoom window in shadow")
axs[0].set_ylabel("metres north"); axs[1].set_ylabel(""); axs[2].set_ylabel("")
fig.savefig(OUT / "03-shadows-solstices-equinox.png"); plt.close(fig)

# 4: model check
months = range(1, 13)
mon = pd.DataFrame({"pvlib isotropic, PVGIS TMY": poa.groupby(poa.index.month).sum().values / 1000,
                    "pvlib Perez, PVGIS TMY": poa_perez.groupby(poa_perez.index.month).sum().values / 1000,
                    f"PVGIS SARAH {hourly.index.year.min()}-{hourly.index.year.max()} mean":
                        hourly["poa"].groupby(hourly.index.month).sum().values / 1000 / n_years,
                    "Ineichen clear sky": cs_poa.groupby(cs_poa.index.month).sum().values / 1000}, index=months)
fig, axs = plt.subplots(1, 2, figsize=(12, 4.0), layout="constrained", gridspec_kw={"width_ratios": [1.6, 1]})
w = 0.27
cols = ["#3288bd", "#66c2a5", "#fdae61"]
for i, (lab, col) in enumerate(zip(mon.columns[:3], cols)):
    axs[0].bar(np.array(months) + (i - 1) * w, mon[lab], w, color=col, label=f"{lab} ({mon[lab].sum():.0f})")
axs[0].plot(months, mon["Ineichen clear sky"], "k--", lw=1,
            label=f"Ineichen clear sky ({mon['Ineichen clear sky'].sum():.0f})")
axs[0].set_xticks(months, ["J", "F", "M", "A", "M", "J", "J", "A", "S", "O", "N", "D"])
axs[0].set_ylabel("kWh/m2 per month"); axs[0].legend(fontsize=7, title="annual kWh/m2", title_fontsize=7)
axs[0].set_title("35 deg south-facing plane in Levenshulme, no shading", loc="left")
grid = [[float(sum(annual_irradiance(np.array([tilt]), np.array([az]), np.zeros((72, 1)), np.ones(1)))[0])
         for az in np.arange(90, 271, 15)] for tilt in np.arange(0, 91, 10)]
im = axs[1].imshow(grid, origin="lower", cmap="Spectral_r", aspect="auto", extent=[82.5, 277.5, -5, 95])
cs_ = axs[1].contour(np.arange(90, 271, 15), np.arange(0, 91, 10), grid, levels=[700, 800, 900, 1000], colors="k",
                     linewidths=0.6)
axs[1].clabel(cs_, fontsize=7, fmt="%d")
fig.colorbar(im, ax=axs[1], label="kWh/m2/yr")
axs[1].set_xticks([90, 135, 180, 225, 270], ["E", "SE", "S", "SW", "W"])
axs[1].set_xlabel("roof aspect"); axs[1].set_ylabel("roof slope (degrees)")
axs[1].set_title("Unshaded annual irradiation by slope and aspect", loc="left")
fig.savefig(OUT / "04-irradiance-model-check.png"); plt.close(fig)
best = np.unravel_index(np.argmax(grid), np.shape(grid))
print(f"best unshaded orientation on the grid: slope {best[0] * 10} deg, aspect {90 + best[1] * 15} deg, "
      f"{np.max(grid):.0f} kWh/m2/yr | flat {grid[0][0]:.0f} | vertical south wall {grid[9][6]:.0f}")

# 5: annual irradiation map
fig, axs = plt.subplots(1, 2, figsize=(12, 5.2), layout="constrained")
g_irr = put(irr)
for ax, zoom in zip(axs, (False, True)):
    ax.imshow(hill, cmap="gray", extent=EXT, vmin=0, vmax=1)
    im = ax.imshow(g_irr, cmap="Spectral_r", vmin=300, vmax=1100, extent=EXT, interpolation="nearest")
    if zoom:
        ax.set_xlim(ZOOM[0], ZOOM[1]); ax.set_ylim(ZOOM[2], ZOOM[3])
        off(bld.geometry).boundary.plot(ax=ax, color="k", lw=0.4)
    ax.set_xlabel("metres east")
axs[0].set_ylabel("metres north")
axs[0].set_title("Annual irradiation on roofs, with shading (PVGIS TMY)", loc="left")
axs[1].set_title("Zoom: south-west and north-east pitches of the same terraces", loc="left")
fig.colorbar(im, ax=axs, shrink=0.8, label="kWh/m2/yr on the roof surface")
axs[1].set_ylabel("")
fig.savefig(OUT / "05-annual-irradiation-roofs.png"); plt.close(fig)

# 6: building ranking
fig, axs = plt.subplots(1, 2, figsize=(12, 5.0), layout="constrained", gridspec_kw={"width_ratios": [1.3, 1]})
axs[0].imshow(hill, cmap="gray", extent=EXT, vmin=0, vmax=1)
kcm = ListedColormap(plt.get_cmap("Spectral")(np.linspace(0.05, 0.95, 6)))
knorm = BoundaryNorm(np.arange(-0.5, 6.5), kcm.N)
b_off = bld.assign(geometry=off(bld.geometry),
                   kclass=np.digitize(bld["kwh_yr"], [1, 1500, HOUSEHOLD_KWH, 5000, 20000]))
b_off.plot(ax=axs[0], column="kclass", cmap=kcm.reversed(), norm=knorm, lw=0)
sm = plt.cm.ScalarMappable(cmap=kcm.reversed(), norm=knorm)
cb = fig.colorbar(sm, ax=axs[0], shrink=0.8, ticks=range(6))
cb.ax.set_yticklabels(["0", "1-1,500", "1,500-2,700", "2,700-5,000", "5,000-20,000", "20,000+"])
cb.set_label("estimated kWh per year")
for _, r in b_off[b_off["rank"] <= 5].iterrows():
    p = r.geometry.representative_point()
    axs[0].annotate(str(r["rank"]), (p.x, p.y), fontsize=8, ha="center", va="center", weight="bold",
                    bbox=dict(boxstyle="circle,pad=0.15", fc="white", ec="k", lw=0.5))
axs[0].set_xlim(0, W); axs[0].set_ylim(0, H); axs[0].set_xlabel("metres east"); axs[0].set_ylabel("metres north")
axs[0].set_title("Estimated PV output per building (top 5 numbered)", loc="left")
hk = houses["kwh_yr"]
axs[1].hist(hk, bins=np.arange(0, 8001, 250), color="#3288bd")
axs[1].axvline(HOUSEHOLD_KWH, color="#d53e4f", ls="--", lw=1.2)
axs[1].annotate(f"typical household\n{HOUSEHOLD_KWH:,} kWh/yr", (HOUSEHOLD_KWH, axs[1].get_ylim()[1] * 0.85),
                xytext=(6, 0), textcoords="offset points", color="#d53e4f", fontsize=8,
                bbox=dict(fc="white", ec="none", alpha=0.85))
axs[1].set_xlabel("estimated kWh per year"); axs[1].set_ylabel("buildings")
axs[1].set_title(f"House-sized roofs (30-120 m2, n = {len(houses):,})", loc="left")
fig.savefig(OUT / "06-roof-ranking-kwh.png"); plt.close(fig)
print("figures written to", OUT.resolve())