make_maps.py

import argparse
import json
import math
import time
from pathlib import Path

import branca.colormap as bcm
import folium
import folium.plugins as fplugins
import geopandas as gpd
import matplotlib
matplotlib.use("Agg")
import matplotlib.colors as mcolors
import matplotlib.patches as mpatches
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import pandas as pd
import rasterio
from rasterio.features import geometry_mask
from rasterio.warp import calculate_default_transform, reproject, Resampling
from shapely.ops import unary_union

try:
    from sklearn.cluster import KMeans
    from sklearn.ensemble import RandomForestRegressor
    from sklearn.preprocessing import StandardScaler
    _SK = True
except ImportError:
    _SK = False

try:
    import cupy as cp
    cp.cuda.Device(0).use()
    _GPU = True
except Exception:
    _GPU = False

_RAW  = Path("data/charlotte/raw")
_PROC = Path("data/charlotte/processed")
_OUT  = Path("outputs/charlotte")

CITY  = "Charlotte, NC"
DPI   = 300
FW, FH = 11.5, 9.5


def _to_wgs84(src_path, ref_transform=None, ref_shape=None, categorical=False):
    res = Resampling.nearest if categorical else Resampling.bilinear
    with rasterio.open(str(src_path)) as src:
        nd = src.nodata if src.nodata is not None else -9999.0
        if ref_transform is None:
            dst_transform, width, height = calculate_default_transform(
                src.crs, "EPSG:4326", src.width, src.height, *src.bounds)
        else:
            dst_transform = ref_transform
            height, width = ref_shape
        dst = np.full((height, width), nd, dtype="float32")
        reproject(source=rasterio.band(src, 1), destination=dst,
                  src_transform=src.transform, src_crs=src.crs,
                  dst_transform=dst_transform, dst_crs="EPSG:4326",
                  resampling=res, src_nodata=nd, dst_nodata=nd)
    arr = np.where(np.abs(dst - nd) < 0.5, np.nan, dst).astype("float32")
    west = dst_transform.c; north = dst_transform.f
    east = west + dst_transform.a * width
    south = north + dst_transform.e * height
    return arr, dst_transform, (height, width), (west, east, south, north)


def load_rasters(raw, proc):
    lst_raw, ref_transform, ref_shape, extent = _to_wgs84(raw / "lst.tif")
    lst  = np.where((lst_raw > -999) & np.isfinite(lst_raw), lst_raw, np.nan)
    ndvi, _, _, _ = _to_wgs84(raw  / "ndvi.tif", ref_transform, ref_shape)
    uhi,  _, _, _ = _to_wgs84(proc / "uhi_intensity.tif", ref_transform, ref_shape)
    nlcd_f, _, _, _ = _to_wgs84(raw / "nlcd.tif", ref_transform, ref_shape, categorical=True)
    nlcd = np.where(np.isnan(nlcd_f), 0, nlcd_f).astype("int16")
    return lst, ndvi, uhi, nlcd, extent, ref_transform, ref_shape


def load_trend(out):
    p = out / "lst_trend.tif"
    if not p.exists():
        return None
    arr, _, _, _ = _to_wgs84(p)
    return arr


def load_gdf(proc):
    gdf = gpd.read_file(str(proc / "tracts_solutions.geojson")).to_crs("EPSG:4326")
    num_cols = ["mean_lst", "mean_ndvi", "mean_uhi", "heat_risk_score", "ndvi_deficit",
                "svi_score", "median_income", "pct_minority", "pct_poverty",
                "pct_elderly", "pct_under5", "pct_old_housing", "total_pop"]
    for c in num_cols:
        if c in gdf.columns:
            gdf[c] = pd.to_numeric(gdf[c], errors="coerce")
    if "equity_priority" in gdf.columns:
        gdf["equity_priority"] = gdf["equity_priority"].fillna(False).astype(bool)
    return gdf


def county_mask(gdf, transform, shape):
    county  = unary_union(gdf.to_crs("EPSG:4326").geometry)
    outside = geometry_mask([county.__geo_interface__],
                            transform=transform, invert=False, out_shape=shape)
    return ~outside


def _north_arrow(ax):
    ax.annotate("", xy=(0.955, 0.955), xytext=(0.955, 0.890),
                xycoords="axes fraction", textcoords="axes fraction",
                arrowprops=dict(arrowstyle="-|>", color="black", lw=2.2, mutation_scale=14))
    ax.text(0.955, 0.965, "N", transform=ax.transAxes,
            ha="center", va="bottom", fontsize=11, fontweight="bold")


def _scale_bar(ax, west, east, south, north):
    lat   = (south + north) / 2.0
    lon_w = 10.0 / (111.32 * math.cos(math.radians(lat)))
    sx, sy = east - west, north - south
    x0, y0 = west + sx * 0.05, south + sy * 0.04
    x1, tk = x0 + lon_w, sy * 0.011
    ax.plot([x0, x1], [y0, y0],       "k-", lw=3.5, solid_capstyle="butt", zorder=8)
    ax.plot([x0, x0], [y0-tk, y0+tk], "k-", lw=1.8, zorder=8)
    ax.plot([x1, x1], [y0-tk, y0+tk], "k-", lw=1.8, zorder=8)
    ax.text((x0+x1)/2, y0+sy*0.017, "10 km", ha="center", va="bottom",
            fontsize=8, zorder=8,
            bbox=dict(facecolor="white", alpha=0.8, edgecolor="none",
                      boxstyle="round,pad=0.1"))


def _source(fig, txt):
    fig.text(0.012, 0.006, f"Data: {txt}", fontsize=6.5, color="#666666", va="bottom")


def _style(ax, title, fontsize=11):
    ax.set_title(title, fontsize=fontsize, fontweight="bold", pad=8, wrap=True)
    ax.set_xlabel("Longitude", fontsize=8.5, labelpad=2)
    ax.set_ylabel("Latitude",  fontsize=8.5, labelpad=2)
    ax.tick_params(labelsize=7)
    for sp in ax.spines.values():
        sp.set_linewidth(0.4)


def _tracts(ax, gdf):
    ep = gdf.get("equity_priority", pd.Series(False, index=gdf.index)).astype(bool)
    gdf[~ep].boundary.plot(ax=ax, color="#888888", lw=0.18, alpha=0.4, zorder=4)
    if ep.any():
        gdf[ep].boundary.plot(ax=ax, color="#d73027", lw=0.85, alpha=0.88, zorder=5)


def _county_outline(ax, gdf):
    gpd.GeoDataFrame(geometry=[unary_union(gdf.geometry)], crs=gdf.crs
                     ).boundary.plot(ax=ax, color="#111111", lw=1.2, zorder=6)


def _save(fig, out, name):
    out.mkdir(parents=True, exist_ok=True)
    fig.savefig(str(out / name), dpi=DPI, bbox_inches="tight")
    plt.close(fig)


