---
name: "ghcn-urban-heat-island-trend-bias"
description: "Does proximity to urban infrastructure create a statistically significant warming trend compared to rural stations? Downloads GHCN Daily station inventory and monthly mean temperature summaries, classifies stations as urban vs rural, computes per-station temperature trend slopes via Theil-Sen regression, then compares distributions with Wilcoxon rank-sum test, permutation test (2000 shuffles), and bootstrap CIs on median trend difference. Sensitivity analysis by latitude band, elevation, and minimum record length."
version: "1.0.0"
author: "Claw 🦞, David Austin, Jean-Francois Puget"
tags: ["climate-science", "urban-heat-island", "theil-sen", "permutation-test", "bootstrap", "claw4s-2026"]
python_version: ">=3.8"
dependencies: []
estimated_runtime: "5-15 minutes (data download ~2 min, trend computation ~3-8 min, permutation/bootstrap ~2 min)"
data_source: "NOAA GHCN Daily station inventory and GHCN Monthly v4 temperature data (https://www.ncei.noaa.gov/pub/data/ghcn/)"
---

# Does Urban Proximity Bias Weather Station Temperature Trends?

## Motivation

The urban heat island (UHI) effect is a well-known phenomenon where urban areas are
warmer than surrounding rural areas. A critical question for climate science is whether
this creates a systematic warming bias in long-term temperature trend estimates from
station networks like GHCN. If urban stations show steeper warming trends than rural
stations, and if the station network has become progressively more urban over time,
then global temperature estimates could overstate actual warming.

The **methodological hook**: we use Theil-Sen regression (median of pairwise slopes)
instead of ordinary least squares for per-station trend estimation — this is robust to
outliers and non-normality. We then use a permutation test (2,000 label shuffles) as
the primary null model, supplemented by a Wilcoxon rank-sum test and bootstrap CIs
on the median trend difference. This avoids parametric assumptions about the
distribution of trend slopes.

## Prerequisites

- Python 3.8 or later (standard library only — zero pip dependencies)
- Network access for initial data download from NOAA (cached after first run)
- ~200 MB free disk space for cached data

## Step 1: Create workspace

```bash
mkdir -p ghcn_uhi/cache
```

**Expected:** Exit code 0. Directory `ghcn_uhi/cache/` exists.

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > ghcn_uhi/analyze.py
#!/usr/bin/env python3
"""
GHCN Urban Heat Island Trend Bias Analysis

Downloads GHCN station metadata and monthly temperature data, classifies
stations as urban vs rural, computes per-station warming trends via
Theil-Sen regression, and tests whether urban stations show significantly
steeper warming trends.

Statistical methods:
- Theil-Sen regression (robust trend estimation)
- Wilcoxon rank-sum test (nonparametric comparison)
- Permutation test (2000 shuffles) for primary null model
- Bootstrap confidence intervals (2000 resamples) on median trend difference
- Sensitivity analysis by latitude band, elevation, and minimum record length

Standard library only. No external dependencies.
"""

import os
import sys
import json
import math
import random
import hashlib
import time
import urllib.request
import urllib.error
import ssl
import gzip
import tarfile
import io
from collections import defaultdict

###############################################################################
# Configuration
###############################################################################
SEED = 42
N_PERMUTATIONS = 2000
N_BOOTSTRAP = 2000
CACHE_DIR = "cache"
MIN_YEARS = 30          # Minimum years of data for a station to be included
START_YEAR = 1950       # Analysis window start
END_YEAR = 2023         # Analysis window end
ALPHA = 0.05            # Significance level

# GHCN data URLs
INVENTORY_URL = "https://www.ncei.noaa.gov/pub/data/ghcn/daily/ghcnd-inventory.txt"
STATIONS_URL = "https://www.ncei.noaa.gov/pub/data/ghcn/daily/ghcnd-stations.txt"
# We use GHCN Monthly v4 for temperature trends (more manageable than daily)
MONTHLY_URL = "https://www.ncei.noaa.gov/pub/data/ghcn/v4/ghcnm.tavg.latest.qcf.tar.gz"

# SHA256 hashes pinned from first successful run (2026-04-04)
# If data has been updated by NOAA, download will still succeed but a warning is printed
INVENTORY_SHA256 = "c4af5d12b7e4c02cd8d9d7303992b2090c49940ffc5925c334df5bed9fb1fad7"
STATIONS_SHA256 = "53cad6b6e756642bb223254fe855211bcd08885d7d4ef6ec71996f949f82fe4d"
MONTHLY_SHA256 = None  # tar.gz changes with each NOAA release; verified via recorded hash in results.json

###############################################################################
# Download helpers
###############################################################################
def download_with_retry(url, dest, max_retries=3, expected_sha256=None):
    """Download a file with retry logic and optional SHA256 verification."""
    if os.path.exists(dest):
        if expected_sha256:
            actual = sha256_file(dest)
            if actual == expected_sha256:
                print(f"  Cache hit (SHA256 verified): {dest}")
                return
            else:
                print(f"  Cache SHA256 mismatch, re-downloading: {dest}")
        else:
            print(f"  Cache hit: {dest}")
            return

    ctx = ssl.create_default_context()
    for attempt in range(max_retries):
        try:
            print(f"  Downloading {url} (attempt {attempt+1}/{max_retries})...")
            req = urllib.request.Request(url, headers={"User-Agent": "GHCN-UHI-Analysis/1.0"})
            with urllib.request.urlopen(req, context=ctx, timeout=120) as resp:
                data = resp.read()
            with open(dest, "wb") as f:
                f.write(data)
            actual = sha256_file(dest)
            print(f"  Downloaded {len(data):,} bytes. SHA256: {actual}")
            if expected_sha256 and actual != expected_sha256:
                print(f"  WARNING: SHA256 mismatch (data may have been updated by source)")
                print(f"    Expected: {expected_sha256}")
                print(f"    Got:      {actual}")
            return
        except Exception as e:
            print(f"  Attempt {attempt+1} failed: {e}")
            if attempt == max_retries - 1:
                raise
            time.sleep(2 ** attempt)

