Source code: Land cover accuracy assessment in Python: confusion matrix, Olofsson area estimates and sample size

"""Land cover map accuracy assessment in Python: ESA WorldCover 2021 against LUCAS 2022 field
survey points in Ireland, with Olofsson et al. (2014) stratified estimators, error-adjusted
area estimates and 95% confidence intervals, and sample-size planning for a new design.

TechnicalGIS tutorial: https://www.technicalgis.com/tutorials/land-cover-accuracy-assessment-python/
Data: ESA WorldCover 2021 v200 (CC BY 4.0), Eurostat LUCAS 2022 microdata and LUCAS master grid
(Eurostat, reuse authorised with acknowledgement), GISCO NUTS 2021 boundaries (EuroGeographics).
"""
import sys
import warnings
from importlib.metadata import version
from pathlib import Path

import numpy as np
import pandas as pd
import geopandas as gpd
import requests
import rasterio
from rasterio.features import geometry_mask
from rasterio.windows import Window, from_bounds
from rasterio.enums import Resampling
from shapely.geometry import box
from pyproj import Geod
import matplotlib
matplotlib.use("Agg")   # render to files only; no GUI backend
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, BoundaryNorm
from matplotlib.patches import Patch
from matplotlib.lines import Line2D

warnings.filterwarnings("ignore")
plt.rcParams.update({"figure.dpi": 150, "savefig.bbox": "tight", "font.size": 9,
                     "figure.facecolor": "white", "axes.facecolor": "white"})

DATA = Path("data"); DATA.mkdir(exist_ok=True)
OUT = Path("output"); OUT.mkdir(exist_ok=True)
SEED = 42
rng = np.random.default_rng(SEED)
N_BOOT = 1000

WC_URL = "https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/ESA_WorldCover_10m_2021_v200_{}_Map.tif"
LUCAS_URL = "https://ec.europa.eu/eurostat/documents/205002/17561401/IE_LUCAS_2022.csv"
GRID_URL = "https://ec.europa.eu/eurostat/documents/d/lucas/lucas-master-grid-2016"
NUTS_URL = "https://gisco-services.ec.europa.eu/distribution/v2/nuts/geojson/NUTS_RG_01M_2021_4326_LEVL_0.geojson"

# ESA WorldCover legend: code -> (short name, official colour)
WC = {10: ("Tree cover", "#006400"), 20: ("Shrubland", "#ffbb22"), 30: ("Grassland", "#ffff4c"),
      40: ("Cropland", "#f096ff"), 50: ("Built-up", "#fa0000"), 60: ("Bare / sparse", "#b4b4b4"),
      70: ("Snow and ice", "#f0f0f0"), 80: ("Water", "#0064c8"), 90: ("Herbaceous wetland", "#0096a0"),
      95: ("Mangroves", "#00cf75"), 100: ("Moss and lichen", "#fae6a0")}


def download(url, path):
    if not path.exists():
        r = requests.get(url, timeout=900, headers={"User-Agent": "TechnicalGIS-tutorial"})
        r.raise_for_status()
        path.write_bytes(r.content)
    return path


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


section("versions")
print(f"python {sys.version.split()[0]} | numpy {version('numpy')} | pandas {version('pandas')} | "
      f"rasterio {version('rasterio')} | geopandas {version('geopandas')} | matplotlib {version('matplotlib')}")

# ---------------------------------------------------------------- Step 1: the map and the study area
section("study area and WorldCover tiles")
nuts = gpd.read_file(download(NUTS_URL, DATA / "NUTS_RG_01M_2021_4326_LEVL_0.geojson"))
ireland = nuts[nuts.NUTS_ID == "IE"].to_crs("EPSG:4326")
ie_geom = ireland.geometry.iloc[0]
x0, y0, x1, y1 = ie_geom.bounds
tiles = [f"N{lat:02d}W{-lon:03d}" for lat in range(int(np.floor(y0 / 3) * 3), int(y1) + 1, 3)
         for lon in range(int(np.floor(x0 / 3) * 3), int(x1) + 1, 3)]
tiles = [t for t in tiles if box(-int(t[4:]), int(t[1:3]), -int(t[4:]) + 3, int(t[1:3]) + 3).intersects(ie_geom)]
paths = {t: download(WC_URL.format(t), DATA / f"WorldCover_2021_{t}.tif") for t in tiles}
geod = Geod(ellps="WGS84")
print(f"Ireland (NUTS IE, 1:1M): bounds {x0:.2f}..{x1:.2f} E, {y0:.2f}..{y1:.2f} N, "
      f"geodesic area {abs(geod.geometry_area_perimeter(ie_geom)[0]) / 1e6:,.0f} km2")
print(f"WorldCover 2021 v200 tiles: {', '.join(tiles)}")


