Source code: Land surface temperature from Landsat in Python: mapping London's urban heat island

"""Urban heat island mapping with Landsat land surface temperature in Python: London, 18 July 2022.

Streams a Landsat 9 Collection 2 Level-2 scene from the Microsoft Planetary Computer, converts the
ST_B10 surface temperature band to degrees Celsius, masks cloud with QA_PIXEL, and then measures
the surface urban heat island: LST by ESA WorldCover land cover and by NDVI, urban against rural
reference land, a west-east transect through central London, the hottest clusters overlaid on
OpenStreetMap, and a check against air temperatures measured at weather stations.
Run with:  python lst_urban_heat_island.py
"""
import json
import platform
import time
from importlib.metadata import version
from pathlib import Path

import geopandas as gpd
import matplotlib

matplotlib.use("Agg")
import matplotlib.patheffects as pe
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import numpy as np
import osmnx as ox
import pandas as pd
import planetary_computer
import pystac_client
import rasterio
import rasterio.features
import requests
from pyproj import Transformer
from rasterio.enums import Resampling
from rasterio.vrt import WarpedVRT
from rasterio.windows import from_bounds
from scipy import ndimage

HERE = Path(__file__).parent
OUT, DATA, CACHE = HERE / "output", HERE / "data", HERE / "cache"
for d in (OUT, DATA, CACHE):
    d.mkdir(exist_ok=True)
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9})
pd.set_option("display.width", 140)
HEADERS = {"User-Agent": "TechnicalGIS tutorial (https://www.technicalgis.com)"}

print("=== versions ===")
print(f"python {platform.python_version()} | pystac-client {pystac_client.__version__} | "
      f"planetary-computer {version('planetary-computer')} | rasterio {rasterio.__version__} | "
      f"osmnx {ox.__version__} | numpy {np.__version__} | pandas {pd.__version__}")


def section(title: str) -> None:
    print(f"\n=== {title} ===")


# --- Step 1: find a clear Landsat 8/9 scene on a hot day ------------------------------------
STAC_URL = "https://planetarycomputer.microsoft.com/api/stac/v1"
AOI = (-0.95, 51.15, 0.65, 51.85)    # lon/lat: Greater London plus a rural margin on every side
CENTRE = (-0.1276, 51.5072)           # Trafalgar Square / Charing Cross, the traditional centre of London

catalog = pystac_client.Client.open(STAC_URL, modifier=planetary_computer.sign_inplace)
search = catalog.search(collections=["landsat-c2-l2"], bbox=AOI, datetime="2022-06-01/2022-08-31",
                        query={"platform": {"in": ["landsat-8", "landsat-9"]}, "eo:cloud_cover": {"lt": 10}})
summer = sorted(search.items(), key=lambda i: i.datetime)
section("STAC search: clear Landsat 8/9 scenes over London, summer 2022")
for it in summer:
    p = it.properties
    print(f"{it.datetime:%Y-%m-%d %H:%M} UTC  {p['platform']}  path/row {p['landsat:wrs_path']}/{p['landsat:wrs_row']}"
          f"  cloud {p['eo:cloud_cover']:5.2f}%")
item = next(it for it in summer if it.id == "LC09_L2SP_201024_20220718_02_T1")
p = item.properties
print(f"\nusing {item.id}: {item.datetime:%d %B %Y %H:%M} UTC, sun elevation {p['view:sun_elevation']:.1f} deg, "
      f"{p['proj:code']}")
st = item.assets["lwir11"].extra_fields["raster:bands"][0]
print(f"ST_B10 (asset 'lwir11'): scale {st['scale']}, offset {st['offset']}, unit {st['unit']}, nodata {st['nodata']}")

# --- Step 2: read the bands, apply scale and offset, convert to degrees Celsius --------------------
CRS = p["proj:code"]                  # EPSG:32631, the scene's native UTM zone
to_utm = Transformer.from_crs("EPSG:4326", CRS, always_xy=True)
xs, ys = to_utm.transform([AOI[0], AOI[2], AOI[0], AOI[2]], [AOI[1], AOI[1], AOI[3], AOI[3]])
BOUNDS = (min(xs), min(ys), max(xs), max(ys))
BANDS = ["lwir11", "qa_pixel", "qa", "red", "green", "blue", "nir08"]