def _raster_fig(data, extent, cmap, vmin, vmax, title, cbar_label,
                source_txt, gdf, out, fname, mask=None):
    west, east, south, north = extent
    if mask is not None:
        data = np.where(mask, data, np.nan)
    cm = plt.get_cmap(cmap).copy()
    cm.set_bad(color="white", alpha=1.0)
    fig, ax = plt.subplots(figsize=(FW, FH))
    im = ax.imshow(data, cmap=cm, vmin=vmin, vmax=vmax,
                   extent=[west, east, south, north], aspect="auto",
                   interpolation="bilinear")
    _tracts(ax, gdf)
    ax.set_xlim(west, east); ax.set_ylim(south, north)
    cb = plt.colorbar(im, ax=ax, shrink=0.62, pad=0.015, aspect=22,
                      location="right")
    cb.set_label(cbar_label, fontsize=10)
    cb.ax.tick_params(labelsize=9)
    _style(ax, title, fontsize=11)
    _north_arrow(ax)
    _scale_bar(ax, west, east, south, north)
    _source(fig, source_txt)
    plt.tight_layout(rect=[0, 0.03, 1, 0.97])
    _save(fig, out, fname)


def _choropleth(gdf, col, cmap, cbar_label, title, source_txt, fname, out,
                vmin=None, vmax=None):
    gdf = gdf.copy()
    gdf[col] = pd.to_numeric(gdf[col], errors="coerce")
    b = gdf.total_bounds
    fig, ax = plt.subplots(figsize=(FW, FH))
    gdf.plot(column=col, cmap=cmap, ax=ax, vmin=vmin, vmax=vmax, legend=True,
             legend_kwds={"label": cbar_label, "orientation": "vertical",
                          "shrink": 0.58, "pad": 0.015, "aspect": 20},
             missing_kwds={"color": "#d9d9d9", "label": "No data"},
             linewidth=0.2, edgecolor="white")
    _county_outline(ax, gdf)
    ep = gdf.get("equity_priority", pd.Series(False, index=gdf.index)).astype(bool)
    if ep.any():
        gdf[ep].boundary.plot(ax=ax, color="#d73027", lw=0.85, alpha=0.88, zorder=5)
    ax.set_xlim(b[0], b[2]); ax.set_ylim(b[1], b[3])
    _style(ax, title, fontsize=11)
    _north_arrow(ax)
    _scale_bar(ax, b[0], b[2], b[1], b[3])
    _source(fig, source_txt)
    plt.tight_layout(rect=[0, 0.03, 1, 0.97])
    _save(fig, out, fname)


def map_01_3d_lst(lst, extent, gdf, out, mask=None):
    west, east, south, north = extent
    data = np.where(mask, lst, np.nan) if mask is not None else lst.copy()
    step = max(1, data.shape[0] // 280)
    ds   = data[::step, ::step]
    if _GPU:
        ds = cp.asnumpy(cp.asarray(ds))

    rows, cols = ds.shape
    X, Y = np.meshgrid(np.linspace(west, east, cols), np.linspace(north, south, rows))
    Z    = ds.copy()
    col_med  = np.nanmedian(Z, axis=0)
    fallback = float(np.nanmedian(Z)) if not np.all(np.isnan(Z)) else 35.0
    for i in range(Z.shape[0]):
        bad = np.isnan(Z[i])
        if bad.any():
            Z[i][bad] = np.where(np.isnan(col_med[bad]), fallback, col_med[bad])

    vmin = float(np.nanpercentile(Z, 1))
    vmax = float(np.nanpercentile(Z, 99))
    norm = mcolors.Normalize(vmin=vmin, vmax=vmax)
    cmap = plt.get_cmap("inferno")

    fig = plt.figure(figsize=(13, 9), facecolor="#0d0d0d")
    ax  = fig.add_subplot(111, projection="3d")
    ax.set_facecolor("#0d0d0d")
    ax.plot_surface(X, Y, Z, facecolors=cmap(norm(Z)),
                    rstride=1, cstride=1, linewidth=0, antialiased=True,
                    shade=False, alpha=0.97)

    for pane in [ax.xaxis.pane, ax.yaxis.pane, ax.zaxis.pane]:
        pane.fill = False
        pane.set_edgecolor("#333333")

    ax.set_xlabel("Longitude", fontsize=8, color="#cccccc", labelpad=8)
    ax.set_ylabel("Latitude",  fontsize=8, color="#cccccc", labelpad=8)
    ax.set_zlabel("LST (°C)",  fontsize=8, color="#cccccc", labelpad=8)
    ax.tick_params(colors="#aaaaaa", labelsize=6.5)
    ax.view_init(elev=30, azim=-115)

    fig.text(0.5, 0.96,
             f"{CITY}  -  Land Surface Temperature: 3D Heat Landscape",
             ha="center", va="top", fontsize=13, fontweight="bold", color="white")
    fig.text(0.5, 0.925,
             "Summer 2023  |  Landsat 8/9  |  Mecklenburg County  |  "
             "Red outlines = equity-priority tracts",
             ha="center", va="top", fontsize=9, color="#aaaaaa")

    sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm)
    sm.set_array([])
    cb = fig.colorbar(sm, ax=ax, shrink=0.45, pad=0.06, aspect=16,
                      location="right")
    cb.set_label("Land Surface Temperature (°C)", fontsize=8.5, color="white")
    cb.ax.tick_params(colors="white", labelsize=7.5)
    cb.outline.set_edgecolor("white")

    Z_clipped = np.clip(Z, vmin, vmax)
    peak = np.unravel_index(np.nanargmax(Z_clipped), Z_clipped.shape)
    ax.text(X[peak], Y[peak], Z_clipped[peak] + 1.5,
            f"Peak\n{Z_clipped[peak]:.1f}°C",
            fontsize=9, color="#ff9966", fontweight="bold", ha="center",
            bbox=dict(facecolor="#0d0d0d", alpha=0.6, edgecolor="none",
                      boxstyle="round,pad=0.2"))

    fig.text(0.012, 0.01,
             "Data: NASA Landsat 8/9 C2 L2 | Microsoft Planetary Computer | Summer 2023",
             fontsize=6.5, color="#888888", va="bottom")
    _save(fig, out, "map_01_3d_lst.png")