# ---------------------------------------------------------------- Step 2: map class areas (pixel counting)
def row_areas(transform, row0, nrows):
    """Ellipsoidal area (m2) of one pixel in each row: 10 m pixels in EPSG:4326 shrink towards the pole."""
    lat_top = transform.f + transform.e * np.arange(row0, row0 + nrows)
    lat_bot = lat_top + transform.e
    lon0, lon1 = 0.0, transform.a
    return np.array([abs(geod.polygon_area_perimeter([lon0, lon1, lon1, lon0], [t, t, b, b])[0])
                     for t, b in zip(lat_top, lat_bot)])


section("map class areas from pixel counts")
STEP = 3000                                     # rows per window: 3000 x 36000 pixels
area_m2 = np.zeros(256)
pix = np.zeros(256, dtype=np.int64)
win_counts = []                                  # per-window class counts, for sampling in Step 8
for t, p in paths.items():
    with rasterio.open(p) as src:
        for r0 in range(0, src.height, STEP):
            win = Window(0, r0, src.width, min(STEP, src.height - r0))
            wb = rasterio.windows.bounds(win, src.transform)
            part = ie_geom.intersection(box(*wb))
            if part.is_empty:
                continue
            wt = rasterio.windows.transform(win, src.transform)
            inside = geometry_mask([part], out_shape=(win.height, win.width), transform=wt, invert=True)
            arr = src.read(1, window=win)
            rows = np.nonzero(inside)[0]
            cls = arr[inside]
            counts = np.bincount(rows * 256 + cls, minlength=win.height * 256).reshape(win.height, 256)
            area_m2 += (counts * row_areas(src.transform, r0, win.height)[:, None]).sum(0)
            pix += counts.sum(0)
            win_counts.append((t, r0, win.height, counts.sum(0)))
# Moss and lichen covers 199 pixels (0.02 km2) in Ireland, too rare to assess: merge it into bare / sparse
area_m2[60] += area_m2[100]; pix[60] += pix[100]; area_m2[100] = pix[100] = 0
for _, _, _, cnt in win_counts:
    cnt[60] += cnt[100]; cnt[100] = 0
classes = [c for c in WC if pix[c] > 0]
A_total = area_m2[classes].sum()
map_area = pd.DataFrame({"class": [WC[c][0] for c in classes], "pixels": pix[classes],
                         "area_km2": area_m2[classes] / 1e6}, index=classes)
map_area["W"] = map_area.area_km2 / map_area.area_km2.sum()
print(f"{pix[classes].sum():,} pixels inside Ireland; mapped area {A_total / 1e6:,.0f} km2 "
      f"(nodata pixels: {pix[0]:,})")
print(map_area.to_string(formatters={"pixels": "{:,}".format, "area_km2": "{:,.1f}".format, "W": "{:.4f}".format}))

# ---------------------------------------------------------------- Step 3: LUCAS 2022 reference points
section("LUCAS 2022 reference data")
lucas = pd.read_csv(download(LUCAS_URL, DATA / "IE_LUCAS_2022.csv"), low_memory=False)
obs_names = {1: "field, <= 100 m", 2: "field, > 100 m", 3: "photo-interpreted in the field",
             7: "photo-interpreted in the office"}
sdate = pd.to_datetime(lucas.SURVEY_DATE, format="%m/%d/%Y %H:%M")
print(f"{len(lucas):,} LUCAS 2022 points in Ireland, surveyed {sdate.min():%d %b %Y} .. {sdate.max():%d %b %Y}")
for k, v in lucas.SURVEY_OBS_TYPE.value_counts().sort_index().items():
    print(f"  observation type {k} ({obs_names.get(k, '?')}): {v:,}")


def lucas_to_worldcover(lc1, lu1):
    """Crosswalk from the LUCAS 2022 land cover legend (LC1) to the ESA WorldCover legend."""
    a, two = lc1[0], lc1[:2]
    if a == "A":
        return 50                                 # roofed and non-roofed artificial surfaces
    if a == "B":
        if lc1 == "B55":
            return 30                             # temporary grassland (sown ley): WorldCover sees grass
        if two == "B7":
            return 10                             # fruit trees
        if two == "B8" or lc1 == "Bx2":
            return 20                             # vineyards, nurseries and other permanent woody crops
        return 40                                 # arable crops, incl. Bx1 (arable land, photo-interpreted)
    if a == "C":
        return 10                                 # broadleaved, coniferous and mixed woodland
    if a == "D":
        return 20
    if a == "E":
        return 30
    if a == "F":
        if lc1 == "F40" and str(lu1).startswith("U11"):
            return 40                             # bare arable soil between crops (agricultural use)
        return 60                                 # rocks, sand, lichens and moss (F30), other bare soil
    if a == "G":
        return 70 if lc1 == "G50" else 80
    if a == "H":
        return 60 if lc1 in ("H22", "H23") else 90   # salines and intertidal flats are unvegetated
    raise ValueError(lc1)