raw_file = CACHE / f"{item.id}.npz"
if not raw_file.exists():
    raw = {}
    for b in BANDS:
        with rasterio.open(item.assets[b].href) as src:
            win = from_bounds(*BOUNDS, src.transform).round_offsets().round_lengths()
            raw[b] = src.read(1, window=win)
            transform = src.window_transform(win)
    raw["transform"] = np.array(transform)[:6]
    np.savez_compressed(raw_file, **raw)
raw = dict(np.load(raw_file))
TRANSFORM = rasterio.Affine(*raw.pop("transform"))
HEIGHT, WIDTH = raw["lwir11"].shape
x0, y1 = TRANSFORM.c, TRANSFORM.f
x1, y0 = x0 + WIDTH * TRANSFORM.a, y1 + HEIGHT * TRANSFORM.e
EXTENT = (x0, x1, y0, y1)
PX_KM2 = 30 * 30 / 1e6

dn = raw["lwir11"]
lst_k = np.where(dn > 0, dn * 0.00341802 + 149.0, np.nan)       # Collection 2 Level-2 ST scale and offset
lst = lst_k - 273.15                                             # kelvin to degrees Celsius
st_unc = np.where(raw["qa"] >= 0, raw["qa"] * 0.01, np.nan)     # ST_QA: uncertainty in kelvin


def sr(band: str) -> np.ndarray:
    """Collection 2 surface reflectance: DN * 0.0000275 - 0.2."""
    v = raw[band]
    return np.where(v > 0, v * 2.75e-5 - 0.2, np.nan)


red, nir = sr("red"), sr("nir08")
with np.errstate(divide="ignore", invalid="ignore"):
    ndvi = (nir - red) / (nir + red)
section("surface temperature")
print(f"grid {WIDTH} x {HEIGHT} pixels of 30 m ({(x1 - x0) / 1000:.0f} x {(y1 - y0) / 1000:.0f} km)")
print(f"DN {np.nanmin(dn[dn > 0])}..{dn.max()} -> {np.nanmin(lst_k):.1f}..{np.nanmax(lst_k):.1f} K -> "
      f"{np.nanmin(lst):.1f}..{np.nanmax(lst):.1f} degC")

# --- Step 3: mask cloud, shadow and water with QA_PIXEL ----------------------------------------------
qa = raw["qa_pixel"]
bit = lambda b: (qa >> b) & 1 == 1   # noqa: E731
QA_BITS = {"fill": 0, "dilated cloud": 1, "cirrus": 2, "cloud": 3, "cloud shadow": 4, "snow": 5, "clear": 6, "water": 7}
for name, b in QA_BITS.items():
    print(f"QA_PIXEL bit {b} ({name}): {bit(b).mean():.2%}")
bad = bit(0) | bit(1) | bit(2) | bit(3) | bit(4) | bit(5)
bad = ndimage.binary_dilation(bad, iterations=3)                  # 90 m buffer: thermal pixels are really 100 m
water = bit(7)
ok = ~bad & ~water & np.isfinite(lst)
print(f"usable land pixels: {ok.mean():.1%} | ST_QA uncertainty over land: median {np.nanmedian(st_unc[ok]):.2f} K, "
      f"95th pct {np.nanpercentile(st_unc[ok], 95):.2f} K")
with rasterio.open(OUT / "lst_london_2022-07-18.tif", "w", driver="GTiff", width=WIDTH, height=HEIGHT, count=1,
                   dtype="float32", crs=CRS, transform=TRANSFORM, nodata=np.nan, compress="deflate") as dst:
    dst.write(np.where(ok, lst, np.nan).astype("float32"), 1)

# --- Step 4: Greater London boundary and ESA WorldCover land cover ---------------------------------------
gl_file = DATA / "greater_london.gpkg"
if not gl_file.exists():
    ox.settings.cache_folder = CACHE / "osmnx"
    ox.geocode_to_gdf("Greater London, England")[["geometry"]].to_file(gl_file, driver="GPKG")