def map_02_canopy_income(gdf, out):
    gdf = gdf.copy()
    for c in ["mean_ndvi", "median_income", "mean_lst", "total_pop"]:
        gdf[c] = pd.to_numeric(gdf[c], errors="coerce")
    sub = gdf.dropna(subset=["mean_ndvi", "median_income", "mean_lst"])

    x   = sub["median_income"].values / 1000.0
    y   = sub["mean_ndvi"].values
    lst = sub["mean_lst"].values
    pop = pd.to_numeric(sub["total_pop"], errors="coerce").fillna(3000).values
    ep  = sub.get("equity_priority",
                  pd.Series(False, index=sub.index)).astype(bool).values

    coeffs = np.polyfit(x, y, 1)
    x_line = np.linspace(x.min(), x.max(), 200)
    y_line = np.polyval(coeffs, x_line)
    ss_res = np.sum((y - np.polyval(coeffs, x)) ** 2)
    r2     = 1 - ss_res / np.sum((y - y.mean()) ** 2)

    corr_coeffs = np.polyfit(y, lst, 1)
    lst_per_ndvi = corr_coeffs[0]
    ndvi_per_10k = coeffs[0] * 10
    lst_per_10k  = ndvi_per_10k * lst_per_ndvi

    norm = mcolors.Normalize(vmin=np.nanpercentile(lst, 2),
                             vmax=np.nanpercentile(lst, 98))
    cmap_lst = plt.get_cmap("plasma")
    sizes = 25 + 180 * (pop - pop.min()) / (pop.max() - pop.min() + 1)

    fig, ax = plt.subplots(figsize=(11, 7.5))

    sc = ax.scatter(x[~ep], y[~ep], c=lst[~ep], cmap=cmap_lst, norm=norm,
                    s=sizes[~ep], alpha=0.60, edgecolors="none", zorder=3,
                    label=f"Standard tracts ({(~ep).sum()})")
    ax.scatter(x[ep], y[ep], c=lst[ep], cmap=cmap_lst, norm=norm,
               s=sizes[ep] + 35, alpha=0.95, edgecolors="#d73027",
               linewidths=1.2, zorder=4, marker="D",
               label=f"Equity-priority ({ep.sum()})  - low income, low NDVI, high heat")

    ax.plot(x_line, y_line, color="#222222", lw=2.2, linestyle="--", zorder=5,
            label=f"OLS: R² = {r2:.3f}  (income strongly predicts canopy cover)")

    cb = plt.colorbar(sc, ax=ax, shrink=0.72, pad=0.01, aspect=22)
    cb.set_label("Mean LST (°C)", fontsize=9)
    cb.ax.tick_params(labelsize=8)

    ax.text(0.97, 0.52,
            f"For every $10K drop in median income:\n"
            f"  -> NDVI falls ~{abs(ndvi_per_10k):.3f}  (less tree canopy)\n"
            f"  -> Surface temperature rises ~{abs(lst_per_10k):.2f}°C\n\n"
            f"Poverty doesn't cause heat directly -\n"
            f"it causes disinvestment in tree canopy,\n"
            f"and bare land absorbs solar radiation.",
            transform=ax.transAxes, va="center", ha="right", fontsize=9.5,
            bbox=dict(facecolor="#fffbe6", edgecolor="#f0c040",
                      boxstyle="round,pad=0.55", alpha=0.97))

    ax.axhline(0.25, color="#d73027", lw=1.2, linestyle=":", alpha=0.75)
    ax.text(x.max() * 0.97, 0.257, "NDVI < 0.25\n(street tree trigger)",
            ha="right", fontsize=8, color="#d73027")

    ax.set_xlabel("Median Household Income ($K per tract)", fontsize=11)
    ax.set_ylabel("Mean NDVI  (tree canopy density proxy)", fontsize=11)
    ax.set_title(
        f"{CITY} - The Canopy-Poverty Link: Why Low-Income Tracts Overheat\n"
        f"Income predicts tree cover (R²={r2:.3f}), tree cover drives surface heat  "
        f"|  dot size = tract population",
        fontsize=11, fontweight="bold", pad=10)
    ax.legend(loc="lower right", fontsize=9, framealpha=0.92,
              edgecolor="#cccccc")
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.grid(alpha=0.18, linestyle="--")
    _source(fig, "Landsat 8/9 NDVI + LST | US Census ACS 5-Year 2022 B19013 | 305 Mecklenburg Co. tracts")
    plt.tight_layout(rect=[0, 0.03, 1, 1])
    _save(fig, out, "map_02_canopy_income.png")


def map_03_trend(slope, extent, gdf, out, mask=None):
    if slope is None:
        return
    finite = slope[np.isfinite(slope)]
    absmax = float(np.percentile(np.abs(finite), 97)) if len(finite) else 0.3
    _raster_fig(
        data=slope, extent=extent, cmap="RdBu_r",
        vmin=-absmax, vmax=absmax,
        title=f"{CITY} - Surface Temperature Trend  2019–2023\nred = warming faster  |  blue = cooling  |  red outlines = equity-priority tracts",
        cbar_label="LST trend (°C per year)",
        source_txt="Pixel-wise OLS regression on 5 summers of Landsat 8/9 LST | Microsoft Planetary Computer",
        gdf=gdf, out=out, fname="map_03_trend.png", mask=mask)


_BV = {
    (0,0):"#e8e8e8", (0,1):"#b0d5df", (0,2):"#64acbe",
    (1,0):"#e4d9ac", (1,1):"#ad9ea5", (1,2):"#627f8c",
    (2,0):"#c85a5a", (2,1):"#985356", (2,2):"#574249",
}

def map_04_bivariate(gdf, out):
    gdf = gdf.copy()
    for c in ["heat_risk_score", "svi_score"]:
        gdf[c] = pd.to_numeric(gdf[c], errors="coerce")
    gdf["h_bin"] = pd.qcut(gdf["heat_risk_score"].rank(method="first"), 3,
                            labels=[0,1,2]).astype(int)
    gdf["v_bin"] = pd.qcut(gdf["svi_score"].rank(method="first"), 3,
                            labels=[0,1,2]).astype(int)
    gdf["bv_color"] = [_BV.get((h,v), "#e8e8e8")
                       for h,v in zip(gdf["h_bin"], gdf["v_bin"])]

    b = gdf.total_bounds
    fig, ax = plt.subplots(figsize=(FW, FH))
    for color in _BV.values():
        sub = gdf[gdf["bv_color"] == color]
        if not sub.empty:
            sub.plot(ax=ax, color=color, linewidth=0.22, edgecolor="white")

    _county_outline(ax, gdf)
    ep = gdf.get("equity_priority", pd.Series(False, index=gdf.index)).astype(bool)
    if ep.any():
        gdf[ep].boundary.plot(ax=ax, color="#d73027", lw=0.85, alpha=0.88, zorder=5)
    ax.set_xlim(b[0], b[2]); ax.set_ylim(b[1], b[3])

    leg_ax = ax.inset_axes([0.015, 0.03, 0.17, 0.17])
    for (hr, sv), color in _BV.items():
        leg_ax.fill_between([sv, sv+1], [hr, hr], [hr+1, hr+1], color=color)
    leg_ax.set_xlim(0, 3); leg_ax.set_ylim(0, 3)
    leg_ax.set_xticks([0,1.5,3]); leg_ax.set_yticks([0,1.5,3])
    leg_ax.set_xticklabels(["Low","","High"], fontsize=6)
    leg_ax.set_yticklabels(["Low","","High"], fontsize=6)
    leg_ax.set_xlabel("Vulnerability ->", fontsize=6.5, labelpad=1)
    leg_ax.set_ylabel("Heat ->", fontsize=6.5, labelpad=1)
    leg_ax.tick_params(length=0)
    for sp in leg_ax.spines.values():
        sp.set_linewidth(0.5)
    leg_ax.text(2.5, 2.5, "*", ha="center", va="center",
                fontsize=10, color="white", fontweight="bold")

    n_worst = int(((gdf["h_bin"]==2) & (gdf["v_bin"]==2)).sum())
    ax.text(b[2] - (b[2]-b[0])*0.01, b[3] - (b[3]-b[1])*0.10,
            f"{n_worst} tracts:\nhigh heat + high vulnerability",
            ha="right", va="top", fontsize=8.5,
            bbox=dict(facecolor="#574249", alpha=0.92, edgecolor="none",
                      boxstyle="round,pad=0.40"), color="white", fontweight="bold")

    _style(ax, f"{CITY} - Bivariate: Heat Risk x Social Vulnerability\n"
               "Dark purple = highest heat + highest vulnerability  |  Red outlines = equity-priority tracts",
               fontsize=10)
    _north_arrow(ax)
    _scale_bar(ax, b[0], b[2], b[1], b[3])
    _source(fig, "Heat Risk Score + SVI. Landsat 8/9 + US Census ACS 2022.")
    plt.tight_layout(rect=[0, 0.03, 1, 0.97])
    _save(fig, out, "map_04_bivariate.png")