lucas["ref"] = [lucas_to_worldcover(a, b) for a, b in zip(lucas.SURVEY_LC1, lucas.SURVEY_LU1)]
xwalk = lucas.groupby(["SURVEY_LC1"]).ref.first().reset_index()
print("crosswalk LUCAS LC1 -> WorldCover:")
for code, grp in xwalk.groupby("ref"):
    print(f"  {WC[code][0]:<18} <- {' '.join(grp.SURVEY_LC1)}")
moved = lucas[(lucas.SURVEY_LC1 == "F40") & (lucas.ref == 40)]
print(f"F40 bare soil on agricultural land counted as cropland: {len(moved)}; "
      f"B55 temporary grassland counted as grassland: {(lucas.SURVEY_LC1 == 'B55').sum()}")

# ---------------------------------------------------------------- Step 4: sampling design weights
section("LUCAS 2022 design weights from the master grid")
grid_csv = DATA / "lucas_master_grid_2016_IE.csv"
if not grid_csv.exists():
    import openpyxl
    xlsx = download(GRID_URL, DATA / "lucas_master_grid_2016_v2.xlsx")
    ws = openpyxl.load_workbook(xlsx, read_only=True)["Master_web"]
    rows = ws.iter_rows(values_only=True)
    hdr = next(rows)
    keep = [r for r in rows if r[1] == "IE"]
    pd.DataFrame(keep, columns=hdr).to_csv(grid_csv, index=False)
grid = pd.read_csv(grid_csv)
grid = grid[grid.NUTS0_21 == "IE"]
grid["cell"] = grid.NUTS2_21 + "_S" + grid.STR18.astype(str)
cells = grid.groupby("cell").agg(N=("POINT_ID", "size"), n=("LUCAS_FLAG", "sum"))
print(f"master grid: {len(grid):,} points on the 2 km grid in Ireland, {grid.LUCAS_FLAG.sum():,} in the 2022 sample")
str18 = {1: "arable", 2: "permanent crops", 3: "grass", 4: "wooded", 5: "shrubs", 6: "bare/rare vegetation",
         7: "artificial", 8: "inland water", 9: "transitional/coastal water", 10: "impossible to PI"}
rates = grid.groupby("STR18").agg(N=("POINT_ID", "size"), n=("LUCAS_FLAG", "sum")).astype(int)
rates["rate"] = rates.n / rates.N
print("selection rate by 2018 photo-interpretation stratum (STR18):")
for s, rr in rates.iterrows():
    print(f"  {s:>2} {str18[s]:<27} {int(rr.n):>5,} of {int(rr.N):>5,}  ({rr.rate:.0%})")
cells["w"] = cells.N / cells.n.replace(0, np.nan)
lucas = lucas.merge(grid[["POINT_ID", "cell", "STR18", "LUCAS_FLAG"]], on="POINT_ID", how="left")
lucas["w"] = lucas.cell.map(cells.w)
print(f"design cells (NUTS 2 x STR18): {len(cells)}; weights {lucas.w.min():.2f} .. {lucas.w.max():.2f}; "
      f"sum of weights {lucas.w.sum():,.0f} = grid points represented "
      f"(x 4 km2 = {lucas.w.sum() * 4:,.0f} km2)")
print(f"grid points in cells with no sample (not represented): {int(cells.N[cells.n == 0].sum())}")

# ---------------------------------------------------------------- Step 5: read the map at each point
section("map labels at the LUCAS points")
lucas["tile"] = [f"N{int(np.floor(la / 3) * 3):02d}W{int(-np.floor(lo / 3) * 3):03d}"
                 for la, lo in zip(lucas.POINT_LAT, lucas.POINT_LONG)]
lucas["mapc"] = 0
lucas["pure3x3"] = 0.0
for t, grp in lucas.groupby("tile"):
    with rasterio.open(paths[t]) as src:
        rc = [src.index(lo, la) for lo, la in zip(grp.POINT_LONG, grp.POINT_LAT)]
        maps, pure = [], []
        for r, c in rc:
            w = src.read(1, window=Window(c - 1, r - 1, 3, 3), boundless=True, fill_value=0)
            maps.append(w[1, 1]); pure.append((w == w[1, 1]).mean())
        lucas.loc[grp.index, "mapc"] = maps
        lucas.loc[grp.index, "pure3x3"] = pure
lucas.loc[lucas.mapc == 100, "mapc"] = 60
pts = lucas[lucas.mapc.isin(classes) & lucas.w.notna()].copy()
print(f"points used: {len(pts):,} (dropped {len(lucas) - len(pts)} with map nodata or no design weight)")
print(f"points whose 3x3 pixel block is all one class: {(pts.pure3x3 == 1).mean():.1%}")
K = len(classes)
idx = {c: i for i, c in enumerate(classes)}
m = pts.mapc.map(idx).to_numpy()
r = pts.ref.map(idx).to_numpy()
w = pts.w.to_numpy()
W = map_area.W.to_numpy()
names = [WC[c][0] for c in classes]