gl = gpd.read_file(gl_file).to_crs(CRS)
in_gl = rasterio.features.geometry_mask(gl.geometry, (HEIGHT, WIDTH), TRANSFORM, invert=True)

wc_file = CACHE / "worldcover_2021_30m.npy"
if not wc_file.exists():
    wc = np.zeros((HEIGHT, WIDTH), dtype="uint8")
    for wi in catalog.search(collections=["esa-worldcover"], bbox=AOI,
                             query={"esa_worldcover:product_version": {"eq": "2.0.0"}}).items():
        with rasterio.open(wi.assets["map"].href) as src, WarpedVRT(src, crs=CRS, transform=TRANSFORM, width=WIDTH,
                                                                   height=HEIGHT, resampling=Resampling.mode) as vrt:
            tile = vrt.read(1)
        wc = np.where(tile > 0, tile, wc)
    np.save(wc_file, wc)
wc = np.load(wc_file)
WC = {10: "Tree cover", 20: "Shrubland", 30: "Grassland", 40: "Cropland", 50: "Built-up", 60: "Bare / sparse",
      80: "Water", 90: "Wetland"}

section("LST by land cover inside Greater London (ESA WorldCover 2021)")
rows = []
for code, name in WC.items():
    m = ok & in_gl & (wc == code)
    if m.sum() * PX_KM2 < 5:          # skip classes covering less than 5 km2 of London
        continue
    v = lst[m]
    rows.append({"land cover": name, "km2": round(m.sum() * PX_KM2, 1), "median_C": np.median(v),
                 "p10_C": np.percentile(v, 10), "p90_C": np.percentile(v, 90), "median_NDVI": np.median(ndvi[m])})
by_lc = pd.DataFrame(rows).sort_values("median_C", ascending=False).round(2)
print(by_lc.to_string(index=False))
by_lc.to_csv(OUT / "lst_by_landcover.csv", index=False)

# --- Step 5: urban against rural: the surface urban heat island intensity ----------------------------------
built_1km = ndimage.uniform_filter((wc == 50).astype("float32"), size=33)          # built-up share within ~1 km
dist_gl_km = ndimage.distance_transform_edt(~in_gl) * 0.03
rural = ok & ~in_gl & (dist_gl_km > 5) & np.isin(wc, [30, 40]) & (built_1km < 0.05)
urban = ok & in_gl & (wc == 50)
suhi = np.median(lst[urban]) - np.median(lst[rural])
section("surface urban heat island")
print(f"urban: built-up pixels in Greater London, {urban.sum() * PX_KM2:,.0f} km2, median LST {np.median(lst[urban]):.1f} C")
print(f"rural: grass/cropland > 5 km outside London with < 5% built-up within 1 km, {rural.sum() * PX_KM2:,.0f} km2, "
      f"median LST {np.median(lst[rural]):.1f} C")
print(f"surface urban heat island intensity (median urban - median rural): {suhi:.1f} C")
inner = ok & in_gl & (np.hypot(*np.meshgrid(np.arange(WIDTH) * 30 + x0 + 15 - to_utm.transform(*CENTRE)[0],
                                             y1 - 15 - np.arange(HEIGHT)[:, None] * 30 - to_utm.transform(*CENTRE)[1],
                                             sparse=True)) < 5000)
print(f"all land within 5 km of Charing Cross: median LST {np.median(lst[inner]):.1f} C")

gl_ok = ok & in_gl
bins = np.arange(-0.1, 0.95, 0.1)
idx = np.digitize(ndvi[gl_ok], bins)
ndvi_tab = pd.DataFrame({"ndvi": ndvi[gl_ok], "lst": lst[gl_ok], "bin": idx})
ndvi_tab = ndvi_tab[(ndvi_tab["bin"] > 0) & (ndvi_tab["bin"] < len(bins))]
summary = ndvi_tab.groupby("bin").agg(ndvi_from=("ndvi", "min"), n=("lst", "size"), median_C=("lst", "median")).round(2)
summary["ndvi_from"] = bins[summary.index - 1].round(1)
slope, intercept = np.polyfit(ndvi[gl_ok], lst[gl_ok], 1)
r = np.corrcoef(ndvi[gl_ok], lst[gl_ok])[0, 1]
print("\nLST by NDVI band inside Greater London:")
print(summary.to_string())
print(f"linear fit: LST = {intercept:.1f} {slope:+.1f} x NDVI  (r = {r:.2f}, {slope / 10:+.1f} C per 0.1 NDVI)")