def map_05_health_risk(gdf, out):
    gdf = gdf.copy()
    for c in ["mean_lst","total_pop","pct_elderly","pct_under5"]:
        gdf[c] = pd.to_numeric(gdf[c], errors="coerce").fillna(0)

    safe_threshold = 35.0
    gdf["heat_excess"] = np.maximum(0, gdf["mean_lst"] - safe_threshold)
    gdf["vuln_weight"] = 1.0 + (gdf["pct_elderly"]/100)*2.0 + (gdf["pct_under5"]/100)*1.5
    gdf["hhe"] = gdf["total_pop"] * gdf["heat_excess"] * gdf["vuln_weight"]
    mx = gdf["hhe"].max()
    gdf["hhe_norm"] = gdf["hhe"] / mx if mx > 0 else 0

    total_exposed = int(gdf.loc[gdf["heat_excess"] > 0, "total_pop"].sum())

    _choropleth(
        gdf, "hhe_norm", "OrRd",
        "Heat Health Exposure Index  (0 = low  |  1 = maximum)",
        f"{CITY} - Population Heat Health Exposure Index\n"
        f"pop x (LST-35°C) x vulnerability weight  |  "
        f"{total_exposed:,} residents in tracts above 35°C threshold",
        "Landsat 8/9 LST + ACS 2022 elderly/youth demographics + population. "
        "35°C threshold: CDC heat illness guideline. Red outlines = equity-priority tracts.",
        "map_05_health_risk.png", out)


def map_06_ndvi(ndvi, extent, gdf, out, mask=None):
    _raster_fig(
        data=ndvi, extent=extent, cmap="RdYlGn",
        vmin=0.05, vmax=0.75,
        title=f"{CITY} - Vegetation Cover: NDVI  (Summer 2023)\nred = bare / impervious  |  green = dense forest canopy",
        cbar_label="NDVI",
        source_txt="NASA Landsat 8/9 NIR Band 5 and Red Band 4 | Microsoft Planetary Computer | Summer 2023",
        gdf=gdf, out=out, fname="map_06_ndvi.png", mask=mask)


def map_07_landcover(nlcd, extent, gdf, out, mask=None):
    west, east, south, north = extent
    NLCD_PALETTE = {
        11: ("#486DA2","Open Water"),
        21: ("#E8D1D1","Developed, Open Space"),
        22: ("#E29E8C","Developed, Low Intensity"),
        23: ("#FF0000","Developed, Med. Intensity"),
        24: ("#B50000","Developed, High Intensity"),
        41: ("#85C77E","Deciduous Forest"),
        42: ("#38814E","Evergreen Forest"),
        43: ("#D4E7B0","Mixed Forest"),
        52: ("#DCCA8F","Shrub / Scrub"),
        71: ("#FDE9AA","Grassland"),
        81: ("#FBF65D","Pasture / Hay"),
        82: ("#CA9146","Cultivated Crops"),
        90: ("#C8E6F8","Woody Wetlands"),
        95: ("#64B3D5","Emergent Wetlands"),
    }
    rgb = np.full((*nlcd.shape, 3), 0.92, dtype="float32")
    for code, (hx, _) in NLCD_PALETTE.items():
        r, g, b = [int(hx[i:i+2], 16)/255.0 for i in (1,3,5)]
        rgb[nlcd == code] = [r, g, b]
    if mask is not None:
        rgb[~mask] = [1.0, 1.0, 1.0]

    fig, ax = plt.subplots(figsize=(FW, FH))
    ax.imshow(rgb, extent=[west, east, south, north], aspect="auto")
    _tracts(ax, gdf)
    ax.set_xlim(west, east); ax.set_ylim(south, north)

    patches = [mpatches.Patch(color=hx, label=lbl)
               for code, (hx, lbl) in NLCD_PALETTE.items() if np.any(nlcd == code)]
    ax.legend(handles=patches, loc="center left",
              bbox_to_anchor=(1.01, 0.5),
              fontsize=8, title="Land Cover (NLCD 2021)", title_fontsize=9,
              framealpha=0.95, handlelength=1.4, edgecolor="#cccccc",
              ncol=1, borderpad=0.6)

    _style(ax, f"{CITY} - Land Cover Classification  (NLCD 2021)\n"
               "Red/brown = impervious urban surface (heat source)  |  Green = vegetation (natural cooling)",
               fontsize=10)
    _north_arrow(ax)
    _scale_bar(ax, west, east, south, north)
    _source(fig, "USGS MRLC National Land Cover Database 2021")
    plt.tight_layout(rect=[0, 0.03, 0.80, 0.97])
    _save(fig, out, "map_07_landcover.png")