# ---------------------------------------------------------------- Step 6: the confusion matrix
section("sample confusion matrix (counts; rows = map, columns = reference)")
cm = np.bincount(m * K + r, minlength=K * K).reshape(K, K)
cm_df = pd.DataFrame(cm, index=names, columns=names)
print(cm_df.to_string())
oa_naive = np.trace(cm) / cm.sum()
ua_naive = np.diag(cm) / cm.sum(1)
pa_naive = np.diag(cm) / cm.sum(0)
print(f"\nnaive overall accuracy (count of agreeing points / all points): {oa_naive:.3f}")
pe = (cm.sum(1) * cm.sum(0)).sum() / cm.sum() ** 2
kappa = (oa_naive - pe) / (1 - pe)
print(f"Cohen's kappa: {kappa:.3f} (chance agreement {pe:.3f})")


def estimate(m, r, w, W, K):
    """Olofsson et al. (2014) stratified estimators with the map classes as (post-)strata.
    Within each map class the reference proportions are design-weighted ratio estimates."""
    sw = np.bincount(m, weights=w, minlength=K)
    nij = np.bincount(m * K + r, weights=w, minlength=K * K).reshape(K, K)
    prop = np.divide(nij, sw[:, None], out=np.zeros_like(nij), where=sw[:, None] > 0)  # P(ref j | map i)
    p = W[:, None] * prop                          # area-proportion error matrix
    oa = np.trace(p)
    ua = np.diag(p) / p.sum(1)
    pa = np.diag(p) / p.sum(0)
    return p, oa, ua, pa, p.sum(0), prop


def kish(m, w, K):
    return np.bincount(m, weights=w, minlength=K) ** 2 / np.bincount(m, weights=w ** 2, minlength=K)


def olofsson_se(prop, W, n, p):
    """Standard errors from Olofsson et al. (2014), equations 5 to 7 and 10, with n the per-stratum sample size.
    A stratum with a single point has no variance estimate; n - 1 is floored at 1 so it contributes ~0."""
    K = len(W)
    n = np.maximum(n, 2.0)
    ua = np.diag(prop)
    se_oa = np.sqrt(np.sum(W ** 2 * ua * (1 - ua) / (n - 1)))
    se_ua = np.sqrt(ua * (1 - ua) / (n - 1))
    se_area = np.sqrt(((W[:, None] ** 2) * prop * (1 - prop) / (n[:, None] - 1)).sum(0))
    Nj = p.sum(0)
    pa = np.diag(p) / Nj
    se_pa = np.zeros(K)
    for j in range(K):
        other = sum(W[i] ** 2 * prop[i, j] * (1 - prop[i, j]) / (n[i] - 1) for i in range(K) if i != j)
        se_pa[j] = np.sqrt((W[j] ** 2 * (1 - pa[j]) ** 2 * ua[j] * (1 - ua[j]) / (n[j] - 1)
                            + pa[j] ** 2 * other) / Nj[j] ** 2)
    return se_oa, se_ua, se_pa, se_area


# ---------------------------------------------------------------- Step 7: area-adjusted estimates
section("area-adjusted error matrix (estimated % of Ireland; rows = map, columns = reference)")
p, oa, ua, pa, pref, prop = estimate(m, r, w, W, K)
print(pd.DataFrame(100 * p, index=names, columns=names).round(2).to_string())
n_i = np.bincount(m, minlength=K)
n_eff = kish(m, w, K)
se_oa, se_ua, se_pa, se_area = olofsson_se(prop, W, n_eff, p)

# stratified bootstrap over the LUCAS design cells, as a check on the analytical standard errors
cell_codes, cell_idx = np.unique(pts.cell.to_numpy(), return_inverse=True)
members = [np.nonzero(cell_idx == k)[0] for k in range(len(cell_codes))]
boot = {"oa": [], "ua": [], "pa": [], "area": []}
for b in range(N_BOOT):
    s = np.concatenate([rng.choice(g, size=len(g), replace=True) for g in members])
    _, oab, uab, pab, prb, _ = estimate(m[s], r[s], w[s], W, K)
    boot["oa"].append(oab); boot["ua"].append(uab); boot["pa"].append(pab); boot["area"].append(prb)
boot = {k: np.array(v) for k, v in boot.items()}

# what you get if you ignore LUCAS's unequal selection rates
p_u, oa_u, ua_u, pa_u, pref_u, _ = estimate(m, r, np.ones_like(w), W, K)