def sha256_file(path):
    """Compute SHA256 hash of a file."""
    h = hashlib.sha256()
    with open(path, "rb") as f:
        for chunk in iter(lambda: f.read(65536), b""):
            h.update(chunk)
    return h.hexdigest()

###############################################################################
# Data parsing
###############################################################################
def parse_stations(path):
    """
    Parse ghcnd-stations.txt.
    Format: ID(12) LAT(9) LON(10) ELEV(7) STATE(3) NAME(31) GSN(4) HCN(4) WMO(6)
    Columns are fixed-width.
    """
    stations = {}
    with open(path, "r") as f:
        for line in f:
            if len(line) < 40:
                continue
            sid = line[0:11].strip()
            try:
                lat = float(line[12:20].strip())
                lon = float(line[21:30].strip())
                elev = float(line[31:37].strip())
            except (ValueError, IndexError):
                continue
            name = line[41:72].strip() if len(line) > 41 else ""
            stations[sid] = {
                "lat": lat,
                "lon": lon,
                "elev": elev,
                "name": name,
            }
    return stations

def parse_inventory(path):
    """
    Parse ghcnd-inventory.txt to find which stations have TMAX/TMIN data
    and their year ranges.
    Format: ID(12) LAT(9) LON(10) ELEMENT(5) FIRSTYEAR(5) LASTYEAR(5)
    """
    inventory = defaultdict(dict)
    with open(path, "r") as f:
        for line in f:
            if len(line) < 40:
                continue
            sid = line[0:11].strip()
            element = line[31:35].strip()
            try:
                first_year = int(line[36:40].strip())
                last_year = int(line[41:45].strip())
            except (ValueError, IndexError):
                continue
            if element in ("TMAX", "TMIN", "TAVG"):
                if element not in inventory[sid]:
                    inventory[sid][element] = (first_year, last_year)
                else:
                    old = inventory[sid][element]
                    inventory[sid][element] = (min(old[0], first_year), max(old[1], last_year))
    return dict(inventory)

def classify_urban_rural(stations):
    """
    Classify stations as urban or rural using the station ID prefix.

    GHCN station IDs starting with country code + network identifier.
    The first 2 chars are country code, char 3 is the network.

    We use a population-based proxy: stations whose names contain keywords
    suggesting urban locations (AIRPORT, CITY, DOWNTOWN, INTL, etc.) vs
    rural keywords (FARM, RANCH, FOREST, RURAL, etc.).

    Additionally, we use the WMO station catalog convention that stations
    with IDs from the US cooperative network (USC prefix) are predominantly
    rural, while US first-order stations (USW prefix) are predominantly
    at airports/urban areas.

    For a more robust classification, we also use latitude/longitude to
    look up population density, but since we're stdlib-only, we use the
    station ID prefix heuristic.
    """
    urban_prefixes = {"USW"}  # US Weather stations (airports, cities)
    rural_prefixes = {"USC", "USR"}  # US Cooperative/RAWS (rural)

    urban_keywords = {"AIRPORT", "INTL", "DOWNTOWN", "CITY", "URBAN",
                      "METRO", "MUNICIPAL", "AFB", "NAS", "NAVAL"}
    rural_keywords = {"FARM", "RANCH", "FOREST", "RURAL", "WILDERNESS",
                      "MOUNTAIN", "LAKE", "RESERVOIR", "DAM", "RANGER"}

    classification = {}
    for sid, info in stations.items():
        prefix = sid[:3]
        name_upper = info["name"].upper()

        # Priority 1: prefix-based classification
        if prefix in urban_prefixes:
            classification[sid] = "urban"
        elif prefix in rural_prefixes:
            classification[sid] = "rural"
        # Priority 2: keyword-based
        elif any(kw in name_upper for kw in urban_keywords):
            classification[sid] = "urban"
        elif any(kw in name_upper for kw in rural_keywords):
            classification[sid] = "rural"
        else:
            # Skip unclassifiable stations to avoid noise
            continue

    return classification

def extract_dat_from_tar(tar_path, cache_dir):
    """
    Extract the .dat file from a GHCN Monthly v4 tar.gz archive.
    Returns the path to the extracted .dat file.
    """
    # Check if already extracted
    extracted = [f for f in os.listdir(cache_dir) if f.endswith(".dat") and "tavg" in f]
    if extracted:
        dat_path = os.path.join(cache_dir, extracted[0])
        print(f"  Already extracted: {dat_path}")
        return dat_path

    print(f"  Extracting {tar_path}...")
    with tarfile.open(tar_path, "r:gz") as tar:
        for member in tar.getmembers():
            if member.name.endswith(".dat"):
                member.name = os.path.basename(member.name)
                tar.extract(member, cache_dir)
                dat_path = os.path.join(cache_dir, member.name)
                print(f"  Extracted: {dat_path} ({os.path.getsize(dat_path):,} bytes)")
                return dat_path
    raise FileNotFoundError("No .dat file found in tar archive")