def map_08_ml_importance(gdf, out):
    if not _SK:
        return

    gdf = gdf.copy()
    feature_cols = ["mean_ndvi","pct_minority","pct_poverty","median_income",
                    "svi_score","pct_elderly","pct_under5","pct_old_housing"]
    labels = ["NDVI (Vegetation Cover)","% Minority Population",
              "% Below Poverty Line","Median Household Income ($)",
              "Social Vulnerability Index","% Elderly (65+)",
              "% Children Under 5","% Pre-1940 Housing"]

    feat_df = gdf[feature_cols].copy()
    for c in feature_cols:
        feat_df[c] = pd.to_numeric(feat_df[c], errors="coerce")
    feat_df = feat_df.fillna(feat_df.median())
    target  = pd.to_numeric(gdf["mean_lst"], errors="coerce")
    target  = target.fillna(target.median())

    rf = RandomForestRegressor(n_estimators=400, random_state=42, n_jobs=-1)
    rf.fit(feat_df.values, target.values)
    importances = rf.feature_importances_
    r2 = rf.score(feat_df.values, target.values)

    order  = np.argsort(importances)
    s_lab  = [labels[i] for i in order]
    s_imp  = importances[order]
    colors = ["#d73027" if v > 0.15 else "#fc8d59" if v > 0.08 else "#fee090"
              for v in s_imp]

    fig, ax = plt.subplots(figsize=(10, 6.5))
    bars = ax.barh(range(len(s_lab)), s_imp * 100, color=colors,
                   edgecolor="white", linewidth=0.5, height=0.62)
    for bar, v in zip(bars, s_imp):
        ax.text(bar.get_width() + 0.4, bar.get_y() + bar.get_height()/2,
                f"{v*100:.1f}%", va="center", ha="left", fontsize=9.5,
                fontweight="bold", color="#222222")

    ax.set_yticks(range(len(s_lab)))
    ax.set_yticklabels(s_lab, fontsize=10)
    ax.set_xlabel("Feature Importance (%)", fontsize=10.5)
    ax.set_xlim(0, s_imp.max() * 100 * 1.28)
    ax.set_title(
        f"{CITY} - What Drives Urban Heat?  (Random Forest, 400 trees)\n"
        f"Predicting Land Surface Temperature from environmental + social factors  |  R² = {r2:.3f}",
        fontsize=11, fontweight="bold", pad=10)
    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.axvline(0, color="#333333", lw=0.8)

    top_feat = s_lab[-1]
    ax.text(s_imp.max()*100*0.58, 0.6,
            f"Key finding:\n{top_feat}\nis the #1 driver of\nheat (not race or poverty)",
            fontsize=9, style="italic",
            bbox=dict(facecolor="#fff3cd", edgecolor="#ffc107",
                      boxstyle="round,pad=0.45", alpha=0.95))

    patches = [mpatches.Patch(color="#d73027", label=">15% importance"),
               mpatches.Patch(color="#fc8d59", label="8–15%"),
               mpatches.Patch(color="#fee090", label="<8%")]
    ax.legend(handles=patches, loc="lower right", fontsize=8.5,
              framealpha=0.9, edgecolor="#cccccc")
    _source(fig, "RandomForestRegressor (scikit-learn) | 400 trees | Landsat 8/9 LST + Census ACS 2022")
    plt.tight_layout(rect=[0, 0.04, 1, 1])
    _save(fig, out, "map_08_ml_importance.png")


_CLUSTER_META = [
    ("Cool Wealthy Suburb",     "#2166ac", "Low heat | High income | Dense canopy"),
    ("Moderate Mixed",          "#74add1", "Average heat | Mixed income"),
    ("Green Transitional",      "#4dac26", "Moderate heat | Rising NDVI"),
    ("Hot Vulnerable Core",     "#fc8d59", "High heat | High poverty | Low canopy"),
    ("Extreme Heat Zone",       "#d73027", "Highest heat | Critical vulnerability"),
]

def map_09_ml_clusters(gdf, out):
    if not _SK:
        return

    gdf = gdf.copy()
    features = ["mean_lst","mean_ndvi","heat_risk_score","svi_score",
                "median_income","pct_minority","pct_poverty"]
    feat_df  = gdf[features].copy()
    for c in features:
        feat_df[c] = pd.to_numeric(feat_df[c], errors="coerce")
    feat_df = feat_df.fillna(feat_df.median())

    X = StandardScaler().fit_transform(feat_df.values)
    labels = KMeans(n_clusters=5, random_state=42, n_init=30).fit_predict(X)
    mean_lst = [feat_df["mean_lst"].values[labels==k].mean() for k in range(5)]
    remap  = {old: new for new, old in enumerate(np.argsort(mean_lst))}
    gdf["cluster"] = np.array([remap[l] for l in labels])

    b = gdf.total_bounds
    fig, ax = plt.subplots(figsize=(FW + 1.5, FH))
    for k, (name, color, desc) in enumerate(_CLUSTER_META):
        sub = gdf[gdf["cluster"] == k]
        if not sub.empty:
            sub.plot(ax=ax, color=color, linewidth=0.22, edgecolor="white")

    _county_outline(ax, gdf)
    ep = gdf.get("equity_priority", pd.Series(False, index=gdf.index)).astype(bool)
    if ep.any():
        gdf[ep].boundary.plot(ax=ax, color="#333333", lw=0.6, alpha=0.6, zorder=6)

    ax.set_xlim(b[0], b[2]); ax.set_ylim(b[1], b[3])

    legend_handles = []
    for k, (name, color, desc) in enumerate(_CLUSTER_META):
        sub = gdf[gdf["cluster"] == k]
        n = len(sub)
        legend_handles.append(
            mpatches.Patch(color=color,
                           label=f"C{k+1}: {name} ({n} tracts)\n         {desc}"))

    leg = ax.legend(handles=legend_handles, loc="center left", bbox_to_anchor=(1.01, 0.5),
                    fontsize=8.5, title="Neighborhood Archetypes  (K-Means, k=5)",
                    title_fontsize=9.5, framealpha=0.96, edgecolor="#cccccc",
                    handlelength=1.4, borderpad=0.8)

    _style(ax, f"{CITY} - ML Neighborhood Archetypes\n"
               "K-Means clustering on: LST | NDVI | Heat Risk | SVI | Income | Demographics",
               fontsize=10)
    _north_arrow(ax)
    _scale_bar(ax, b[0], b[2], b[1], b[3])
    _source(fig, "K-Means (k=5, 30 inits) on standardized features. Landsat 8/9 + ACS 2022.")
    plt.tight_layout(rect=[0, 0.03, 0.78, 0.97])
    _save(fig, out, "map_09_ml_clusters.png")


_COOLING_DELTA = {
    "rec_street_trees":    2.5,
    "rec_pocket_parks":    3.0,
    "rec_cool_roofs":      1.8,
    "rec_cool_pavement":   1.2,
    "rec_weatherization":  2.2,
    "rec_misting_stations":0.5,
    "rec_green_corridors": 2.0,
    "rec_shade_structures":0.8,
}