A_km2 = A_total / 1e6
res = pd.DataFrame({
    "class": names, "n": n_i, "n_eff": n_eff.round(0).astype(int),
    "UA_naive": ua_naive, "UA": ua, "UA_ci": 1.96 * se_ua,
    "PA_naive": pa_naive, "PA": pa, "PA_ci": 1.96 * se_pa,
    "map_km2": W * A_km2, "adj_km2": pref * A_km2, "adj_ci_km2": 1.96 * se_area * A_km2,
    "boot_ci_km2": 1.96 * np.nanstd(boot["area"], 0) * A_km2, "unweighted_km2": pref_u * A_km2,
})
fmt = {"UA_naive": "{:.3f}".format, "UA": "{:.3f}".format, "UA_ci": "{:.3f}".format,
       "PA_naive": "{:.3f}".format, "PA": "{:.3f}".format, "PA_ci": "{:.3f}".format,
       "map_km2": "{:,.0f}".format, "adj_km2": "{:,.0f}".format, "adj_ci_km2": "{:,.0f}".format,
       "boot_ci_km2": "{:,.0f}".format, "unweighted_km2": "{:,.0f}".format}
print()
print(res.to_string(index=False, formatters=fmt))
print(f"\noverall accuracy (area-weighted, design-weighted): {oa:.3f} +/- {1.96 * se_oa:.3f} "
      f"(bootstrap +/- {1.96 * np.nanstd(boot['oa']):.3f}); ignoring the design weights: {oa_u:.3f}; "
      f"naive count-based: {oa_naive:.3f}")
q_dis = 0.5 * np.abs(p.sum(0) - p.sum(1)).sum()
a_dis = (1 - oa) - q_dis
pe_p = (p.sum(0) * p.sum(1)).sum()
print(f"Pontius & Millones: quantity disagreement {q_dis:.3f}, allocation disagreement {a_dis:.3f} "
      f"(total {1 - oa:.3f}); kappa on the area-adjusted matrix {(oa - pe_p) / (1 - pe_p):.3f}")
res.to_csv(OUT / "accuracy_and_area.csv", index=False)
print(f"note: map shrubland has n = {n_i[names.index('Shrubland')]} point, so its user's accuracy is not estimable")

# ---------------------------------------------------------------- Step 8: point vs pixel and other caveats
section("sensitivity: observation type, homogeneity and the Copernicus module")


def group_oa(mask):
    return estimate(m[mask], r[mask], w[mask], W, K)[1], int(mask.sum())


sens = []
ot = pts.SURVEY_OBS_TYPE.to_numpy()
for k, label in [(1, "Field, <= 100 m"), (2, "Field, > 100 m"), (3, "PI in the field"), (7, "PI in the office")]:
    sens.append(("Observation type", label, *group_oa(ot == k)))
pure = pts.pure3x3.to_numpy()
sens.append(("3x3 map window", "All 9 pixels agree", *group_oa(pure == 1)))
sens.append(("3x3 map window", "Mixed", *group_oa(pure < 1)))
ext = pts[["SURVEY_CPRN_LC1N", "SURVEY_CPRNC_LC1E", "SURVEY_CPRNC_LC1S", "SURVEY_CPRNC_LC1W"]].min(1).to_numpy()
cop = ~np.isnan(ext)
sens.append(("Copernicus module", "Same cover >= 20 m all round", *group_oa(cop & (ext >= 20))))
sens.append(("Copernicus module", "Same cover < 20 m", *group_oa(cop & (ext < 20))))
sens = pd.DataFrame(sens, columns=["group", "subset", "OA", "n"])
print(sens.to_string(index=False, formatters={"OA": "{:.3f}".format}))
year = pd.to_datetime(pts.SURVEY_DATE, format="%m/%d/%Y %H:%M").dt.year
print(f"survey year: 2022 {int((year == 2022).sum()):,}, 2023 {int((year == 2023).sum()):,} points "
      f"(WorldCover 2021 describes 2021)")
grass_as_crop = ((pts.mapc == 40) & (pts.ref == 30))
print(f"map cropland, LUCAS grassland: {grass_as_crop.sum()} points, of which B55 temporary grassland "
      f"{((pts.SURVEY_LC1 == 'B55') & grass_as_crop).sum()}")
bog = pts.SURVEY_LC1 == "H12"
print("LUCAS peat bog (H12) points by map class: " +
      ", ".join(f"{WC[c][0]} {v}" for c, v in pts[bog].mapc.value_counts().items()))

# ---------------------------------------------------------------- Step 9: sample-size planning
section("sample size for a new stratified random sample of the map")
# conjectured user's accuracies from this assessment; 0.5 (the worst case) where n < 30 points
U_plan = np.where(n_i >= 30, np.clip(ua, 0.5, 0.95), 0.5)
print("planning user's accuracies: " + ", ".join(f"{names[i]} {U_plan[i]:.2f}" for i in range(K)))
S = np.sqrt(U_plan * (1 - U_plan))
target_se = 0.01
n_total = int(np.ceil((np.sum(W * S) / target_se) ** 2))
print(f"Cochran (1977) eq. 5.25 with target SE(OA) = {target_se}: n = {n_total}")


def expected_se(n_alloc):
    """Expected SE of OA and of each UA for an allocation, from the conjectured user's accuracies."""
    v = U_plan * (1 - U_plan) / (np.maximum(n_alloc, 2) - 1)
    return np.sqrt(np.sum(W ** 2 * v)), np.sqrt(v)