# --- Step 6: a west-east transect through central London ------------------------------------------------------
section("west-east transect along 51.507 N")
cx, cy = to_utm.transform(*CENTRE)
row = int((y1 - cy) // 30)
tx = x0 + 15 + np.arange(WIDTH) * 30
t_lst = pd.Series(np.where(ok[row], lst[row], np.nan))
t_ndvi = pd.Series(np.where(ok[row], ndvi[row], np.nan))
t_lst_s = t_lst.rolling(31, center=True, min_periods=10).median()            # ~1 km running median
t_ndvi_s = t_ndvi.rolling(31, center=True, min_periods=10).median()
km = (tx - cx) / 1000
inside_row = in_gl[row]
edges_gl = km[np.flatnonzero(np.diff(inside_row.astype(int)))]
print(f"transect length {km[-1] - km[0]:.0f} km; Greater London boundary crossed at {', '.join(f'{e:+.1f}' for e in edges_gl)} km")
for lo, hi, name in ((-60, -40, "rural west (Berkshire)"), (-5, 5, "central London"), (40, 60, "rural east (Essex)")):
    sel = (km >= lo) & (km < hi)
    print(f"{name:<24} {lo:+} to {hi:+} km: median LST {np.nanmedian(t_lst[sel]):.1f} C, NDVI {np.nanmedian(t_ndvi[sel]):.2f}")

# --- Step 7: hotspots, named with OpenStreetMap ---------------------------------------------------------------------
section("hotspots (top 2% of Greater London land pixels)")
thr = np.percentile(lst[gl_ok], 98)
hot = gl_ok & (lst >= thr)
hot = ndimage.binary_closing(hot, iterations=1) & gl_ok
lab, n = ndimage.label(hot)
sizes = ndimage.sum(hot, lab, index=np.arange(1, n + 1))
print(f"threshold {thr:.1f} C; {n:,} patches, {np.sum(sizes * PX_KM2 * 100 >= 5)} of at least 5 ha")
order = np.argsort(sizes)[::-1][:10] + 1
to_ll = Transformer.from_crs(CRS, "EPSG:4326", always_xy=True)
geo_cache = CACHE / "nominatim.json"
names = json.loads(geo_cache.read_text()) if geo_cache.exists() else {}
hs = []
for k, lid in enumerate(order, 1):
    m = lab == lid
    rr, cc = np.nonzero(m)
    j = np.argmax(lst[rr, cc])
    hx, hy = x0 + (cc[j] + 0.5) * 30, y1 - (rr[j] + 0.5) * 30
    lon, lat = to_ll.transform(hx, hy)
    key = f"{lat:.5f},{lon:.5f}"
    if key not in names:
        resp = requests.get("https://nominatim.openstreetmap.org/reverse", headers=HEADERS, timeout=30,
                            params={"lat": lat, "lon": lon, "format": "jsonv2", "zoom": 17})
        names[key] = resp.json()
        time.sleep(1.1)                                    # Nominatim usage policy: at most one request a second
    a = names[key].get("address", {})
    place = next((a[f] for f in ("suburb", "town", "village", "neighbourhood") if f in a), "")
    feature = next((a[f] for f in ("industrial", "retail", "commercial", "leisure", "park", "amenity") if f in a), "")
    hs.append({"rank": k, "area_ha": round(m.sum() * PX_KM2 * 100, 1), "max_C": round(float(lst[m].max()), 1),
               "mean_C": round(float(lst[m].mean()), 1), "place": place,
               "borough": a.get("city_district", "").replace("London Borough of ", ""),
               "near": feature or a.get("road", ""), "built_%": round(100 * float((wc[m] == 50).mean())),
               "grass_%": round(100 * float(np.isin(wc[m], [30, 40]).mean())),
               "tree_%": round(100 * float((wc[m] == 10).mean())), "NDVI": round(float(np.median(ndvi[m])), 2),
               "x": hx, "y": hy, "lat": round(lat, 5), "lon": round(lon, 5)})
geo_cache.write_text(json.dumps(names))
hs = pd.DataFrame(hs)
print(hs.drop(columns=["x", "y", "lat", "lon"]).to_string(index=False))
hs.drop(columns=["x", "y"]).to_csv(OUT / "hotspots.csv", index=False)
all_hot = hot & (np.isin(lab, order))
print(f"ten largest hotspots together: {all_hot.sum() * PX_KM2 * 100:,.0f} ha, "
      f"{100 * (wc[all_hot] == 50).mean():.0f}% built-up, {100 * np.isin(wc[all_hot], [30, 40]).mean():.0f}% grass or crops, "
      f"{100 * (wc[all_hot] == 10).mean():.0f}% trees (WorldCover)")

# OpenStreetMap context for a zoomed map around the largest hotspot.
zx, zy = hs.loc[0, "x"], hs.loc[0, "y"]
ZOOM = (zx - 4000, zy - 3000, zx + 4000, zy + 3000)
zoom_ll = to_ll.transform([ZOOM[0], ZOOM[2]], [ZOOM[1], ZOOM[3]])
osm_file = DATA / "osm_zoom.gpkg"
if not osm_file.exists():
    ox.settings.cache_folder = CACHE / "osmnx"
    tags = {"landuse": ["residential", "industrial", "commercial", "retail", "railway"], "leisure": ["park", "golf_course"],
            "highway": ["motorway", "trunk", "primary"]}
    osm = None
    for attempt in range(4):              # public Overpass servers are often busy: retry with a mirror
        for url in ("https://overpass-api.de/api", "https://overpass.kumi.systems/api"):
            ox.settings.overpass_url = url
            try:
                osm = ox.features_from_bbox((zoom_ll[0][0], zoom_ll[1][0], zoom_ll[0][1], zoom_ll[1][1]), tags)
                break
            except Exception as exc:
                print(f"{url} failed: {str(exc).splitlines()[0][:70]}")
        if osm is not None:
            break
        time.sleep(30 * (attempt + 1))
    keep = [c for c in ("landuse", "leisure", "highway", "name") if c in osm.columns]
    osm = osm.reset_index()[keep + ["geometry"]]
    osm.to_file(osm_file, driver="GPKG")
osm = gpd.read_file(osm_file).to_crs(CRS)
zr = (slice(int((y1 - ZOOM[3]) // 30), int((y1 - ZOOM[1]) // 30)), slice(int((ZOOM[0] - x0) // 30), int((ZOOM[2] - x0) // 30)))
zone_ok = ok[zr]
zone_hot = hot[zr]
print(f"\nOpenStreetMap features in the 8 x 6 km box around hotspot 1: {len(osm):,}")
print(f"{'OSM tag':<22} {'km2':>5} {'median LST':>11} {'share of hotspot pixels':>24}")
for col, vals in (("landuse", ["residential", "industrial", "commercial", "retail"]), ("leisure", ["park", "golf_course"])):
    if col not in osm.columns:
        continue
    for v in vals:
        polys = osm[(osm[col] == v) & osm.geom_type.isin(["Polygon", "MultiPolygon"])]
        if len(polys):
            mm = rasterio.features.geometry_mask(polys.geometry, (HEIGHT, WIDTH), TRANSFORM, invert=True)[zr] & zone_ok
            if mm.sum() > 50:
                print(f"{col + '=' + v:<22} {mm.sum() * PX_KM2:5.2f} {np.median(lst[zr][mm]):9.1f} C "
                      f"{100 * (mm & zone_hot).sum() / zone_hot.sum():23.0f}%")

# --- Step 8: land surface temperature is not air temperature ------------------------------------------------------
section("LST against air temperature at weather stations (Meteostat hourly, 11:00 UTC)")
STATIONS = {"03772": ("Heathrow", 51.4833, -0.45), "03672": ("Northolt", 51.55, -0.4167),
            "03781": ("Kenley", 51.30, -0.0833), "EGKB0": ("Biggin Hill", 51.3167, 0.0333)}
met = []
for sid, (name, lat, lon) in STATIONS.items():
    f = CACHE / f"meteostat_{sid}_2022.csv.gz"
    if not f.exists():
        f.write_bytes(requests.get(f"https://bulk.meteostat.net/v2/hourly/2022/{sid}.csv.gz", headers=HEADERS,
                                   timeout=60).content)
    obs = pd.read_csv(f, header=None, usecols=[0, 1, 2], names=["date", "hour", "temp"])
    t = obs[(obs["date"] == "2022-07-18") & (obs["hour"].isin([10, 11]))]["temp"]
    sx, sy = to_utm.transform(lon, lat)
    rr, cc = int((y1 - sy) // 30), int((sx - x0) // 30)
    win = (slice(rr - 16, rr + 17), slice(cc - 16, cc + 17))      # ~1 km box: station coordinates are rounded
    v = lst[win][ok[win]]
    grass = lst[win][ok[win] & np.isin(wc[win], [30, 40])]
    met.append({"station": name, "air_10UTC_C": t.iloc[0], "air_11UTC_C": t.iloc[1],
                "LST_median_1km_C": round(float(np.median(v)), 1),
                "LST_grass_median_C": round(float(np.median(grass)), 1) if grass.size > 20 else np.nan})
met = pd.DataFrame(met)
met["LST_minus_air"] = (met["LST_median_1km_C"] - met["air_11UTC_C"]).round(1)
print(met.to_string(index=False))

# --- Figures --------------------------------------------------------------------------------------------------------
LST_CMAP = plt.get_cmap("inferno").copy()
LST_CMAP.set_bad("#bfc7cc")
lst_show = np.where(ok, lst, np.nan)
vmin, vmax = np.nanpercentile(lst_show, [1, 99.5])
halo = [pe.withStroke(linewidth=2.5, foreground="white")]


def frame(ax, xlim=(x0, x1), ylim=(y0, y1), bar_km=10):
    ax.set_xticks([]), ax.set_yticks([])
    ax.set_xlim(*xlim), ax.set_ylim(*ylim)
    xr, yb = xlim[1] - (xlim[1] - xlim[0]) * 0.04, ylim[0] + (ylim[1] - ylim[0]) * 0.05
    ax.plot([xr - bar_km * 1000, xr], [yb, yb], color="black", lw=3, solid_capstyle="butt")
    ax.text(xr - bar_km * 500, yb + (ylim[1] - ylim[0]) * 0.015, f"{bar_km} km", ha="center", va="bottom",
            fontsize=8, path_effects=halo)


true = np.dstack([sr("red"), sr("green"), sr("blue")])
true = np.clip(np.nan_to_num(true) * 3.2, 0, 1) ** (1 / 1.3)
fig, axes = plt.subplots(1, 2, figsize=(13, 4.6), constrained_layout=True)
axes[0].imshow(true, extent=EXTENT)
gl.boundary.plot(ax=axes[0], color="white", lw=0.8)
frame(axes[0])
axes[0].set_title(f"Landsat 9 true colour, {item.datetime:%d %B %Y %H:%M} UTC", fontsize=9)
im = axes[1].imshow(lst_show, cmap=LST_CMAP, vmin=vmin, vmax=vmax, extent=EXTENT, interpolation="nearest")
gl.boundary.plot(ax=axes[1], color="#2a9df4", lw=0.8)
frame(axes[1])
axes[1].set_title("Land surface temperature (ST_B10), grey = water or masked", fontsize=9)
fig.colorbar(im, ax=axes[1], shrink=0.85, label="LST (°C)")
fig.savefig(OUT / "01-true-colour-lst.png", dpi=110)
plt.close(fig)

fig, axes = plt.subplots(1, 2, figsize=(13, 4.4), constrained_layout=True, gridspec_kw={"width_ratios": [1.1, 1]})
order_lc = by_lc["land cover"].tolist()
code_of = {v: k for k, v in WC.items()}
data = [lst[ok & in_gl & (wc == code_of[n])] for n in order_lc]
lc_cols = {"Built-up": "#c92a2a", "Bare / sparse": "#a0785a", "Cropland": "#e8b923", "Grassland": "#94c973",
           "Shrubland": "#6c9a3d", "Tree cover": "#2b7a3b", "Wetland": "#4aa3a2", "Water": "#1c7ed6"}
bp = axes[0].boxplot(data, vert=False, showfliers=False, patch_artist=True, widths=0.6,
                     tick_labels=[f"{n} ({by_lc.loc[by_lc['land cover'] == n, 'km2'].iloc[0]:,.0f} km²)" for n in order_lc])
for patch, n in zip(bp["boxes"], order_lc):
    patch.set_facecolor(lc_cols.get(n, "grey"))
    patch.set_alpha(0.85)
for med in bp["medians"]:
    med.set_color("black")
axes[0].invert_yaxis()
axes[0].set_xlabel("LST (°C), Greater London, 18 July 2022")
axes[0].set_title("LST by ESA WorldCover class (box = 25th–75th percentile)", fontsize=9)
axes[0].spines[["top", "right"]].set_visible(False)
hb = axes[1].hexbin(ndvi[gl_ok], lst[gl_ok], gridsize=60, cmap="Blues", mincnt=5, bins="log", extent=(-0.2, 0.9, 20, 55))
xx = np.linspace(-0.1, 0.85, 10)
axes[1].plot(xx, intercept + slope * xx, color="#c92a2a", lw=1.5,
             label=f"fit: {slope / 10:+.1f} °C per 0.1 NDVI (r = {r:.2f})")
axes[1].set_xlabel("NDVI")
axes[1].set_ylabel("LST (°C)")
axes[1].set_title("Greener pixels are cooler (Greater London)", fontsize=9)
axes[1].legend(loc="upper right", fontsize=8)
fig.colorbar(hb, ax=axes[1], shrink=0.85, label="pixels (log scale)")
fig.savefig(OUT / "02-lst-landcover-ndvi.png")
plt.close(fig)

fig, ax = plt.subplots(figsize=(11, 4.4), constrained_layout=True)
ax.plot(km, t_lst, color="#f59f00", lw=0.4, alpha=0.6, label="LST, 30 m pixels")
ax.plot(km, t_lst_s, color="#c92a2a", lw=1.8, label="LST, 1 km running median")
for e in edges_gl:
    ax.axvline(e, color="grey", lw=0.8, ls="--")
if len(edges_gl):
    ax.axvspan(edges_gl[0], edges_gl[-1], color="grey", alpha=0.08, lw=0)
    ax.text((edges_gl[0] + edges_gl[-1]) / 2, 0.97, "Greater London", transform=ax.get_xaxis_transform(),
            ha="center", va="top", fontsize=8)
ax.text(0, 0.03, "Charing Cross", transform=ax.get_xaxis_transform(), ha="center", fontsize=8)
ax.set_xlabel("km east of Charing Cross along latitude 51.507° N")
ax.set_ylabel("LST (°C)")
ax2 = ax.twinx()
ax2.plot(km, t_ndvi_s, color="#2b8a3e", lw=1.3, label="NDVI, 1 km running median")
ax2.set_ylabel("NDVI", color="#2b8a3e")
ax2.set_ylim(-0.2, 1.4)
h1, l1 = ax.get_legend_handles_labels()
h2, l2 = ax2.get_legend_handles_labels()
ax.legend(h1 + h2, l1 + l2, loc="upper left", fontsize=8, framealpha=0.9)
ax.spines["top"].set_visible(False)
ax2.spines["top"].set_visible(False)
fig.savefig(OUT / "03-urban-rural-transect.png")
plt.close(fig)

fig, axes = plt.subplots(1, 2, figsize=(13, 5.0), constrained_layout=True, gridspec_kw={"width_ratios": [1.25, 1]})
gx0, gy0, gx1, gy1 = gl.total_bounds
axes[0].imshow(lst_show, cmap=LST_CMAP, vmin=vmin, vmax=vmax, extent=EXTENT, interpolation="nearest")
axes[0].imshow(np.ma.masked_where(~hot, hot), cmap="cool", extent=EXTENT, interpolation="nearest", alpha=1)
gl.boundary.plot(ax=axes[0], color="#2a9df4", lw=0.8)
for _, h in hs.iterrows():
    axes[0].annotate(str(h["rank"]), (h["x"], h["y"]), xytext=(5, 5), textcoords="offset points", fontsize=8,
                     fontweight="bold", path_effects=halo)
axes[0].add_patch(plt.Rectangle((ZOOM[0], ZOOM[1]), ZOOM[2] - ZOOM[0], ZOOM[3] - ZOOM[1], fill=False, ec="#2a9df4", lw=1.2))
frame(axes[0], (gx0 - 2000, gx1 + 2000), (gy0 - 2000, gy1 + 2000))
axes[0].set_title(f"Hotspots: top 2% of Greater London land pixels (≥ {thr:.1f} °C, cyan), ten largest numbered",
                  fontsize=9)
ax = axes[1]
ax.imshow(lst_show, cmap=LST_CMAP, vmin=vmin, vmax=vmax, extent=EXTENT, interpolation="nearest", alpha=0.9)
polys = osm[osm.geom_type.isin(["Polygon", "MultiPolygon"])]
if "landuse" in polys.columns:
    polys[polys["landuse"].isin(["industrial", "commercial", "retail"])].boundary.plot(ax=ax, color="white", lw=0.8)
if "leisure" in polys.columns:
    polys[polys["leisure"].isin(["park", "golf_course"])].boundary.plot(ax=ax, color="#69db7c", lw=1.0)
if "highway" in osm.columns:
    osm[osm["highway"].isin(["motorway", "trunk", "primary"]) & osm.geom_type.isin(["LineString", "MultiLineString"])
        ].plot(ax=ax, color="#4dabf7", lw=0.9)
ax.contour(hot, levels=[0.5], extent=(x0, x1, y1, y0), colors="cyan", linewidths=0.9)
frame(ax, (ZOOM[0], ZOOM[2]), (ZOOM[1], ZOOM[3]), bar_km=1)
ax.legend(handles=[Line2D([], [], color="cyan", label=f"hotspot (≥ {thr:.1f} °C)"),
                   Line2D([], [], color="white", label="OSM industrial / retail"),
                   Line2D([], [], color="#69db7c", label="OSM park / golf"),
                   Line2D([], [], color="#4dabf7", label="OSM main road")],
          loc="lower left", fontsize=7, framealpha=0.85, facecolor="#555555", labelcolor="white")
ax.set_title(f"Zoom on hotspot 1: {hs.loc[0, 'place']} ({hs.loc[0, 'borough']}), 8 × 6 km", fontsize=9)
fig.savefig(OUT / "04-hotspots-osm.png")
plt.close(fig)

fig, ax = plt.subplots(figsize=(8, 4.0), constrained_layout=True)
yy = np.arange(len(met))
ax.barh(yy + 0.2, met["air_11UTC_C"], height=0.38, color="#1c7ed6", label="air temperature, 11:00 UTC (station)")
ax.barh(yy - 0.2, met["LST_median_1km_C"], height=0.38, color="#e8590c", label="LST, median within ~1 km (Landsat 9, 10:52 UTC)")
for i, rrow in met.iterrows():
    ax.text(rrow["LST_median_1km_C"] + 0.4, i - 0.2, f"{rrow['LST_median_1km_C']:.1f}", va="center", fontsize=8)
    ax.text(rrow["air_11UTC_C"] + 0.4, i + 0.2, f"{rrow['air_11UTC_C']:.1f}", va="center", fontsize=8)
ax.set_yticks(yy, met["station"])
ax.invert_yaxis()
ax.set_xlabel("°C on 18 July 2022")
ax.set_xlim(0, max(met["LST_median_1km_C"].max(), met["air_11UTC_C"].max()) * 1.18)
ax.legend(loc="upper center", bbox_to_anchor=(0.5, -0.14), ncol=2, fontsize=8, frameon=False)
ax.spines[["top", "right"]].set_visible(False)
fig.savefig(OUT / "05-lst-vs-air-temperature.png")
plt.close(fig)

section("outputs")
for f in sorted(OUT.iterdir()):
    if f.name != "console.txt":
        print(f"{f.name:<32} {f.stat().st_size / 1024:8.0f} KB")