def map_10_cooling(gdf, out):
    gdf = gdf.copy()
    cool = np.zeros(len(gdf))
    for col, delta in _COOLING_DELTA.items():
        if col in gdf.columns:
            cool += gdf[col].fillna(False).astype(float) * delta
    gdf["cooling_potential"] = np.where(cool > 0, cool, np.nan)

    ep_mask = gdf.get("equity_priority", pd.Series(False, index=gdf.index)).astype(bool)
    avg_ep  = float(gdf.loc[ep_mask, "cooling_potential"].mean()) if ep_mask.any() else 0
    n_ep    = int(ep_mask.sum())
    n_cool  = int((cool > 0).sum())

    b = gdf.total_bounds
    vmax_cool = float(np.nanpercentile(gdf["cooling_potential"].dropna(), 98))

    fig, ax = plt.subplots(figsize=(FW, FH))
    gdf.plot(column="cooling_potential", cmap="YlGnBu", ax=ax,
             vmin=0.01, vmax=vmax_cool, legend=True,
             legend_kwds={"label": "Projected LST reduction (°C)",
                          "orientation": "vertical", "shrink": 0.58, "pad": 0.015, "aspect": 22},
             missing_kwds={"color": "#e8e8e8", "label": "No interventions needed"},
             linewidth=0.2, edgecolor="white")
    _county_outline(ax, gdf)
    if ep_mask.any():
        gdf[ep_mask].boundary.plot(ax=ax, color="#d73027", lw=0.9, alpha=0.92, zorder=5)
    ax.set_xlim(b[0], b[2]); ax.set_ylim(b[1], b[3])

    ax.text(0.03, 0.97,
            f"{n_cool} tracts receive green interventions\n"
            f"{n_ep} equity-priority tracts (red outline)\n"
            f"Avg projected cooling: -{avg_ep:.1f}°C\n"
            f"Best case: street trees + pocket parks = -5.5°C",
            transform=ax.transAxes, va="top", fontsize=9,
            bbox=dict(facecolor="white", edgecolor="#2b8cbe",
                      boxstyle="round,pad=0.5", alpha=0.93))

    _style(ax,
           f"{CITY} - Projected Cooling Potential per Tract\n"
           f"Gray = no interventions recommended  |  Red outlines = {n_ep} equity-priority tracts",
           fontsize=10)
    _north_arrow(ax)
    _scale_bar(ax, b[0], b[2], b[1], b[3])
    _source(fig, "Literature cooling estimates: street trees -2.5°C | pocket parks -3.0°C | "
                 "cool roofs -1.8°C | weatherization -2.2°C | green corridors -2.0°C.")
    plt.tight_layout(rect=[0, 0.03, 1, 0.97])
    _save(fig, out, "map_10_cooling.png")
    return cool


_INT_DATA = [
    ("rec_street_trees",    "Street Tree\nCanopy",      8,   2.5, 1.3, "#4dac26"),
    ("rec_pocket_parks",    "Pocket Parks /\nGreen Space",250,3.0, 1.5, "#1a9641"),
    ("rec_cool_roofs",      "Cool / Green\nRoof",      150, 1.8, 1.1, "#fdae61"),
    ("rec_cool_pavement",   "Reflective\nPavement",    400, 1.2, 1.0, "#d7191c"),
    ("rec_weatherization",  "Housing\nWeatherization", 300, 2.2, 2.0, "#2b83ba"),
    ("rec_misting_stations","Misting\nStations",        40, 0.5, 1.2, "#018571"),
    ("rec_green_corridors", "Green\nCorridor",         100, 2.0, 1.4, "#74c476"),
    ("rec_shade_structures","Bus Stop /\nShade",       100, 0.8, 1.3, "#7b3294"),
]

def map_11_roi(gdf, out):
    gdf = gdf.copy()
    records = []
    for col, label, cost_k, cooling, eq_mult, color in _INT_DATA:
        if col not in gdf.columns:
            continue
        mask = gdf[col].fillna(False).astype(bool)
        n    = int(mask.sum())
        if n == 0:
            continue
        pop = float(gdf.loc[mask, "total_pop"].sum()) \
              if "total_pop" in gdf.columns else n * 3000
        impact = cooling * eq_mult * math.sqrt(pop)
        records.append({"label":label,"cost_k":cost_k,"cooling":cooling,
                         "eq_mult":eq_mult,"color":color,"n":n,
                         "impact":impact,"pop":pop})
    if not records:
        return

    _LABEL_OFFSETS = {
        "Street Tree\nCanopy":        (  0,   26, "center"),
        "Pocket Parks /\nGreen Space":(-70,  -38, "right"),
        "Cool / Green\nRoof":         (  0,  -34, "center"),
        "Reflective\nPavement":       (  0,  -34, "center"),
        "Housing\nWeatherization":    ( 70,  -38, "left"),
        "Misting\nStations":          (  0,   32, "center"),
        "Green\nCorridor":            (  0,   26, "center"),
        "Bus Stop /\nShade":          (  0,  -34, "center"),
    }

    fig, ax = plt.subplots(figsize=(12, 7.5))
    max_n = max(r["n"] for r in records)
    for r in records:
        sz = 220 + 1600 * (r["n"] / max_n)
        yval = r["cooling"] * r["eq_mult"]
        ax.scatter(r["cost_k"], yval, s=sz, color=r["color"],
                   alpha=0.84, edgecolors="white", linewidths=1.5, zorder=3)
        dx, dy, ha = _LABEL_OFFSETS.get(r["label"], (0, 24, "center"))
        ax.annotate(r["label"],
                    xy=(r["cost_k"], yval),
                    xytext=(dx, dy), textcoords="offset points",
                    ha=ha, va="bottom" if dy >= 0 else "top",
                    fontsize=9, fontweight="bold",
                    arrowprops=dict(arrowstyle="-", color="#999999", lw=0.9,
                                   shrinkA=4, shrinkB=4),
                    bbox=dict(facecolor="white", alpha=0.92, edgecolor="#cccccc",
                              boxstyle="round,pad=0.25"))

    ax.set_xscale("log")
    ax.set_xlabel("Estimated cost per tract  ($K, log scale)", fontsize=11)
    ax.set_ylabel("Cooling impact score\n(°C reduction x equity multiplier)", fontsize=11)
    ax.set_title(
        f"{CITY} - Green Infrastructure ROI: Cost vs. Cooling Benefit per Tract\n"
        "Bubble size = # tracts where intervention is recommended\n"
        "Equity multiplier: weatherization 2.0x, pocket parks 1.5x, green corridors 1.4x",
        fontsize=10, fontweight="bold", pad=10)

    xlim = ax.get_xlim(); ylim = ax.get_ylim()
    costs  = [r["cost_k"] for r in records]
    impacts = [r["cooling"]*r["eq_mult"] for r in records]
    med_c = float(np.median(costs)); med_i = float(np.median(impacts))
    ax.axvline(med_c, color="#bbbbbb", lw=1.0, linestyle="--", alpha=0.7, zorder=1)
    ax.axhline(med_i, color="#bbbbbb", lw=1.0, linestyle="--", alpha=0.7, zorder=1)

    ax.text(xlim[0]*1.15, ylim[1]*0.97,
            "* Best ROI\nhigh impact, low cost",
            va="top", fontsize=9, color="#1a9641", fontweight="bold",
            bbox=dict(facecolor="#f0fff4", edgecolor="#1a9641",
                      boxstyle="round,pad=0.3", alpha=0.85))
    ax.text(xlim[1]*0.75, ylim[0] + (ylim[1]-ylim[0])*0.04,
            "x Avoid\nlow impact, high cost",
            va="bottom", ha="right", fontsize=9, color="#d73027",
            bbox=dict(facecolor="#fff0f0", edgecolor="#d73027",
                      boxstyle="round,pad=0.3", alpha=0.85))

    for nv, lbl in [(5,"5 tracts"),(40,"40 tracts"),(max_n,f"{max_n} tracts")]:
        ax.scatter([],[], s=220+1600*(nv/max_n), color="#aaaaaa",
                   alpha=0.65, edgecolors="white", label=lbl)
    ax.legend(title="# Tracts Benefiting", loc="lower left",
              fontsize=8.5, title_fontsize=9, framealpha=0.92)

    ax.spines["top"].set_visible(False)
    ax.spines["right"].set_visible(False)
    ax.grid(alpha=0.2, linestyle="--")
    _source(fig, "Cost estimates: peer-reviewed green infrastructure literature. "
                 "Cooling deltas: Bowler et al. 2010, Santamouris 2014.")
    plt.tight_layout(rect=[0, 0.04, 1, 1])
    _save(fig, out, "map_11_roi.png")