allocs = {}
allocs["Proportional"] = np.maximum(np.round(n_total * W), 2)
allocs["Equal"] = np.full(K, np.round(n_total / K))
comp = np.full(K, 100.0)                       # at least 100 per class, the rest in proportion to area
rest = n_total - comp.sum()
big = W * n_total > 100
comp[big] = np.round(100 + rest * W[big] / W[big].sum())
allocs["Compromise (min 100)"] = comp
for name, a in allocs.items():
    so, su = expected_se(a)
    print(f"{name:<21} n = {int(a.sum()):>5}: SE(OA) {so:.4f}; SE(UA) " +
          ", ".join(f"{names[i].split()[0]} {su[i]:.3f}" for i in range(K)))
print("compromise allocation: " + ", ".join(f"{names[i]} {int(comp[i])}" for i in range(K)))

# draw the compromise sample: a stratified random sample of map pixels, without replacement
wc_df = pd.DataFrame([(t, r0, h, *cnt[classes]) for t, r0, h, cnt in win_counts],
                     columns=["tile", "row0", "rows", *classes])
draws = []
for i, c in enumerate(classes):
    per_win = rng.multivariate_hypergeometric(wc_df[c].to_numpy(), int(comp[i]))
    for (t, r0, h), k in zip(wc_df[["tile", "row0", "rows"]].itertuples(index=False), per_win):
        if k:
            draws.append((t, r0, h, c, k))
sample = []
for (t, r0, h), grp in pd.DataFrame(draws, columns=["t", "r0", "h", "c", "k"]).groupby(["t", "r0", "h"]):
    with rasterio.open(paths[t]) as src:
        win = Window(0, r0, src.width, h)
        wt = rasterio.windows.transform(win, src.transform)
        part = ie_geom.intersection(box(*rasterio.windows.bounds(win, src.transform)))
        inside = geometry_mask([part], out_shape=(h, src.width), transform=wt, invert=True)
        arr = src.read(1, window=win)
        for c, k in zip(grp.c, grp.k):
            rr, cc = np.nonzero(inside & (arr == c))
            pick = rng.choice(len(rr), size=k, replace=False)
            xs, ys = rasterio.transform.xy(wt, rr[pick], cc[pick])
            sample += [(c, WC[c][0], x, y) for x, y in zip(xs, ys)]
new = gpd.GeoDataFrame(pd.DataFrame(sample, columns=["map_code", "map_class", "lon", "lat"]),
                       geometry=gpd.points_from_xy([s[2] for s in sample], [s[3] for s in sample]), crs="EPSG:4326")
new["reference"] = ""                          # to be filled in by an interpreter, blind to map_class
new = new.sample(frac=1, random_state=SEED).reset_index(drop=True)   # shuffle so interpreters can't guess
new.to_file(OUT / "worldcover_ireland_validation_sample.gpkg", driver="GPKG")
print(f"wrote {len(new):,} stratified random sample points to output/worldcover_ireland_validation_sample.gpkg")
print(new.map_class.value_counts().to_string())

# ---------------------------------------------------------------- Figures
section("figures")
cmap = ListedColormap([WC[c][1] for c in classes])
norm = BoundaryNorm(np.arange(K + 1) - 0.5, K)

# Figure 1: the map and the reference points
res_deg = 1 / 400                                    # 250 m: majority (mode) of 25 x 25 pixels
mx0, my0, mx1, my1 = -10.75, 51.35, -5.85, 55.45
W_px, H_px = int(round((mx1 - mx0) / res_deg)), int(round((my1 - my0) / res_deg))
mosaic = np.zeros((H_px, W_px), dtype=np.uint8)
mtrans = rasterio.transform.from_origin(mx0, my1, res_deg, res_deg)
for t, pth in paths.items():
    with rasterio.open(pth) as src:
        b = src.bounds
        ix0, ix1 = max(b.left, mx0), min(b.right, mx1)
        iy0, iy1 = max(b.bottom, my0), min(b.top, my1)
        if ix0 >= ix1 or iy0 >= iy1:
            continue
        sw = from_bounds(ix0, iy0, ix1, iy1, src.transform).round_offsets().round_lengths()
        dw = from_bounds(ix0, iy0, ix1, iy1, mtrans).round_offsets().round_lengths()
        mosaic[dw.row_off:dw.row_off + dw.height, dw.col_off:dw.col_off + dw.width] = src.read(
            1, window=sw, out_shape=(dw.height, dw.width), resampling=Resampling.mode)
inside_lr = geometry_mask([ie_geom], out_shape=mosaic.shape, transform=mtrans, invert=True)
lut = np.full(256, np.nan); lut[classes] = np.arange(K)
img = np.where(inside_lr, lut[mosaic], np.nan)
fig, axes = plt.subplots(1, 2, figsize=(11, 6.6), constrained_layout=True)
ext_ = (mx0, mx1, my0, my1)
axes[0].imshow(img, cmap=cmap, norm=norm, extent=ext_, interpolation="nearest")
axes[0].set_title("ESA WorldCover 2021 (10 m, shown at 250 m)")
handles = [Patch(facecolor=WC[c][1], edgecolor="0.4", lw=0.4, label=f"{WC[c][0]} ({100 * map_area.W[c]:.1f}%)")
           for c in classes]