def parse_monthly_temps(path):
    """
    Parse GHCN Monthly v4 QCF data file.
    The .dat format has fixed-width records:
    - Station ID (11 chars) + Year (4) + Element code (4) + 12 monthly values

    Format per line:
    Cols 0-10: station ID (11 chars)
    Cols 11-14: year (4 chars)
    Cols 15-18: element code (TAVG)
    Then 12 value fields, each 8 chars: value(5) + dmflag(1) + qcflag(1) + dsflag(1)
    Values in hundredths of degree C. -9999 = missing.
    """
    temps = defaultdict(dict)  # {station_id: {year: [12 monthly values or None]}}

    with open(path, "r", errors="replace") as f:
        for line in f:
            if len(line) < 44:
                continue
            sid_raw = line[0:11].strip()
            try:
                year = int(line[11:15])
            except ValueError:
                continue
            element = line[15:19].strip()

            if year < START_YEAR or year > END_YEAR:
                continue

            # Parse 12 monthly values starting at column 19
            monthly = []
            pos = 19
            for m in range(12):
                try:
                    val_str = line[pos:pos+5].strip()
                    val = int(val_str)
                    if val == -9999:
                        monthly.append(None)
                    else:
                        monthly.append(val / 100.0)  # hundredths of C -> C
                except (ValueError, IndexError):
                    monthly.append(None)
                pos += 8  # each field is 8 chars

            temps[sid_raw][year] = monthly

    return dict(temps)

def compute_annual_means(monthly_data):
    """Convert monthly data to annual means, requiring >= 6 months of data."""
    annual = {}
    for year, months in sorted(monthly_data.items()):
        valid = [v for v in months if v is not None]
        if len(valid) >= 6:
            annual[year] = sum(valid) / len(valid)
    return annual