def _safe_f(val, default=0.0):
    try:
        v = float(val); return default if (v!=v) else v
    except (TypeError, ValueError):
        return default


def make_html_map(gdf, out):
    gdf = gdf.copy().to_crs("EPSG:4326")
    num_cols = ["mean_lst","mean_ndvi","heat_risk_score","svi_score",
                "median_income","pct_minority","pct_poverty","pct_elderly","total_pop"]
    for c in num_cols:
        if c not in gdf.columns: gdf[c] = 0.0
        gdf[c] = pd.to_numeric(gdf[c], errors="coerce").fillna(0.0)
    for c in ["equity_priority"]:
        if c not in gdf.columns: gdf[c] = False
        gdf[c] = gdf[c].fillna(False).astype(bool)
    for c in ["priority_tier","NAME"]:
        if c not in gdf.columns: gdf[c] = ""
        gdf[c] = gdf[c].fillna("").astype(str)
    rec_cols = [c for c in gdf.columns if c.startswith("rec_")]
    for c in rec_cols:
        gdf[c] = gdf[c].fillna(False).astype(bool)

    cool = np.zeros(len(gdf))
    for col, delta in _COOLING_DELTA.items():
        if col in gdf.columns:
            cool += gdf[col].astype(float) * delta
    gdf["cooling_potential"] = cool

    gdf["heat_excess"] = np.maximum(0, gdf["mean_lst"] - 35.0)
    gdf["hhe"] = gdf["total_pop"] * gdf["heat_excess"] * \
                 (1.0 + (gdf["pct_elderly"]/100)*2 + (gdf["pct_elderly"]/100)*1.5)
    mx = gdf["hhe"].max()
    gdf["hhe_norm"] = gdf["hhe"] / mx if mx > 0 else 0.0

    center = [gdf.geometry.centroid.y.mean(), gdf.geometry.centroid.x.mean()]
    m = folium.Map(location=center, zoom_start=11, tiles="CartoDB positron",
                   prefer_canvas=True)
    folium.TileLayer("CartoDB dark_matter", name="Dark Basemap", show=False).add_to(m)
    folium.TileLayer("Esri.WorldImagery",   name="Satellite",    show=False).add_to(m)

    m.get_root().html.add_child(folium.Element("""
    <div style="position:fixed;top:12px;left:55px;z-index:9999;background:white;
                padding:10px 16px;border-radius:7px;border:1px solid
                box-shadow:2px 3px 8px rgba(0,0,0,0.22);max-width:320px">
      <b style="font-size:14px;font-family:Arial,sans-serif;color:#c0392b">
        Charlotte - Urban Heat Mitigation
      </b><br>
      <span style="font-size:10.5px;color:#555;font-family:Arial,sans-serif">
        TSA Geospatial Technology 2025–2026<br>
        Red outlines = 98 equity-priority tracts for green investment.<br>
        <b>Toggle layers   |  Click any tract for details.</b>
      </span>
    </div>"""))

    def _popup(row):
        ep = bool(row.get("equity_priority", False))
        badge = ('<span style="background:#c0392b;color:white;padding:1px 6px;'
                 'border-radius:3px;font-size:10px;margin-left:4px">EQUITY PRIORITY</span>'
                 if ep else "")
        recs = [c.replace("rec_","").replace("_"," ").title()
                for c in rec_cols if row.get(c, False)]
        rec_txt = ", ".join(recs) if recs else "None"
        cool_v = _safe_f(row.get("cooling_potential", 0))
        lst_v  = _safe_f(row.get("mean_lst", 0))
        inc_v  = _safe_f(row.get("median_income", 0))
        ndvi_v = _safe_f(row.get("mean_ndvi", 0))
        return f"""
<div style="font-family:Arial,sans-serif;font-size:12px;width:265px;line-height:1.6">
  <b style="font-size:13px">{str(row.get("NAME",""))[:42]}</b>{badge}
  <hr style="margin:5px 0;border-color:#eee">
  <div style="display:grid;grid-template-columns:1fr 1fr;gap:2px 10px;font-size:11.5px">
    <div><b>Temp:</b> {lst_v:.1f}°C</div>
    <div><b>NDVI:</b> {ndvi_v:.3f}</div>
    <div><b>Income:</b> ${inc_v:,.0f}</div>
    <div><b>Risk:</b> {_safe_f(row.get("heat_risk_score")):.2f}</div>
    <div><b>Minority:</b> {_safe_f(row.get("pct_minority")):.0f}%</div>
    <div><b>Tier:</b> {str(row.get("priority_tier",""))[:12]}</div>
  </div>
  <hr style="margin:5px 0;border-color:#eee">
  <div style="background:#f0fff4;padding:5px 7px;border-radius:4px;font-size:11.5px">
    <b style="color:#1a7a3c">Projected cooling: -{cool_v:.1f}°C</b> if all interventions applied
  </div>
  <div style="margin-top:5px;font-size:11px"><b>Actions:</b> {rec_txt}</div>
</div>"""

    gj = json.loads(gdf.to_json())
    for feat, (_, row) in zip(gj["features"], gdf.iterrows()):
        feat["properties"]["_risk"]  = float(row["heat_risk_score"])
        feat["properties"]["_cool"]  = float(row["cooling_potential"])
        feat["properties"]["_hhe"]   = float(row["hhe_norm"])
        feat["properties"]["_ep"]    = bool(row["equity_priority"])
        feat["properties"]["_pop"]   = _popup(row)

    rq2  = float(gdf["heat_risk_score"].quantile(0.02))
    rq98 = float(gdf["heat_risk_score"].quantile(0.98))
    rcm  = bcm.LinearColormap(["#ffffcc","#fd8d3c","#800026"],
                               vmin=rq2, vmax=rq98,
                               caption="Composite Heat Risk Score")
    rcm.add_to(m)
    folium.GeoJson(gj, name="Heat Risk Score  *  (default)",
        style_function=lambda f: {
            "fillColor":   rcm(max(rq2, min(rq98, f["properties"]["_risk"]))),
            "fillOpacity": 0.78,
            "color":       "#d73027" if f["properties"]["_ep"] else "#888888",
            "weight":      2.0       if f["properties"]["_ep"] else 0.25,
        },
        highlight_function=lambda _: {"fillOpacity":0.93,"weight":1.8},
        tooltip=folium.GeoJsonTooltip(["NAME","heat_risk_score","priority_tier"],
                                      ["Tract","Heat Risk","Tier"], localize=True),
        popup=folium.GeoJsonPopup(fields=["_pop"], labels=False),
        show=True).add_to(m)

    cq2  = float(gdf["cooling_potential"].quantile(0.02))
    cq98 = float(gdf["cooling_potential"].quantile(0.98))
    if cq2 == cq98: cq98 = cq2 + 1.0
    ccm = bcm.LinearColormap(["#f7fcf5","#74c476","#005a32"],
                              vmin=cq2, vmax=cq98,
                              caption="Projected Cooling (°C)")
    ccm.add_to(m)
    cj = json.loads(gdf.to_json())
    for feat, (_, row) in zip(cj["features"], gdf.iterrows()):
        feat["properties"]["_cv"] = float(row["cooling_potential"])
        feat["properties"]["_ep"] = bool(row["equity_priority"])
        feat["properties"]["_pop"] = _popup(row)
    folium.GeoJson(cj, name="Projected Cooling Potential (°C)",
        style_function=lambda f: {
            "fillColor":   ccm(max(cq2, min(cq98, f["properties"]["_cv"]))),
            "fillOpacity": 0.80,
            "color":       "#d73027" if f["properties"]["_ep"] else "#666666",
            "weight":      2.0       if f["properties"]["_ep"] else 0.25,
        },
        highlight_function=lambda _: {"fillOpacity":0.92,"weight":1.8},
        tooltip=folium.GeoJsonTooltip(["NAME","cooling_potential"],
                                      ["Tract","Cooling Potential (°C)"], localize=True),
        popup=folium.GeoJsonPopup(fields=["_pop"], labels=False),
        show=False).add_to(m)

    lq2  = float(gdf["mean_lst"].quantile(0.02))
    lq98 = float(gdf["mean_lst"].quantile(0.98))
    lcm  = bcm.LinearColormap(["#2b83ba","#ffffbf","#d7191c"],
                               vmin=lq2, vmax=lq98,
                               caption="Mean Land Surface Temp (°C)")
    lcm.add_to(m)
    lj = json.loads(gdf.to_json())
    for feat, (_, row) in zip(lj["features"], gdf.iterrows()):
        feat["properties"]["_lv"] = float(row["mean_lst"])
        feat["properties"]["_ep"] = bool(row["equity_priority"])
        feat["properties"]["_pop"] = _popup(row)
    folium.GeoJson(lj, name="Land Surface Temperature (°C)",
        style_function=lambda f: {
            "fillColor":   lcm(max(lq2, min(lq98, f["properties"]["_lv"]))),
            "fillOpacity": 0.75,
            "color":       "#d73027" if f["properties"]["_ep"] else "#666666",
            "weight":      2.0       if f["properties"]["_ep"] else 0.25,
        },
        highlight_function=lambda _: {"fillOpacity":0.92,"weight":1.8},
        tooltip=folium.GeoJsonTooltip(["NAME","mean_lst"],["Tract","LST (°C)"], localize=True),
        popup=folium.GeoJsonPopup(fields=["_pop"], labels=False),
        show=False).add_to(m)

    hj = json.loads(gdf.to_json())
    hcm = bcm.LinearColormap(["#fff5f0","#fb6a4a","#67000d"],
                              vmin=0, vmax=1,
                              caption="Heat Health Exposure (0–1)")
    hcm.add_to(m)
    for feat, (_, row) in zip(hj["features"], gdf.iterrows()):
        feat["properties"]["_hv"] = float(row["hhe_norm"])
        feat["properties"]["_ep"] = bool(row["equity_priority"])
        feat["properties"]["_pop"] = _popup(row)
    folium.GeoJson(hj, name="Heat Health Exposure Index",
        style_function=lambda f: {
            "fillColor":   hcm(min(1.0, max(0.0, f["properties"]["_hv"]))),
            "fillOpacity": 0.78,
            "color":       "#d73027" if f["properties"]["_ep"] else "#666666",
            "weight":      2.0       if f["properties"]["_ep"] else 0.25,
        },
        highlight_function=lambda _: {"fillOpacity":0.92,"weight":1.8},
        tooltip=folium.GeoJsonTooltip(["NAME","hhe_norm"],
                                      ["Tract","Health Exposure (0–1)"], localize=True),
        popup=folium.GeoJsonPopup(fields=["_pop"], labels=False),
        show=False).add_to(m)

    tier_pal = {"Tier 1 - Critical":"#d73027","Tier 2 - High":"#fc8d59",
                "Tier 3 - Medium":"#fee08b","Tier 4 - Low":"#91bfdb"}
    tfg   = folium.FeatureGroup(name="Priority Tiers (1–4)", show=False)
    tjson = json.loads(gdf.to_json())
    for feat, (_, row) in zip(tjson["features"], gdf.iterrows()):
        feat["properties"]["_t"]   = str(row.get("priority_tier",""))
        feat["properties"]["_pop"] = _popup(row)
    folium.GeoJson(tjson,
        style_function=lambda f: {
            "fillColor":   tier_pal.get(f["properties"]["_t"], "#d9d9d9"),
            "fillOpacity": 0.75,
            "color":       "#ffffff", "weight": 0.3,
        },
        tooltip=folium.GeoJsonTooltip(["NAME","priority_tier"],["Tract","Tier"]),
        popup=folium.GeoJsonPopup(fields=["_pop"], labels=False),
    ).add_to(tfg)
    tfg.add_to(m)

    fplugins.MiniMap(toggle_display=True, tile_layer="CartoDB positron").add_to(m)
    fplugins.Fullscreen(position="topright").add_to(m)
    fplugins.MousePosition(position="bottomright").add_to(m)
    folium.LayerControl(collapsed=False, position="topright").add_to(m)

    out.mkdir(parents=True, exist_ok=True)
    m.save(str(out / "uhi_map.html"))


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--raw",  default=str(_RAW))
    parser.add_argument("--proc", default=str(_PROC))
    parser.add_argument("--out",  default=str(_OUT))
    args = parser.parse_args()
    raw, proc, out = Path(args.raw), Path(args.proc), Path(args.out)

    t0 = time.perf_counter()

    lst, ndvi, uhi, nlcd, extent, ref_transform, ref_shape = load_rasters(raw, proc)
    slope = load_trend(out)

    gdf  = load_gdf(proc)
    mask = county_mask(gdf, ref_transform, ref_shape)
    ep_n = int(gdf.get("equity_priority", pd.Series(False)).sum())


    map_01_3d_lst(lst, extent, gdf, out, mask=mask)
    map_02_canopy_income(gdf, out)
    map_03_trend(slope, extent, gdf, out, mask=mask)

    map_04_bivariate(gdf, out)
    map_05_health_risk(gdf, out)

    map_06_ndvi(ndvi, extent, gdf, out, mask=mask)
    map_07_landcover(nlcd, extent, gdf, out, mask=mask)

    map_08_ml_importance(gdf, out)
    map_09_ml_clusters(gdf, out)

    map_10_cooling(gdf, out)
    map_11_roi(gdf, out)

    make_html_map(gdf, out)

    elapsed = time.perf_counter() - t0


main()