axes[0].legend(handles=handles, loc="lower right", fontsize=7, title="Map class (share of area)",
               title_fontsize=7, framealpha=0.95)
axes[1].imshow(np.where(inside_lr, 0.0, np.nan), cmap=ListedColormap(["#e6e6e6"]), extent=ext_)
agree = pts.mapc == pts.ref
for ok, col, lab in [(True, "#0072B2", "Map agrees with LUCAS"), (False, "#D55E00", "Map disagrees")]:
    s = pts[agree == ok]
    axes[1].scatter(s.POINT_LONG, s.POINT_LAT, s=1.6, c=col, lw=0, label=f"{lab} ({len(s):,})")
axes[1].set_title(f"{len(pts):,} LUCAS 2022 survey points")
axes[1].legend(loc="lower right", markerscale=5, fontsize=7.5, framealpha=0.95)
for ax in axes:
    ax.set_xlim(mx0, mx1); ax.set_ylim(my0, my1); ax.set_aspect(1 / np.cos(np.radians(53.4)))
    ax.set_xlabel("Longitude (°)"); ax.set_ylabel("Latitude (°)")
fig.savefig(OUT / "01-worldcover-lucas-map.png"); plt.close(fig)

# Figure 2: confusion matrix (counts) and area-adjusted error matrix
fig, axes = plt.subplots(1, 2, figsize=(12.5, 5.6), constrained_layout=True)
short = [n.replace("Herbaceous wetland", "Herb. wetland").replace("Moss and lichen", "Moss/lichen") for n in names]
for ax, mat, title, fmtf, cm_name in [
        (axes[0], cm, "Sample counts (LUCAS points)", lambda v: f"{int(v):,}", "Blues"),
        (axes[1], 100 * p, "Area-adjusted error matrix (% of Ireland)", lambda v: f"{v:.2f}", "Blues")]:
    im = ax.imshow(np.log10(mat + (1 if mat is cm else 0.01)), cmap=cm_name)
    for i in range(K):
        for j in range(K):
            v = mat[i, j]
            if v > 0:
                dark = im.norm(np.log10(v + (1 if mat is cm else 0.01))) > 0.6
                ax.text(j, i, fmtf(v), ha="center", va="center", fontsize=6.5,
                        color="white" if dark else "black", fontweight="bold" if i == j else None)
    ax.set_xticks(range(K), short, rotation=45, ha="right"); ax.set_yticks(range(K), short)
    ax.set_xlabel("Reference class (LUCAS 2022)"); ax.set_ylabel("Map class (WorldCover 2021)")
    ax.set_title(title)
    cb = fig.colorbar(im, ax=ax, shrink=0.8)
    cb.set_label("log10(points + 1)" if mat is cm else "log10(% of area + 0.01)")
fig.savefig(OUT / "02-confusion-matrices.png"); plt.close(fig)

# Figure 3: user's and producer's accuracy, naive vs stratified with 95% CIs
fig, axes = plt.subplots(1, 2, figsize=(12, 4.8), sharey=True, constrained_layout=True)
y = np.arange(K)
ok = n_i >= 5
for ax, est, ci, naive, title in [(axes[0], np.where(ok, ua, np.nan), 1.96 * se_ua, np.where(ok, ua_naive, np.nan), "User's accuracy (map reliability)"),
                                  (axes[1], pa, 1.96 * se_pa, pa_naive, "Producer's accuracy (omission)")]:
    lo, hi = np.clip(est - ci, 0, 1), np.clip(est + ci, 0, 1)
    ax.errorbar(est, y, xerr=[est - lo, hi - est], fmt="o", color="#0072B2", ms=5, capsize=3,
                label="Stratified estimate (Olofsson) with 95% CI")
    ax.scatter(naive, y + 0.25, marker="x", color="#D55E00", s=28, label="Naive (sample counts)")
    ax.set_xlim(0, 1.02); ax.set_title(title); ax.set_xlabel("Accuracy")
    ax.axvline(oa, color="0.5", ls="--", lw=0.8)
    ax.grid(axis="x", color="0.9")
axes[0].set_yticks(y, [f"{nm} (n={n_i[i]})" for i, nm in enumerate(names)]); axes[0].invert_yaxis()
h, l = axes[1].get_legend_handles_labels()
fig.legend(h, l, loc="outside lower center", ncols=2, fontsize=8)
axes[0].text(oa + 0.01, K - 0.6, f"OA {oa:.3f}", color="0.35", fontsize=7.5)
for i in np.nonzero(~ok)[0]:
    axes[0].text(0.02, i, f"not estimable (n = {n_i[i]})", va="center", fontsize=7.5, color="0.35")
