"""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")