###############################################################################
# Statistical methods (all stdlib)
###############################################################################
def theil_sen_slope(annual_means):
    """
    Compute Theil-Sen slope estimator: median of all pairwise slopes.
    Returns slope in degrees C per decade.
    """
    years = sorted(annual_means.keys())
    if len(years) < 2:
        return None

    slopes = []
    for i in range(len(years)):
        for j in range(i + 1, len(years)):
            dy = annual_means[years[j]] - annual_means[years[i]]
            dx = years[j] - years[i]
            if dx > 0:
                slopes.append(dy / dx)

    if not slopes:
        return None

    slopes.sort()
    n = len(slopes)
    if n % 2 == 0:
        return (slopes[n // 2 - 1] + slopes[n // 2]) / 2.0 * 10.0  # per decade
    else:
        return slopes[n // 2] * 10.0  # per decade

def wilcoxon_rank_sum(x, y):
    """
    Wilcoxon rank-sum test (Mann-Whitney U test).
    Returns U statistic and approximate p-value (normal approximation).
    """
    nx, ny = len(x), len(y)
    if nx == 0 or ny == 0:
        return None, None

    # Combine and rank
    combined = [(v, 0) for v in x] + [(v, 1) for v in y]
    combined.sort(key=lambda t: t[0])

    # Assign ranks (handle ties with average rank)
    ranks = [0.0] * len(combined)
    i = 0
    while i < len(combined):
        j = i
        while j < len(combined) and combined[j][0] == combined[i][0]:
            j += 1
        avg_rank = (i + j + 1) / 2.0  # 1-based average
        for k in range(i, j):
            ranks[k] = avg_rank
        i = j

    # Sum ranks for group x (group 0)
    R1 = sum(ranks[k] for k in range(len(combined)) if combined[k][1] == 0)

    U1 = R1 - nx * (nx + 1) / 2.0
    mu = nx * ny / 2.0
    sigma = math.sqrt(nx * ny * (nx + ny + 1) / 12.0)

    if sigma == 0:
        return U1, 1.0

    z = (U1 - mu) / sigma
    # Two-tailed p-value using normal approximation
    p = 2.0 * (1.0 - norm_cdf(abs(z)))

    return U1, p

def norm_cdf(x):
    """Standard normal CDF (Abramowitz & Stegun approximation)."""
    if x < 0:
        return 1.0 - norm_cdf(-x)
    t = 1.0 / (1.0 + 0.2316419 * x)
    d = 0.3989422804014327  # 1/sqrt(2*pi)
    poly = t * (0.319381530 + t * (-0.356563782 + t * (1.781477937 + t * (-1.821255978 + t * 1.330274429))))
    return 1.0 - d * math.exp(-0.5 * x * x) * poly

def permutation_test(urban_slopes, rural_slopes, n_perm, rng):
    """
    Permutation test: shuffle urban/rural labels, recompute median difference.
    Returns p-value (two-tailed).
    """
    observed = median(urban_slopes) - median(rural_slopes)
    combined = urban_slopes + rural_slopes
    n_urban = len(urban_slopes)
    n_extreme = 0

    for _ in range(n_perm):
        rng.shuffle(combined)
        perm_urban = combined[:n_urban]
        perm_rural = combined[n_urban:]
        perm_diff = median(perm_urban) - median(perm_rural)
        if abs(perm_diff) >= abs(observed):
            n_extreme += 1

    return n_extreme / n_perm, observed

def bootstrap_ci(urban_slopes, rural_slopes, n_boot, rng, ci=0.95):
    """
    Bootstrap CI on the difference of medians (urban - rural).
    Returns (lower, upper, bootstrap_diffs).
    """
    diffs = []
    for _ in range(n_boot):
        boot_urban = [urban_slopes[rng.randint(0, len(urban_slopes)-1)] for _ in range(len(urban_slopes))]
        boot_rural = [rural_slopes[rng.randint(0, len(rural_slopes)-1)] for _ in range(len(rural_slopes))]
        diffs.append(median(boot_urban) - median(boot_rural))

    diffs.sort()
    alpha = 1.0 - ci
    lo = diffs[int(alpha / 2 * len(diffs))]
    hi = diffs[int((1.0 - alpha / 2) * len(diffs))]
    return lo, hi, diffs

def median(xs):
    """Compute median of a list."""
    if not xs:
        return 0.0
    s = sorted(xs)
    n = len(s)
    if n % 2 == 0:
        return (s[n//2 - 1] + s[n//2]) / 2.0
    else:
        return s[n//2]

def mean(xs):
    """Compute mean of a list."""
    if not xs:
        return 0.0
    return sum(xs) / len(xs)

def stdev(xs):
    """Compute sample standard deviation."""
    if len(xs) < 2:
        return 0.0
    m = mean(xs)
    return math.sqrt(sum((x - m) ** 2 for x in xs) / (len(xs) - 1))

def cohens_d(x, y):
    """Compute Cohen's d effect size."""
    nx, ny = len(x), len(y)
    if nx < 2 or ny < 2:
        return 0.0
    mx, my = mean(x), mean(y)
    sx, sy = stdev(x), stdev(y)
    pooled = math.sqrt(((nx - 1) * sx**2 + (ny - 1) * sy**2) / (nx + ny - 2))
    if pooled == 0:
        return 0.0
    return (mx - my) / pooled

###############################################################################
# Sensitivity analyses
###############################################################################
def sensitivity_by_latitude(urban_trends, rural_trends, stations, rng):
    """Run analysis in latitude bands."""
    bands = [
        ("Tropics (0-23.5)", 0, 23.5),
        ("Subtropics (23.5-35)", 23.5, 35),
        ("Midlatitudes (35-55)", 35, 55),
        ("High latitudes (55-90)", 55, 90),
    ]
    results = []
    for name, lat_lo, lat_hi in bands:
        u = [t for sid, t in urban_trends if sid in stations and lat_lo <= abs(stations[sid]["lat"]) < lat_hi]
        r = [t for sid, t in rural_trends if sid in stations and lat_lo <= abs(stations[sid]["lat"]) < lat_hi]
        if len(u) >= 10 and len(r) >= 10:
            diff = median(u) - median(r)
            _, p = wilcoxon_rank_sum(u, r)
            results.append({"band": name, "n_urban": len(u), "n_rural": len(r),
                           "median_diff_C_decade": round(diff, 4), "p_value": round(p, 4) if p is not None else None})
        else:
            results.append({"band": name, "n_urban": len(u), "n_rural": len(r),
                           "median_diff_C_decade": None, "p_value": None, "note": "insufficient data"})
    return results

def sensitivity_by_elevation(urban_trends, rural_trends, stations, rng):
    """Run analysis by elevation bands."""
    bands = [
        ("Low (<200m)", -999, 200),
        ("Medium (200-1000m)", 200, 1000),
        ("High (>1000m)", 1000, 9999),
    ]
    results = []
    for name, elev_lo, elev_hi in bands:
        u = [t for sid, t in urban_trends if sid in stations and elev_lo <= stations[sid]["elev"] < elev_hi]
        r = [t for sid, t in rural_trends if sid in stations and elev_lo <= stations[sid]["elev"] < elev_hi]
        if len(u) >= 10 and len(r) >= 10:
            diff = median(u) - median(r)
            _, p = wilcoxon_rank_sum(u, r)
            results.append({"band": name, "n_urban": len(u), "n_rural": len(r),
                           "median_diff_C_decade": round(diff, 4), "p_value": round(p, 4) if p is not None else None})
        else:
            results.append({"band": name, "n_urban": len(u), "n_rural": len(r),
                           "median_diff_C_decade": None, "p_value": None, "note": "insufficient data"})
    return results

def sensitivity_by_min_years(all_trends, classification, annual_data, thresholds=None):
    """Run analysis with different minimum-year thresholds."""
    if thresholds is None:
        thresholds = [30, 40, 50, 60]
    results = []
    for min_yr in thresholds:
        u = [t for sid, t in all_trends if classification.get(sid) == "urban"
             and sid in annual_data and len(annual_data[sid]) >= min_yr]
        r = [t for sid, t in all_trends if classification.get(sid) == "rural"
             and sid in annual_data and len(annual_data[sid]) >= min_yr]
        if len(u) >= 10 and len(r) >= 10:
            diff = median(u) - median(r)
            _, p = wilcoxon_rank_sum(u, r)
            results.append({"min_years": min_yr, "n_urban": len(u), "n_rural": len(r),
                           "median_diff_C_decade": round(diff, 4), "p_value": round(p, 4) if p is not None else None})
        else:
            results.append({"min_years": min_yr, "n_urban": len(u), "n_rural": len(r),
                           "median_diff_C_decade": None, "p_value": None, "note": "insufficient data"})
    return results

###############################################################################
# ID mapping helper
###############################################################################
def map_monthly_to_daily_id(monthly_id, daily_stations):
    """
    GHCN Monthly v4 IDs are 11 chars: 2-char country + 9-char WMO-based.
    GHCN Daily IDs are 11 chars: 2-char country + 1-char network + 8-char ID.
    We try to match by WMO ID or by fuzzy matching.
    For US stations, Monthly uses USW/USC prefix same as Daily.
    """
    # Direct match
    if monthly_id in daily_stations:
        return monthly_id
    # Try matching by prefix (first 5 chars often overlap)
    # This is best-effort since the ID schemes differ
    return None

###############################################################################
# Main analysis
###############################################################################
def main():
    verify_mode = "--verify" in sys.argv

    rng = random.Random(SEED)
    os.makedirs(CACHE_DIR, exist_ok=True)

    total_sections = 9 if not verify_mode else 10
    section = 0

    # ── Section 1: Download data ──────────────────────────────────────────
    section += 1
    print(f"[{section}/{total_sections}] Downloading GHCN data...")
    stations_file = os.path.join(CACHE_DIR, "ghcnd-stations.txt")
    inventory_file = os.path.join(CACHE_DIR, "ghcnd-inventory.txt")
    monthly_tar_file = os.path.join(CACHE_DIR, "ghcnm.tavg.latest.qcf.tar.gz")

    download_with_retry(STATIONS_URL, stations_file, expected_sha256=STATIONS_SHA256)
    download_with_retry(INVENTORY_URL, inventory_file, expected_sha256=INVENTORY_SHA256)
    download_with_retry(MONTHLY_URL, monthly_tar_file, expected_sha256=MONTHLY_SHA256)

    # Extract .dat from tar.gz
    monthly_dat_file = extract_dat_from_tar(monthly_tar_file, CACHE_DIR)

    # Record SHA256 for reproducibility
    sha_stations = sha256_file(stations_file)
    sha_inventory = sha256_file(inventory_file)
    sha_monthly = sha256_file(monthly_tar_file)
    print(f"  Stations SHA256:  {sha_stations}")
    print(f"  Inventory SHA256: {sha_inventory}")
    print(f"  Monthly SHA256:   {sha_monthly}")

    # ── Section 2: Parse station metadata ─────────────────────────────────
    section += 1
    print(f"\n[{section}/{total_sections}] Parsing station metadata...")
    stations = parse_stations(stations_file)
    print(f"  Parsed {len(stations):,} stations from ghcnd-stations.txt")

    inventory = parse_inventory(inventory_file)
    print(f"  Parsed inventory for {len(inventory):,} stations")

    # ── Section 3: Classify stations ──────────────────────────────────────
    section += 1
    print(f"\n[{section}/{total_sections}] Classifying stations as urban/rural...")
    classification = classify_urban_rural(stations)
    n_urban = sum(1 for v in classification.values() if v == "urban")
    n_rural = sum(1 for v in classification.values() if v == "rural")
    print(f"  Classified: {n_urban:,} urban, {n_rural:,} rural, "
          f"{len(stations) - len(classification):,} unclassifiable (excluded)")

    # ── Section 4: Parse monthly temperature data ─────────────────────────
    section += 1
    print(f"\n[{section}/{total_sections}] Parsing GHCN Monthly temperature data...")
    print(f"  (this may take a minute for the full dataset)")
    monthly_temps = parse_monthly_temps(monthly_dat_file)
    print(f"  Parsed temperature data for {len(monthly_temps):,} stations")

    # ── Section 5: Match and compute annual means ─────────────────────────
    section += 1
    print(f"\n[{section}/{total_sections}] Computing annual means and matching stations...")

    # GHCN Monthly v4 station IDs may differ from Daily IDs.
    # For US stations, both use the same 11-char format with USW/USC prefixes.
    # We match by direct ID overlap.
    matched = 0
    annual_data = {}
    for sid in monthly_temps:
        if sid in classification:
            annuals = compute_annual_means(monthly_temps[sid])
            if len(annuals) >= MIN_YEARS:
                annual_data[sid] = annuals
                matched += 1

    print(f"  Matched {matched:,} stations with classification AND >= {MIN_YEARS} years of data")

    # If we have too few matches, also try to match by constructing IDs
    if matched < 100:
        print(f"  Attempting fuzzy matching to increase sample size...")
        # For non-US stations, try prefix matching
        monthly_ids = set(monthly_temps.keys())
        daily_ids = set(classification.keys())
        # Many international stations share the first 2 country chars
        for mid in monthly_ids:
            if mid not in classification:
                # Try finding a daily station with matching country code and last 5+ chars
                for did in daily_ids:
                    if did[:2] == mid[:2] and did[-5:] == mid[-5:] and did not in annual_data:
                        if classification.get(did):
                            annuals = compute_annual_means(monthly_temps[mid])
                            if len(annuals) >= MIN_YEARS:
                                annual_data[did] = annuals
                                matched += 1
                                break
        print(f"  After fuzzy matching: {matched:,} stations")

    n_urban_matched = sum(1 for sid in annual_data if classification.get(sid) == "urban")
    n_rural_matched = sum(1 for sid in annual_data if classification.get(sid) == "rural")
    print(f"  Urban stations with trends: {n_urban_matched:,}")
    print(f"  Rural stations with trends: {n_rural_matched:,}")

    # ── Section 6: Compute Theil-Sen trends ───────────────────────────────
    section += 1
    print(f"\n[{section}/{total_sections}] Computing Theil-Sen trend slopes...")
    t0 = time.time()

    all_trends = []  # (station_id, slope_C_per_decade)
    urban_slopes = []
    rural_slopes = []
    urban_trends_with_id = []
    rural_trends_with_id = []

    for i, (sid, annuals) in enumerate(annual_data.items()):
        slope = theil_sen_slope(annuals)
        if slope is not None:
            all_trends.append((sid, slope))
            if classification.get(sid) == "urban":
                urban_slopes.append(slope)
                urban_trends_with_id.append((sid, slope))
            elif classification.get(sid) == "rural":
                rural_slopes.append(slope)
                rural_trends_with_id.append((sid, slope))
        if (i + 1) % 500 == 0:
            print(f"  Computed {i+1:,}/{len(annual_data):,} trends...")

    elapsed = time.time() - t0
    print(f"  Computed {len(all_trends):,} trend slopes in {elapsed:.1f}s")
    print(f"  Urban: n={len(urban_slopes):,}, median={median(urban_slopes):.4f} C/decade, "
          f"mean={mean(urban_slopes):.4f} C/decade")
    print(f"  Rural: n={len(rural_slopes):,}, median={median(rural_slopes):.4f} C/decade, "
          f"mean={mean(rural_slopes):.4f} C/decade")

    if len(urban_slopes) < 10 or len(rural_slopes) < 10:
        print("ERROR: Insufficient data for statistical analysis.")
        print(f"  Need >= 10 stations in each group, got {len(urban_slopes)} urban, {len(rural_slopes)} rural")
        sys.exit(1)

    # ── Section 7: Statistical tests ──────────────────────────────────────
    section += 1
    print(f"\n[{section}/{total_sections}] Running statistical tests...")

    # Observed difference
    obs_diff = median(urban_slopes) - median(rural_slopes)
    print(f"  Observed median difference (urban - rural): {obs_diff:.4f} C/decade")

    # Wilcoxon rank-sum test
    U, p_wilcoxon = wilcoxon_rank_sum(urban_slopes, rural_slopes)
    print(f"  Wilcoxon rank-sum: U={U:.0f}, p={p_wilcoxon:.6f}")

    # Effect size (Cohen's d)
    d = cohens_d(urban_slopes, rural_slopes)
    print(f"  Cohen's d: {d:.4f}")

    # Permutation test
    print(f"  Running permutation test ({N_PERMUTATIONS:,} permutations)...")
    t0 = time.time()
    p_perm, obs_diff_perm = permutation_test(urban_slopes, rural_slopes, N_PERMUTATIONS, rng)
    elapsed = time.time() - t0
    p_perm_display = f"< {1/N_PERMUTATIONS}" if p_perm == 0 else f"{p_perm:.6f}"
    print(f"  Permutation p-value: {p_perm_display} (computed in {elapsed:.1f}s)")

    # Bootstrap CI
    print(f"  Running bootstrap ({N_BOOTSTRAP:,} resamples)...")
    t0 = time.time()
    ci_lo, ci_hi, boot_diffs = bootstrap_ci(urban_slopes, rural_slopes, N_BOOTSTRAP, rng)
    elapsed = time.time() - t0
    print(f"  95% Bootstrap CI for median difference: [{ci_lo:.4f}, {ci_hi:.4f}] C/decade")
    print(f"  (computed in {elapsed:.1f}s)")

    # ── Section 8: Sensitivity analyses ───────────────────────────────────
    section += 1
    print(f"\n[{section}/{total_sections}] Running sensitivity analyses...")

    lat_results = sensitivity_by_latitude(urban_trends_with_id, rural_trends_with_id, stations, rng)
    print("  By latitude band:")
    for r in lat_results:
        if r.get("median_diff_C_decade") is not None:
            print(f"    {r['band']}: n_u={r['n_urban']}, n_r={r['n_rural']}, "
                  f"diff={r['median_diff_C_decade']:.4f} C/dec, p={r['p_value']}")
        else:
            print(f"    {r['band']}: insufficient data (n_u={r['n_urban']}, n_r={r['n_rural']})")

    elev_results = sensitivity_by_elevation(urban_trends_with_id, rural_trends_with_id, stations, rng)
    print("  By elevation band:")
    for r in elev_results:
        if r.get("median_diff_C_decade") is not None:
            print(f"    {r['band']}: n_u={r['n_urban']}, n_r={r['n_rural']}, "
                  f"diff={r['median_diff_C_decade']:.4f} C/dec, p={r['p_value']}")
        else:
            print(f"    {r['band']}: insufficient data (n_u={r['n_urban']}, n_r={r['n_rural']})")

    min_yr_results = sensitivity_by_min_years(all_trends, classification, annual_data)
    print("  By minimum record length:")
    for r in min_yr_results:
        if r.get("median_diff_C_decade") is not None:
            print(f"    >= {r['min_years']} years: n_u={r['n_urban']}, n_r={r['n_rural']}, "
                  f"diff={r['median_diff_C_decade']:.4f} C/dec, p={r['p_value']}")
        else:
            print(f"    >= {r['min_years']} years: insufficient data")

    # ── Section 9: Write results ──────────────────────────────────────────
    section += 1
    print(f"\n[{section}/{total_sections}] Writing results...")

    results = {
        "metadata": {
            "analysis": "GHCN Urban Heat Island Trend Bias",
            "start_year": START_YEAR,
            "end_year": END_YEAR,
            "min_years_required": MIN_YEARS,
            "n_permutations": N_PERMUTATIONS,
            "n_bootstrap": N_BOOTSTRAP,
            "seed": SEED,
            "sha256_stations": sha_stations,
            "sha256_inventory": sha_inventory,
            "sha256_monthly": sha_monthly,
        },
        "sample_sizes": {
            "total_stations_parsed": len(stations),
            "classified_urban": n_urban,
            "classified_rural": n_rural,
            "matched_with_temp_data": matched,
            "urban_with_trends": len(urban_slopes),
            "rural_with_trends": len(rural_slopes),
        },
        "primary_results": {
            "urban_median_trend_C_per_decade": round(median(urban_slopes), 4),
            "rural_median_trend_C_per_decade": round(median(rural_slopes), 4),
            "urban_mean_trend_C_per_decade": round(mean(urban_slopes), 4),
            "rural_mean_trend_C_per_decade": round(mean(rural_slopes), 4),
            "median_difference_C_per_decade": round(obs_diff, 4),
            "wilcoxon_U": round(U, 1),
            "wilcoxon_p": round(p_wilcoxon, 6),
            "permutation_p": round(p_perm, 6),
            "bootstrap_95ci_lower": round(ci_lo, 4),
            "bootstrap_95ci_upper": round(ci_hi, 4),
            "cohens_d": round(d, 4),
        },
        "sensitivity": {
            "by_latitude": lat_results,
            "by_elevation": elev_results,
            "by_min_years": min_yr_results,
        },
        "interpretation": {
            "significant_at_005": p_perm < ALPHA,
            "direction": "urban warmer" if obs_diff > 0 else "rural warmer" if obs_diff < 0 else "no difference",
            "effect_size_interpretation": (
                "negligible" if abs(d) < 0.2 else
                "small" if abs(d) < 0.5 else
                "medium" if abs(d) < 0.8 else
                "large"
            ),
        },
    }

    with open("results.json", "w") as f:
        json.dump(results, f, indent=2)
    print(f"  Wrote results.json")

    # Write report.md
    report = f"""# GHCN Urban Heat Island Trend Bias Analysis

## Summary

This analysis compares long-term temperature trends (1950-2023) between urban
and rural weather stations in the GHCN network to quantify potential urban heat
island bias in trend estimates.

## Data

- Source: NOAA GHCN Daily station metadata + GHCN Monthly v4 temperature data
- Stations parsed: {len(stations):,}
- Classification: {n_urban:,} urban, {n_rural:,} rural
- Stations with sufficient temperature data: {len(urban_slopes):,} urban, {len(rural_slopes):,} rural
- Analysis window: {START_YEAR}-{END_YEAR}
- Minimum record length: {MIN_YEARS} years

## Results

### Trend Comparison

| Metric | Urban | Rural |
|--------|-------|-------|
| Median trend (C/decade) | {median(urban_slopes):.4f} | {median(rural_slopes):.4f} |
| Mean trend (C/decade) | {mean(urban_slopes):.4f} | {mean(rural_slopes):.4f} |
| Std dev (C/decade) | {stdev(urban_slopes):.4f} | {stdev(rural_slopes):.4f} |
| N stations | {len(urban_slopes):,} | {len(rural_slopes):,} |

### Statistical Tests

| Test | Statistic | p-value |
|------|-----------|---------|
| Wilcoxon rank-sum | U = {U:.0f} | {p_wilcoxon:.6f} |
| Permutation test ({N_PERMUTATIONS:,} shuffles) | — | {p_perm:.6f} |

**Median trend difference (urban − rural): {obs_diff:.4f} C/decade**
**95% Bootstrap CI: [{ci_lo:.4f}, {ci_hi:.4f}] C/decade**
**Cohen's d: {d:.4f} ({results['interpretation']['effect_size_interpretation']})**

### Sensitivity Analysis

#### By Latitude Band

| Band | N urban | N rural | Diff (C/dec) | p-value |
|------|---------|---------|-------------|---------|
"""
    for r in lat_results:
        if r.get("median_diff_C_decade") is not None:
            report += f"| {r['band']} | {r['n_urban']} | {r['n_rural']} | {r['median_diff_C_decade']:.4f} | {r['p_value']} |\n"
        else:
            report += f"| {r['band']} | {r['n_urban']} | {r['n_rural']} | — | — |\n"

    report += f"""
#### By Elevation Band

| Band | N urban | N rural | Diff (C/dec) | p-value |
|------|---------|---------|-------------|---------|
"""
    for r in elev_results:
        if r.get("median_diff_C_decade") is not None:
            report += f"| {r['band']} | {r['n_urban']} | {r['n_rural']} | {r['median_diff_C_decade']:.4f} | {r['p_value']} |\n"
        else:
            report += f"| {r['band']} | {r['n_urban']} | {r['n_rural']} | — | — |\n"

    report += f"""
#### By Minimum Record Length

| Min Years | N urban | N rural | Diff (C/dec) | p-value |
|-----------|---------|---------|-------------|---------|
"""
    for r in min_yr_results:
        if r.get("median_diff_C_decade") is not None:
            report += f"| {r['min_years']} | {r['n_urban']} | {r['n_rural']} | {r['median_diff_C_decade']:.4f} | {r['p_value']} |\n"
        else:
            report += f"| {r['min_years']} | {r['n_urban']} | {r['n_rural']} | — | — |\n"

    report += f"""
## Limitations

1. **Classification heuristic**: Urban/rural classification is based on station ID
   prefixes (USW vs USC/USR) and name keywords, not actual population density or
   land use data. This introduces misclassification noise.
2. **US bias**: The prefix-based classification works best for US stations. International
   stations are only classified by name keywords, reducing global coverage.
3. **Correlation ≠ causation**: Even if urban stations warm faster, this could reflect
   genuine regional climate change rather than (or in addition to) UHI bias.
4. **Temporal confounding**: Station relocations, instrument changes, and time-of-observation
   bias are not controlled for in this analysis.
5. **Survivorship bias**: Stations with >= {MIN_YEARS} years of data are disproportionately
   well-maintained and may not represent the broader network.
6. **Monthly aggregation**: Using monthly means loses information about diurnal UHI
   patterns (UHI is typically stronger at night).
7. **English-language classification bias**: International stations are classified by
   English keywords in station names. Stations with non-English names are excluded,
   biasing coverage toward Anglophone countries.

## Reproducibility

- All data from NOAA GHCN (publicly available)
- SHA256 hashes recorded in results.json
- Random seed: {SEED}
- Python standard library only — zero external dependencies
- Verification mode: `python3 analyze.py --verify`
"""

    with open("report.md", "w") as f:
        f.write(report)
    print(f"  Wrote report.md")

    print(f"\nANALYSIS COMPLETE")

    # ── Verify mode ───────────────────────────────────────────────────────
    if verify_mode:
        section += 1
        print(f"\n[{section}/{total_sections}] Running verification checks...")
        n_pass = 0
        n_fail = 0

        def check(name, condition, detail=""):
            nonlocal n_pass, n_fail
            if condition:
                print(f"  PASS: {name}")
                n_pass += 1
            else:
                print(f"  FAIL: {name} — {detail}")
                n_fail += 1

        # Load results
        with open("results.json", "r") as f:
            res = json.load(f)

        check("results.json exists and is valid JSON", True)
        check("report.md exists", os.path.exists("report.md"))
        check("Urban stations >= 10",
              res["sample_sizes"]["urban_with_trends"] >= 10,
              f"got {res['sample_sizes']['urban_with_trends']}")
        check("Rural stations >= 10",
              res["sample_sizes"]["rural_with_trends"] >= 10,
              f"got {res['sample_sizes']['rural_with_trends']}")
        check("Permutation p-value is valid",
              0 <= res["primary_results"]["permutation_p"] <= 1,
              f"got {res['primary_results']['permutation_p']}")
        check("Bootstrap CI lower < upper",
              res["primary_results"]["bootstrap_95ci_lower"] < res["primary_results"]["bootstrap_95ci_upper"],
              f"got [{res['primary_results']['bootstrap_95ci_lower']}, {res['primary_results']['bootstrap_95ci_upper']}]")
        check("Cohen's d is finite",
              math.isfinite(res["primary_results"]["cohens_d"]),
              f"got {res['primary_results']['cohens_d']}")
        check("Wilcoxon p-value is valid",
              0 <= res["primary_results"]["wilcoxon_p"] <= 1,
              f"got {res['primary_results']['wilcoxon_p']}")
        check("Sensitivity analysis has latitude results",
              len(res["sensitivity"]["by_latitude"]) >= 3,
              f"got {len(res['sensitivity']['by_latitude'])}")
        check("Sensitivity analysis has elevation results",
              len(res["sensitivity"]["by_elevation"]) >= 2,
              f"got {len(res['sensitivity']['by_elevation'])}")
        check("Sensitivity analysis has min-years results",
              len(res["sensitivity"]["by_min_years"]) >= 3,
              f"got {len(res['sensitivity']['by_min_years'])}")
        check("SHA256 hashes recorded",
              all(res["metadata"].get(k) for k in ["sha256_stations", "sha256_inventory", "sha256_monthly"]),
              "missing SHA256 hashes")

        print(f"\n  Verification: {n_pass} passed, {n_fail} failed")
        if n_fail > 0:
            print("  VERIFICATION FAILED")
            sys.exit(1)
        else:
            print("  ALL CHECKS PASSED")

if __name__ == "__main__":
    main()
SCRIPT_EOF
```

**Expected:** Exit code 0. File `ghcn_uhi/analyze.py` exists and is non-empty.

## Step 3: Run analysis

```bash
cd ghcn_uhi && python3 analyze.py
```

**Expected:** Output shows sections `[1/9]` through `[9/9]`, prints statistical results including Wilcoxon p-value, permutation p-value, bootstrap CI, and sensitivity tables. Ends with `ANALYSIS COMPLETE`. Files `results.json` and `report.md` are created.

## Step 4: Verify results

```bash
cd ghcn_uhi && python3 analyze.py --verify
```

**Expected:** Output shows all verification checks with `PASS`. Ends with `ALL CHECKS PASSED`. Exit code 0.

## Success Criteria

- All 4 steps complete with exit code 0
- `results.json` contains valid JSON with all statistical results
- `report.md` contains a readable summary with tables
- Verification mode passes all 11+ checks
- Permutation test runs 2,000 shuffles
- Bootstrap CI computed from 2,000 resamples
- Sensitivity analysis covers latitude bands, elevation bands, and record-length thresholds

## Failure Conditions

- Any step exits with non-zero code
- Network timeout downloading GHCN data (retry up to 3 times)
- Fewer than 10 stations in either urban or rural group after matching
- `results.json` missing required fields
- Verification check fails