fig.savefig(OUT / "03-users-producers-accuracy.png"); plt.close(fig)

# Figure 4: area estimates
fig, ax = plt.subplots(figsize=(9, 4.8), constrained_layout=True)
order = np.argsort(-res.map_km2.to_numpy())
yy = np.arange(K)
ax.barh(yy - 0.2, res.map_km2.to_numpy()[order], height=0.38, color="#bbbbbb", label="Pixel counting (map area)")
ax.barh(yy + 0.2, res.adj_km2.to_numpy()[order], height=0.38, color="#0072B2",
        xerr=res.adj_ci_km2.to_numpy()[order], capsize=3, error_kw={"lw": 1, "ecolor": "black"},
        label="Error-adjusted estimate with 95% CI")
ax.set_yticks(yy, [names[i] for i in order]); ax.invert_yaxis()
ax.set_xscale("log"); ax.set_xlabel("Area (km², log scale)")
for k, i in enumerate(order):
    ax.text(res.map_km2[i] * 1.12, k - 0.2, f"{res.map_km2[i]:,.0f}", va="center", fontsize=7, color="0.3")
    ax.text((res.adj_km2[i] + res.adj_ci_km2[i]) * 1.12, k + 0.2,
            f"{res.adj_km2[i]:,.0f} ± {res.adj_ci_km2[i]:,.0f}", va="center", fontsize=7)
ax.set_xlim(0.5, 3e5)
fig.legend(loc="outside lower center", ncols=2, fontsize=8); ax.grid(axis="x", color="0.9", which="both"); ax.set_axisbelow(True)
ax.set_title("Area of each land cover class in Ireland: map vs error-adjusted")
fig.savefig(OUT / "04-area-estimates.png"); plt.close(fig)

# Figure 5: sensitivity to observation type, map homogeneity and reference homogeneity
fig, ax = plt.subplots(figsize=(8.5, 4.2), constrained_layout=True)
cols = {"Observation type": "#56B4E9", "3x3 map window": "#E69F00", "Copernicus module": "#009E73"}
yy = np.arange(len(sens))
ax.barh(yy, sens.OA, color=[cols[g] for g in sens.group])
for i, rw in sens.iterrows():
    ax.text(rw.OA + 0.005, i, f"{rw.OA:.3f}  (n={rw.n:,})", va="center", fontsize=7.5)
ax.set_yticks(yy, sens.subset); ax.invert_yaxis(); ax.set_xlim(0.4, 1.0)
ax.axvline(oa, color="0.4", ls="--", lw=0.8)
ax.set_xlabel("Overall accuracy (area- and design-weighted)")
ax.legend(handles=[Patch(color=c, label=g) for g, c in cols.items()] +
          [Line2D([], [], color="0.4", ls="--", lw=0.8, label=f"All points ({oa:.3f})")], loc="lower right", fontsize=7.5)
fig.savefig(OUT / "05-point-vs-pixel.png"); plt.close(fig)

# Figure 6: sample-size planning
fig, axes = plt.subplots(1, 2, figsize=(12, 4.6), constrained_layout=True)
ns = np.arange(100, 3001, 50)
axes[0].plot(ns, np.sum(W * S) / np.sqrt(ns), color="#0072B2")
axes[0].axhline(target_se, color="0.5", ls="--", lw=0.8)
axes[0].axvline(n_total, color="#D55E00", lw=0.8)
axes[0].text(n_total + 40, 0.025, f"n = {n_total}", color="#D55E00")
axes[0].set_xlabel("Total sample size n"); axes[0].set_ylabel("SE of overall accuracy")
axes[0].set_title("Cochran's formula: SE(OA) ≈ Σ W$_i$ S$_i$ / √n"); axes[0].grid(color="0.9")
bw = 0.27
pal = ["#bbbbbb", "#56B4E9", "#0072B2"]
for k, (name, a) in enumerate(allocs.items()):
    su = expected_se(a)[1]
    axes[1].bar(np.arange(K) + (k - 1) * bw, 1.96 * su, width=bw, color=pal[k], label=f"{name}")
    for i in np.nonzero(1.96 * su > 0.3)[0]:
        axes[1].text(i + (k - 1) * bw, 0.285, f"{1.96 * su[i]:.2f}", rotation=90, ha="center", va="top", fontsize=7)
axes[1].set_xticks(range(K), short, rotation=45, ha="right")
axes[1].set_ylabel("Expected 95% CI half-width of user's accuracy")
axes[1].set_title(f"Allocating n = {n_total} between the map classes"); axes[1].legend(fontsize=7.5)
axes[1].set_ylim(0, 0.3); axes[1].grid(axis="y", color="0.9"); axes[1].set_axisbelow(True)
fig.savefig(OUT / "06-sample-size-planning.png"); plt.close(fig)
print("wrote 6 figures to output/")