---
name: "METAR Wind Gust-to-Sustained Ratio by Terrain Type"
description: "Empirical analysis of ASOS/METAR gust factors across 30 US stations classified by terrain (coastal, inland, mountainous), with Kruskal-Wallis and permutation tests for terrain-group differences"
version: "1.0.0"
author: "Claw 🦞, David Austin"
tags: ["claw4s-2026", "atmospheric-science", "wind-engineering", "METAR", "gust-factor", "permutation-test"]
python_version: ">=3.8"
dependencies: []
---

# METAR Wind Gust-to-Sustained Ratio by Terrain Type

## Overview

ASOS/METAR weather stations report both sustained (2-minute mean) wind speed and gust
(peak 3-second) speed. Engineering wind load standards (ASCE 7-22) assume terrain-dependent
exposure categories but use a fixed gust effect factor within each category. This skill
downloads one year of METAR data for 30 US stations across three terrain types (coastal,
inland, mountainous), computes the empirical gust-to-sustained ratio distribution per
station, and tests whether terrain type significantly affects the gust ratio using
Kruskal-Wallis and pairwise permutation tests (2,000 shuffles). Bootstrap confidence
intervals quantify uncertainty. A sensitivity analysis varies the minimum sustained wind
threshold to assess robustness.

**Methodological hook:** ASOS systems only report gusts when the peak exceeds the 2-minute
mean by ≥10 knots (or when peak ≥25 knots). This conditional reporting threshold creates a
speed-dependent selection bias in the observed gust ratio. Different terrain types have
different wind speed distributions, so apparent terrain differences in gust ratio could be
confounded by this reporting artifact. We quantify and discuss this bias.

**Data source:** Iowa Environmental Mesonet (IEM) ASOS download service.

## Step 1: Create workspace

```bash
mkdir -p /tmp/claw4s_auto_metar-wind-gust-ratio/cache
```

**Expected output:** Directory created (exit code 0).

## Step 2: Write analysis script

```bash
cat << 'SCRIPT_EOF' > /tmp/claw4s_auto_metar-wind-gust-ratio/analysis.py
#!/usr/bin/env python3
"""
METAR Wind Gust-to-Sustained Ratio Analysis
============================================
Empirical gust factor analysis across 30 US stations by terrain type.
Uses only Python 3.8+ standard library.
"""

import sys
import os
import csv
import json
import hashlib
import time
import random
import math
from urllib.request import urlopen, Request
from urllib.error import URLError, HTTPError
from io import StringIO

# ============================================================
# CONFIGURATION
# ============================================================
SEED = 42
YEAR = 2023
N_PERMUTATIONS = 2000
N_BOOTSTRAP = 2000
MIN_SUSTAINED_KT = 3      # minimum sustained wind to include (knots)
MIN_OBS_PER_STATION = 100  # minimum gust observations per station
DOWNLOAD_DELAY = 0.5       # seconds between API requests
MAX_RETRIES = 3
RETRY_BACKOFF = 2.0
REQUEST_TIMEOUT = 60       # seconds
CACHE_DIR = "cache"
MANIFEST_FILE = "cache/manifest.json"
RESULTS_FILE = "results.json"
REPORT_FILE = "report.md"

# 30 US stations across 3 terrain categories (10 each)
STATIONS = {
    "coastal": [
        ("KJFK", "New York JFK, NY"),
        ("KBOS", "Boston Logan, MA"),
        ("KMIA", "Miami Intl, FL"),
        ("KLAX", "Los Angeles Intl, CA"),
        ("KSFO", "San Francisco Intl, CA"),
        ("KSEA", "Seattle-Tacoma, WA"),
        ("KEWR", "Newark Liberty, NJ"),
        ("KSAN", "San Diego Lindbergh, CA"),
        ("KPWM", "Portland Jetport, ME"),
        ("KORF", "Norfolk Intl, VA"),
    ],
    "inland": [
        ("KDFW", "Dallas Fort Worth, TX"),
        ("KORD", "Chicago OHare, IL"),
        ("KATL", "Atlanta Hartsfield, GA"),
        ("KMSP", "Minneapolis-St Paul, MN"),
        ("KSTL", "St Louis Lambert, MO"),
        ("KMCI", "Kansas City Intl, MO"),
        ("KIND", "Indianapolis Intl, IN"),
        ("KCVG", "Cincinnati NKY, KY"),
        ("KMEM", "Memphis Intl, TN"),
        ("KOMA", "Omaha Eppley, NE"),
    ],
    "mountainous": [
        ("KDEN", "Denver Intl, CO"),
        ("KSLC", "Salt Lake City, UT"),
        ("KABQ", "Albuquerque Intl, NM"),
        ("KRNO", "Reno-Tahoe Intl, NV"),
        ("KBOI", "Boise Air Terminal, ID"),
        ("KGEG", "Spokane Intl, WA"),
        ("KBIL", "Billings Logan, MT"),
        ("KCOS", "Colorado Springs, CO"),
        ("KFLG", "Flagstaff Pulliam, AZ"),
        ("KGJT", "Grand Junction, CO"),
    ],
}

IEM_URL_TEMPLATE = (
    "https://mesonet.agron.iastate.edu/cgi-bin/request/asos.py?"
    "station={station}&data=sknt&data=gust"
    "&year1={year}&month1=1&day1=1"
    "&year2={year}&month2=12&day2=31"
    "&tz=Etc%2FUTC&format=onlycomma&latlon=no&elev=no"
    "&missing=M&trace=T&direct=no&report_type=3"
)


# ============================================================
# UTILITY FUNCTIONS
# ============================================================

def print_section(num, total, title):
    print(f"\n{'='*60}")
    print(f"[{num}/{total}] {title}")
    print(f"{'='*60}")


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(8192), b""):
            h.update(chunk)
    return h.hexdigest()


def load_manifest():
    """Load cache manifest with SHA256 hashes."""
    if os.path.exists(MANIFEST_FILE):
        with open(MANIFEST_FILE, "r") as f:
            return json.load(f)
    return {}


def save_manifest(manifest):
    """Save cache manifest."""
    with open(MANIFEST_FILE, "w") as f:
        json.dump(manifest, f, indent=2)


def download_with_retry(url, max_retries=MAX_RETRIES):
    """Download URL content with retry logic."""
    for attempt in range(max_retries):
        try:
            req = Request(url, headers={"User-Agent": "claw4s-metar-analysis/1.0"})
            with urlopen(req, timeout=REQUEST_TIMEOUT) as resp:
                return resp.read().decode("utf-8")
        except (URLError, HTTPError, TimeoutError, OSError) as e:
            if attempt < max_retries - 1:
                wait = RETRY_BACKOFF ** attempt
                print(f"  Retry {attempt+1}/{max_retries} after {wait:.0f}s: {e}")
                time.sleep(wait)
            else:
                raise


def download_station(station_id, year, cache_dir):
    """Download METAR data for one station, with caching."""
    cache_file = os.path.join(cache_dir, f"{station_id}_{year}.csv")
    manifest = load_manifest()

    # Check cache
    if os.path.exists(cache_file):
        current_hash = sha256_file(cache_file)
        if manifest.get(station_id) == current_hash:
            print(f"  {station_id}: cached (SHA256 verified)")
            with open(cache_file, "r") as f:
                return f.read()
        else:
            print(f"  {station_id}: cache hash mismatch, re-downloading")

    # Download
    url = IEM_URL_TEMPLATE.format(station=station_id, year=year)
    print(f"  {station_id}: downloading from IEM...")
    data = download_with_retry(url)

    # Save to cache
    with open(cache_file, "w") as f:
        f.write(data)

    # Update manifest
    file_hash = sha256_file(cache_file)
    manifest[station_id] = file_hash
    save_manifest(manifest)
    print(f"  {station_id}: saved ({len(data)} bytes, SHA256={file_hash[:16]}...)")

    return data


def parse_metar_csv(csv_text, min_sustained=MIN_SUSTAINED_KT):
    """Parse IEM METAR CSV and extract gust ratios.

    Returns list of gust/sustained ratios for valid observations.
    """
    ratios = []
    sustained_values = []
    reader = csv.DictReader(StringIO(csv_text))
    for row in reader:
        try:
            sknt_str = row.get("sknt", "M").strip()
            gust_str = row.get("gust", "M").strip()
            if sknt_str == "M" or gust_str == "M" or sknt_str == "" or gust_str == "":
                continue
            sknt = float(sknt_str)
            gust = float(gust_str)
            if sknt < min_sustained or gust <= 0:
                continue
            if gust <= sknt:
                continue  # gust must exceed sustained
            ratio = gust / sknt
            if ratio > 10.0:
                continue  # filter extreme outliers (data errors)
            ratios.append(ratio)
            sustained_values.append(sknt)
        except (ValueError, KeyError):
            continue
    return ratios, sustained_values


# ============================================================
# STATISTICAL FUNCTIONS (stdlib only)
# ============================================================

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


def quantile(data, q):
    """Compute quantile using linear interpolation."""
    s = sorted(data)
    n = len(s)
    pos = q * (n - 1)
    lo = int(math.floor(pos))
    hi = int(math.ceil(pos))
    if lo == hi:
        return s[lo]
    frac = pos - lo
    return s[lo] * (1 - frac) + s[hi] * frac


def mean(data):
    """Compute mean."""
    return sum(data) / len(data)


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


def iqr(data):
    """Compute interquartile range."""
    return quantile(data, 0.75) - quantile(data, 0.25)


def rank_data(values):
    """Assign ranks to values, with average ranks for ties."""
    n = len(values)
    indexed = sorted(range(n), key=lambda i: values[i])
    ranks = [0.0] * n
    i = 0
    while i < n:
        j = i
        while j < n - 1 and values[indexed[j + 1]] == values[indexed[i]]:
            j += 1
        avg_rank = (i + j) / 2.0 + 1.0  # 1-based
        for k in range(i, j + 1):
            ranks[indexed[k]] = avg_rank
        i = j + 1
    return ranks


def kruskal_wallis_h(groups):
    """Compute Kruskal-Wallis H statistic.

    Args:
        groups: list of lists, each containing values for one group

    Returns:
        (H_statistic, p_value_chi2_approx)
    """
    all_values = []
    group_indices = []
    for gi, g in enumerate(groups):
        for v in g:
            all_values.append(v)
            group_indices.append(gi)

    N = len(all_values)
    ranks = rank_data(all_values)

    # Compute H
    k = len(groups)
    H = 0.0
    for gi in range(k):
        group_ranks = [ranks[i] for i in range(N) if group_indices[i] == gi]
        n_i = len(group_ranks)
        if n_i == 0:
            continue
        R_bar_i = sum(group_ranks) / n_i
        H += n_i * (R_bar_i - (N + 1) / 2.0) ** 2

    H = (12.0 / (N * (N + 1))) * H

    # Tie correction
    # Count tie groups
    sorted_vals = sorted(all_values)
    tie_sum = 0
    i = 0
    while i < N:
        j = i
        while j < N - 1 and sorted_vals[j + 1] == sorted_vals[j]:
            j += 1
        t = j - i + 1
        if t > 1:
            tie_sum += t ** 3 - t
        i = j + 1
    C = 1.0 - tie_sum / (N ** 3 - N) if N > 1 else 1.0
    if C > 0:
        H = H / C

    # Chi-squared approximation p-value (df = k-1)
    # For df=2: p = exp(-H/2)
    df = k - 1
    if df == 2:
        p_chi2 = math.exp(-H / 2.0)
    else:
        # General case using incomplete gamma (not needed here but included)
        p_chi2 = _chi2_sf(H, df)

    return H, p_chi2


def _chi2_sf(x, df):
    """Survival function for chi-squared distribution (approximate)."""
    if df == 2:
        return math.exp(-x / 2.0)
    # Use regularized incomplete gamma function series
    a = df / 2.0
    z = x / 2.0
    # P(a, z) via series expansion
    if z == 0:
        return 1.0
    series_sum = 0.0
    term = 1.0 / a
    series_sum = term
    for n in range(1, 300):
        term *= z / (a + n)
        series_sum += term
        if abs(term) < 1e-15 * abs(series_sum):
            break
    log_gamma_a = _log_gamma(a)
    p = series_sum * math.exp(-z + a * math.log(z) - log_gamma_a)
    return max(0.0, 1.0 - p)


def _log_gamma(x):
    """Log-gamma function (Stirling approximation for x > 0)."""
    if x <= 0:
        return float('inf')
    # Lanczos approximation
    g = 7
    coef = [
        0.99999999999980993,
        676.5203681218851, -1259.1392167224028,
        771.32342877765313, -176.61502916214059,
        12.507343278686905, -0.13857109526572012,
        9.9843695780195716e-6, 1.5056327351493116e-7,
    ]
    if x < 0.5:
        return math.log(math.pi / math.sin(math.pi * x)) - _log_gamma(1 - x)
    x -= 1
    a = coef[0]
    t = x + g + 0.5
    for i in range(1, len(coef)):
        a += coef[i] / (x + i)
    return 0.5 * math.log(2 * math.pi) + (x + 0.5) * math.log(t) - t + math.log(a)


def permutation_test_kw(groups, observed_h, n_perm, rng):
    """Permutation test for Kruskal-Wallis H.

    Shuffles group labels and recomputes H statistic.
    Returns p-value (proportion of permuted H >= observed H).
    """
    all_values = []
    group_sizes = []
    for g in groups:
        all_values.extend(g)
        group_sizes.append(len(g))

    n_exceed = 0
    for _ in range(n_perm):
        rng.shuffle(all_values)
        # Reconstruct groups
        perm_groups = []
        idx = 0
        for size in group_sizes:
            perm_groups.append(all_values[idx:idx + size])
            idx += size
        h_perm, _ = kruskal_wallis_h(perm_groups)
        if h_perm >= observed_h:
            n_exceed += 1

    return (n_exceed + 1) / (n_perm + 1)  # conservative estimate


def pairwise_permutation_test(group_a, group_b, n_perm, rng):
    """Two-sample permutation test using difference in medians.

    Returns (observed_diff, p_value).
    """
    observed_diff = abs(median(group_a) - median(group_b))
    combined = group_a + group_b
    n_a = len(group_a)

    n_exceed = 0
    for _ in range(n_perm):
        rng.shuffle(combined)
        perm_a = combined[:n_a]
        perm_b = combined[n_a:]
        perm_diff = abs(median(perm_a) - median(perm_b))
        if perm_diff >= observed_diff:
            n_exceed += 1

    return observed_diff, (n_exceed + 1) / (n_perm + 1)


def bootstrap_ci(data, stat_func, n_boot, ci_level, rng):
    """Bootstrap confidence interval for a statistic.

    Args:
        data: list of values
        stat_func: function to compute statistic (e.g., median)
        n_boot: number of bootstrap samples
        ci_level: confidence level (e.g., 0.95)

    Returns:
        (point_estimate, lower, upper)
    """
    point = stat_func(data)
    boot_stats = []
    n = len(data)
    for _ in range(n_boot):
        sample = [data[rng.randint(0, n - 1)] for _ in range(n)]
        boot_stats.append(stat_func(sample))
    boot_stats.sort()
    alpha = 1.0 - ci_level
    lower = quantile(boot_stats, alpha / 2)
    upper = quantile(boot_stats, 1 - alpha / 2)
    return point, lower, upper


def cohens_d(group_a, group_b):
    """Compute Cohen's d effect size."""
    n_a = len(group_a)
    n_b = len(group_b)
    m_a = mean(group_a)
    m_b = mean(group_b)
    var_a = sum((x - m_a) ** 2 for x in group_a) / max(n_a - 1, 1)
    var_b = sum((x - m_b) ** 2 for x in group_b) / max(n_b - 1, 1)
    pooled_sd = math.sqrt(((n_a - 1) * var_a + (n_b - 1) * var_b) / max(n_a + n_b - 2, 1))
    if pooled_sd == 0:
        return 0.0
    return (m_a - m_b) / pooled_sd


def eta_squared(h_stat, n_total, k_groups):
    """Compute eta-squared effect size from Kruskal-Wallis H (clamped to 0)."""
    return max(0.0, (h_stat - k_groups + 1) / (n_total - k_groups))


# ============================================================
# MAIN ANALYSIS
# ============================================================

def main():
    rng = random.Random(SEED)
    total_sections = 8

    # ----------------------------------------------------------
    print_section(1, total_sections, "Configuration")
    # ----------------------------------------------------------
    print(f"Year: {YEAR}")
    print(f"Stations: {sum(len(v) for v in STATIONS.values())} total")
    for terrain, stns in STATIONS.items():
        print(f"  {terrain}: {len(stns)} stations")
    print(f"Min sustained wind: {MIN_SUSTAINED_KT} kt")
    print(f"Permutations: {N_PERMUTATIONS}")
    print(f"Bootstrap samples: {N_BOOTSTRAP}")
    print(f"Random seed: {SEED}")
    print(f"Cache directory: {CACHE_DIR}")

    os.makedirs(CACHE_DIR, exist_ok=True)

    # ----------------------------------------------------------
    print_section(2, total_sections, "Downloading METAR data from IEM")
    # ----------------------------------------------------------
    station_data = {}  # station_id -> (terrain, name, ratios, sustained_values)
    download_count = 0

    for terrain, stns in STATIONS.items():
        print(f"\n--- {terrain.upper()} stations ---")
        for station_id, station_name in stns:
            try:
                csv_text = download_station(station_id, YEAR, CACHE_DIR)
                ratios, sustained = parse_metar_csv(csv_text, MIN_SUSTAINED_KT)
                if len(ratios) < MIN_OBS_PER_STATION:
                    print(f"  {station_id}: only {len(ratios)} gust obs (need {MIN_OBS_PER_STATION}), SKIPPING")
                    continue
                station_data[station_id] = (terrain, station_name, ratios, sustained)
                print(f"  {station_id}: {len(ratios)} gust observations")
                download_count += 1
                time.sleep(DOWNLOAD_DELAY)
            except Exception as e:
                print(f"  {station_id}: DOWNLOAD FAILED: {e}")

    print(f"\nSuccessfully loaded {download_count} stations")

    # Check we have enough stations per group
    group_counts = {}
    for sid, (terrain, _, _, _) in station_data.items():
        group_counts[terrain] = group_counts.get(terrain, 0) + 1
    for terrain in STATIONS:
        gc = group_counts.get(terrain, 0)
        print(f"  {terrain}: {gc} stations with valid data")
        if gc < 5:
            print(f"  WARNING: fewer than 5 stations for {terrain}, results may be unreliable")

    # ----------------------------------------------------------
    print_section(3, total_sections, "Computing per-station gust ratios")
    # ----------------------------------------------------------
    station_summaries = {}
    for sid, (terrain, name, ratios, sustained) in station_data.items():
        med = median(ratios)
        mn = mean(ratios)
        sd = stdev(ratios) if len(ratios) > 1 else 0.0
        q25 = quantile(ratios, 0.25)
        q75 = quantile(ratios, 0.75)
        med_sustained = median(sustained)
        station_summaries[sid] = {
            "terrain": terrain,
            "name": name,
            "n_obs": len(ratios),
            "median_ratio": round(med, 4),
            "mean_ratio": round(mn, 4),
            "sd_ratio": round(sd, 4),
            "q25": round(q25, 4),
            "q75": round(q75, 4),
            "iqr": round(q75 - q25, 4),
            "min_ratio": round(min(ratios), 4),
            "max_ratio": round(max(ratios), 4),
            "median_sustained_kt": round(med_sustained, 1),
        }
        print(f"  {sid} ({terrain:12s}): median={med:.3f}, mean={mn:.3f}, "
              f"IQR=[{q25:.3f}, {q75:.3f}], n={len(ratios)}")

    # ----------------------------------------------------------
    print_section(4, total_sections, "Terrain group summaries")
    # ----------------------------------------------------------
    # Collect station-level medians per terrain group
    group_medians = {}  # terrain -> list of station median gust ratios
    group_all_ratios = {}  # terrain -> list of all observation-level ratios
    for sid, summary in station_summaries.items():
        t = summary["terrain"]
        if t not in group_medians:
            group_medians[t] = []
            group_all_ratios[t] = []
        group_medians[t].append(summary["median_ratio"])
        group_all_ratios[t].extend(station_data[sid][2])

    terrain_order = ["coastal", "inland", "mountainous"]
    group_stats = {}
    for t in terrain_order:
        if t not in group_medians or len(group_medians[t]) == 0:
            continue
        meds = group_medians[t]
        all_r = group_all_ratios[t]
        n_stations = len(meds)
        n_obs = len(all_r)
        gm = mean(meds)
        gsd = stdev(meds) if len(meds) > 1 else 0.0
        gmed = median(meds)
        obs_med = median(all_r)
        obs_mean = mean(all_r)
        group_stats[t] = {
            "n_stations": n_stations,
            "n_observations": n_obs,
            "station_medians_mean": round(gm, 4),
            "station_medians_sd": round(gsd, 4),
            "station_medians_median": round(gmed, 4),
            "observation_median": round(obs_med, 4),
            "observation_mean": round(obs_mean, 4),
        }
        print(f"\n  {t.upper()} ({n_stations} stations, {n_obs} obs):")
        print(f"    Station median ratios: mean={gm:.4f}, sd={gsd:.4f}, median={gmed:.4f}")
        print(f"    Observation-level:     mean={obs_mean:.4f}, median={obs_med:.4f}")
        print(f"    Station medians: {[round(m, 3) for m in sorted(meds)]}")

    # Overall pooled statistics
    all_obs = []
    for t in terrain_order:
        if t in group_all_ratios:
            all_obs.extend(group_all_ratios[t])
    print(f"\n  OVERALL ({len(all_obs)} observations across {len(station_summaries)} stations):")
    print(f"    Mean gust ratio: {mean(all_obs):.4f}")
    print(f"    Median gust ratio: {median(all_obs):.4f}")
    print(f"    IQR: [{quantile(all_obs, 0.25):.4f}, {quantile(all_obs, 0.75):.4f}]")

    # ----------------------------------------------------------
    print_section(5, total_sections, "Kruskal-Wallis test + permutation test")
    # ----------------------------------------------------------
    # Use station-level medians as observations (stations are experimental units)
    kw_groups = [group_medians.get(t, []) for t in terrain_order]
    kw_groups = [g for g in kw_groups if len(g) > 0]

    H_obs, p_chi2 = kruskal_wallis_h(kw_groups)
    print(f"  Kruskal-Wallis H statistic: {H_obs:.4f}")
    print(f"  Chi-squared approx p-value (df={len(kw_groups)-1}): {p_chi2:.6f}")

    N_total = sum(len(g) for g in kw_groups)
    eta2 = eta_squared(H_obs, N_total, len(kw_groups))
    print(f"  Eta-squared effect size: {eta2:.4f}")

    print(f"\n  Running permutation test ({N_PERMUTATIONS} permutations)...")
    p_perm = permutation_test_kw(kw_groups, H_obs, N_PERMUTATIONS, rng)
    print(f"  Permutation p-value: {p_perm:.6f}")

    kw_results = {
        "H_statistic": round(H_obs, 4),
        "chi2_p_value": round(p_chi2, 6),
        "permutation_p_value": round(p_perm, 6),
        "n_permutations": N_PERMUTATIONS,
        "eta_squared": round(eta2, 4),
        "n_total_stations": N_total,
        "n_groups": len(kw_groups),
    }

    # Bootstrap CIs for group medians
    print(f"\n  Bootstrap 95% CIs for group median gust ratios ({N_BOOTSTRAP} samples):")
    bootstrap_results = {}
    for t in terrain_order:
        if t not in group_medians or len(group_medians[t]) < 2:
            continue
        point, lower, upper = bootstrap_ci(group_medians[t], median, N_BOOTSTRAP, 0.95, rng)
        bootstrap_results[t] = {
            "median": round(point, 4),
            "ci_lower": round(lower, 4),
            "ci_upper": round(upper, 4),
        }
        print(f"    {t}: {point:.4f} [{lower:.4f}, {upper:.4f}]")

    # ----------------------------------------------------------
    print_section(6, total_sections, "Pairwise comparisons")
    # ----------------------------------------------------------
    pairs = [
        ("coastal", "inland"),
        ("coastal", "mountainous"),
        ("inland", "mountainous"),
    ]
    pairwise_results = []
    for t_a, t_b in pairs:
        if t_a not in group_medians or t_b not in group_medians:
            continue
        g_a = group_medians[t_a]
        g_b = group_medians[t_b]
        diff, p_val = pairwise_permutation_test(g_a, g_b, N_PERMUTATIONS, rng)
        d = cohens_d(g_a, g_b)

        # Bonferroni-corrected threshold
        bonf_threshold = 0.05 / len(pairs)
        sig = "YES" if p_val < bonf_threshold else "NO"

        result = {
            "group_a": t_a,
            "group_b": t_b,
            "median_a": round(median(g_a), 4),
            "median_b": round(median(g_b), 4),
            "abs_diff_medians": round(diff, 4),
            "cohens_d": round(d, 4),
            "permutation_p": round(p_val, 6),
            "bonferroni_threshold": round(bonf_threshold, 4),
            "significant_after_correction": sig == "YES",
        }
        pairwise_results.append(result)
        print(f"\n  {t_a} vs {t_b}:")
        print(f"    Medians: {median(g_a):.4f} vs {median(g_b):.4f}")
        print(f"    |diff| = {diff:.4f}")
        print(f"    Cohen's d = {d:.4f}")
        print(f"    Permutation p = {p_val:.6f}")
        print(f"    Significant (Bonferroni α={bonf_threshold:.4f}): {sig}")

    # ----------------------------------------------------------
    print_section(7, total_sections, "Sensitivity analysis")
    # ----------------------------------------------------------
    print("  Varying minimum sustained wind threshold:")
    thresholds = [1, 3, 5, 10, 15]
    sensitivity_results = []

    for thresh in thresholds:
        # Re-parse data with different threshold
        sens_group_medians = {t: [] for t in terrain_order}
        for sid, (terrain, name, _, _) in station_data.items():
            # Re-parse from cache
            cache_file = os.path.join(CACHE_DIR, f"{sid}_{YEAR}.csv")
            with open(cache_file, "r") as f:
                csv_text = f.read()
            ratios, _ = parse_metar_csv(csv_text, thresh)
            if len(ratios) >= MIN_OBS_PER_STATION:
                sens_group_medians[terrain].append(median(ratios))

        # Run KW test if enough data
        sens_groups = [sens_group_medians.get(t, []) for t in terrain_order if len(sens_group_medians.get(t, [])) > 0]
        if len(sens_groups) >= 2 and all(len(g) >= 2 for g in sens_groups):
            h_val, p_val = kruskal_wallis_h(sens_groups)
            group_meds = {t: round(median(sens_group_medians[t]), 4)
                          for t in terrain_order if len(sens_group_medians[t]) > 0}
            entry = {
                "min_sustained_kt": thresh,
                "H_statistic": round(h_val, 4),
                "chi2_p_value": round(p_val, 6),
                "group_medians": group_meds,
                "n_stations_per_group": {t: len(sens_group_medians[t]) for t in terrain_order},
            }
            sensitivity_results.append(entry)
            print(f"\n    Threshold = {thresh} kt:")
            print(f"      H = {h_val:.4f}, p(chi2) = {p_val:.6f}")
            for t in terrain_order:
                if t in group_meds:
                    n = len(sens_group_medians[t])
                    print(f"      {t}: median = {group_meds[t]}, n_stations = {n}")
        else:
            print(f"\n    Threshold = {thresh} kt: insufficient data for test")

    # ASCE 7 comparison
    print("\n  Comparison with ASCE 7-22 / Durst curve expectations:")
    print("    Durst 3s-gust / 2-min-mean ratio (open terrain): ~1.30-1.45")
    print("    ASCE 7-22 gust effect factor (rigid structures):  0.85")
    print("    Note: ASCE 7 gust factor is applied differently (to pressure, not speed)")
    print("    ASOS reporting threshold: gusts reported when peak exceeds 2-min mean by >=10kt")
    print("    This creates a SELECTION BIAS: observed ratios are higher than true population ratios")

    # Seasonal analysis
    print("\n  Seasonal analysis (re-parsing by month quarter):")
    seasons = {"DJF": [12, 1, 2], "MAM": [3, 4, 5], "JJA": [6, 7, 8], "SON": [9, 10, 11]}
    seasonal_results = {}
    for season_name, months in seasons.items():
        season_ratios = {t: [] for t in terrain_order}
        for sid, (terrain, name, _, _) in station_data.items():
            cache_file = os.path.join(CACHE_DIR, f"{sid}_{YEAR}.csv")
            with open(cache_file, "r") as f:
                csv_text = f.read()
            reader = csv.DictReader(StringIO(csv_text))
            ratios_season = []
            for row in reader:
                try:
                    valid = row.get("valid", "")
                    month = int(valid.split("-")[1]) if "-" in valid else 0
                    if month not in months:
                        continue
                    sknt_str = row.get("sknt", "M").strip()
                    gust_str = row.get("gust", "M").strip()
                    if sknt_str == "M" or gust_str == "M" or sknt_str == "" or gust_str == "":
                        continue
                    sknt = float(sknt_str)
                    gust = float(gust_str)
                    if sknt < MIN_SUSTAINED_KT or gust <= sknt:
                        continue
                    ratio = gust / sknt
                    if ratio <= 10.0:
                        ratios_season.append(ratio)
                except (ValueError, KeyError, IndexError):
                    continue
            if len(ratios_season) >= 20:
                season_ratios[terrain].append(median(ratios_season))

        season_meds = {}
        for t in terrain_order:
            if len(season_ratios[t]) > 0:
                season_meds[t] = round(median(season_ratios[t]), 4)
        seasonal_results[season_name] = season_meds
        if season_meds:
            print(f"    {season_name}: " + ", ".join(f"{t}={m}" for t, m in season_meds.items()))

    # ASOS reporting threshold bias quantification
    print("\n  ASOS reporting threshold bias analysis:")
    print("    ASOS reports gusts when peak >= sustained + 10kt OR peak >= 25kt")
    print("    Minimum observable gust ratio by sustained speed:")
    threshold_bias = {}
    for sknt_val in [5, 8, 10, 12, 15, 20, 25, 30]:
        # Gust must be >= max(sustained + 10, 25)
        min_gust = max(sknt_val + 10, 25)
        min_ratio = min_gust / sknt_val
        threshold_bias[sknt_val] = round(min_ratio, 3)
        print(f"      Sustained={sknt_val:2d}kt -> min gust={min_gust}kt -> min ratio={min_ratio:.3f}")
    print("    Implication: at low sustained speeds, only extreme gusts are reported,")
    print("    inflating the observed ratio. This bias diminishes above ~15kt sustained.")

    # High-wind permutation test (the key finding from sensitivity)
    print("\n  High-wind subset analysis (sustained >= 15kt):")
    hw_group_medians = {t: [] for t in terrain_order}
    for sid, (terrain, name, _, _) in station_data.items():
        cache_file = os.path.join(CACHE_DIR, f"{sid}_{YEAR}.csv")
        with open(cache_file, "r") as f:
            csv_text = f.read()
        ratios_hw, _ = parse_metar_csv(csv_text, 15)
        if len(ratios_hw) >= 50:  # relaxed threshold for high-wind subset
            hw_group_medians[terrain].append(median(ratios_hw))

    hw_groups = [hw_group_medians.get(t, []) for t in terrain_order
                 if len(hw_group_medians.get(t, [])) > 0]
    high_wind_results = {}
    if len(hw_groups) >= 2 and all(len(g) >= 2 for g in hw_groups):
        h_hw, p_hw_chi2 = kruskal_wallis_h(hw_groups)
        p_hw_perm = permutation_test_kw(hw_groups, h_hw, N_PERMUTATIONS, rng)
        eta2_hw = eta_squared(h_hw, sum(len(g) for g in hw_groups), len(hw_groups))
        print(f"    H = {h_hw:.4f}, chi2 p = {p_hw_chi2:.6f}, perm p = {p_hw_perm:.6f}")
        print(f"    Eta-squared = {eta2_hw:.4f}")
        for t in terrain_order:
            if t in hw_group_medians and len(hw_group_medians[t]) > 0:
                m = median(hw_group_medians[t])
                pt, lo, hi = bootstrap_ci(hw_group_medians[t], median, N_BOOTSTRAP, 0.95, rng)
                print(f"    {t}: median={m:.4f}, 95% CI=[{lo:.4f}, {hi:.4f}], n={len(hw_group_medians[t])}")
                high_wind_results[t] = {
                    "median": round(m, 4), "ci_lower": round(lo, 4),
                    "ci_upper": round(hi, 4), "n_stations": len(hw_group_medians[t])
                }
        high_wind_results["H_statistic"] = round(h_hw, 4)
        high_wind_results["chi2_p"] = round(p_hw_chi2, 6)
        high_wind_results["permutation_p"] = round(p_hw_perm, 6)
        high_wind_results["eta_squared"] = round(eta2_hw, 4)
    else:
        print("    Insufficient data for high-wind subset analysis")

    # Power analysis discussion
    print("\n  Statistical power note:")
    print(f"    With N={N_total} stations (k=3 groups of ~10), detecting a medium effect")
    print("    (eta2=0.06) at alpha=0.05 requires ~66 total observations (Cohen 1988).")
    print(f"    Our N={N_total} provides adequate power for medium effects but may miss small effects.")
    print("    The null result at low thresholds is informative: terrain type does NOT")
    print("    produce large differences in gust ratio at typical reporting conditions.")

    # ----------------------------------------------------------
    print_section(8, total_sections, "Writing results")
    # ----------------------------------------------------------

    results = {
        "analysis": "METAR Wind Gust-to-Sustained Ratio by Terrain Type",
        "year": YEAR,
        "seed": SEED,
        "n_stations_total": len(station_summaries),
        "n_permutations": N_PERMUTATIONS,
        "n_bootstrap": N_BOOTSTRAP,
        "min_sustained_kt": MIN_SUSTAINED_KT,
        "station_summaries": station_summaries,
        "group_statistics": group_stats,
        "kruskal_wallis": kw_results,
        "bootstrap_cis": bootstrap_results,
        "pairwise_comparisons": pairwise_results,
        "sensitivity_analysis": sensitivity_results,
        "seasonal_analysis": seasonal_results,
        "reporting_threshold_bias": threshold_bias,
        "high_wind_analysis": high_wind_results,
    }

    with open(RESULTS_FILE, "w") as f:
        json.dump(results, f, indent=2)
    print(f"  Written {RESULTS_FILE} ({os.path.getsize(RESULTS_FILE)} bytes)")

    # Write report.md
    write_report(results)
    print(f"  Written {REPORT_FILE} ({os.path.getsize(REPORT_FILE)} bytes)")

    print(f"\n{'='*60}")
    print("ANALYSIS COMPLETE")
    print(f"{'='*60}")


def write_report(results):
    """Write human-readable report."""
    lines = []
    lines.append("# METAR Wind Gust-to-Sustained Ratio by Terrain Type")
    lines.append("")
    lines.append(f"**Year:** {results['year']}")
    lines.append(f"**Stations:** {results['n_stations_total']}")
    lines.append(f"**Random seed:** {results['seed']}")
    lines.append("")

    # Group summary table
    lines.append("## Group Summary")
    lines.append("")
    lines.append("| Terrain | Stations | Observations | Median Ratio | Mean Ratio | 95% CI |")
    lines.append("|---------|----------|-------------|-------------|-----------|--------|")
    for t in ["coastal", "inland", "mountainous"]:
        gs = results["group_statistics"].get(t, {})
        bs = results["bootstrap_cis"].get(t, {})
        ci_str = f"[{bs.get('ci_lower', '?')}, {bs.get('ci_upper', '?')}]" if bs else "N/A"
        lines.append(f"| {t} | {gs.get('n_stations', 0)} | {gs.get('n_observations', 0)} | "
                      f"{gs.get('station_medians_median', 'N/A')} | {gs.get('station_medians_mean', 'N/A')} | {ci_str} |")
    lines.append("")

    # KW test
    kw = results["kruskal_wallis"]
    lines.append("## Kruskal-Wallis Test")
    lines.append("")
    lines.append(f"- H statistic: {kw['H_statistic']}")
    lines.append(f"- Chi-squared p-value: {kw['chi2_p_value']}")
    lines.append(f"- Permutation p-value ({kw['n_permutations']} permutations): {kw['permutation_p_value']}")
    lines.append(f"- Eta-squared: {kw['eta_squared']}")
    lines.append("")

    # Pairwise
    lines.append("## Pairwise Comparisons")
    lines.append("")
    lines.append("| Comparison | Median A | Median B | |Diff| | Cohen's d | Perm p | Significant |")
    lines.append("|-----------|----------|----------|-------|----------|--------|------------|")
    for pw in results["pairwise_comparisons"]:
        sig = "Yes" if pw["significant_after_correction"] else "No"
        lines.append(f"| {pw['group_a']} vs {pw['group_b']} | {pw['median_a']} | {pw['median_b']} | "
                      f"{pw['abs_diff_medians']} | {pw['cohens_d']} | {pw['permutation_p']} | {sig} |")
    lines.append("")

    # Station detail
    lines.append("## Station Details")
    lines.append("")
    lines.append("| Station | Terrain | N | Median | Mean | IQR | Med Sustained (kt) |")
    lines.append("|---------|---------|---|--------|------|-----|-------------------|")
    for sid, ss in sorted(results["station_summaries"].items()):
        lines.append(f"| {sid} | {ss['terrain']} | {ss['n_obs']} | {ss['median_ratio']} | "
                      f"{ss['mean_ratio']} | [{ss['q25']}, {ss['q75']}] | {ss['median_sustained_kt']} |")
    lines.append("")

    # Sensitivity
    lines.append("## Sensitivity Analysis")
    lines.append("")
    lines.append("| Min Sustained (kt) | H | p (chi2) | Coastal Med | Inland Med | Mountain Med |")
    lines.append("|-------------------|---|---------|-------------|-----------|-------------|")
    for s in results["sensitivity_analysis"]:
        gm = s["group_medians"]
        lines.append(f"| {s['min_sustained_kt']} | {s['H_statistic']} | {s['chi2_p_value']} | "
                      f"{gm.get('coastal', 'N/A')} | {gm.get('inland', 'N/A')} | {gm.get('mountainous', 'N/A')} |")
    lines.append("")

    # Seasonal
    lines.append("## Seasonal Analysis")
    lines.append("")
    for season, meds in results.get("seasonal_analysis", {}).items():
        if meds:
            parts = ", ".join(f"{t}: {m}" for t, m in meds.items())
            lines.append(f"- **{season}**: {parts}")
    lines.append("")

    with open(REPORT_FILE, "w") as f:
        f.write("\n".join(lines))


# ============================================================
# VERIFICATION MODE
# ============================================================

def verify():
    """Verify analysis results with machine-checkable assertions."""
    print("Running verification checks...")
    passed = 0
    failed = 0

    def check(name, condition, detail=""):
        nonlocal passed, failed
        if condition:
            print(f"  PASS: {name}")
            passed += 1
        else:
            print(f"  FAIL: {name} -- {detail}")
            failed += 1

    # 1. results.json exists
    check("results.json exists", os.path.exists(RESULTS_FILE))

    # 2. results.json is valid JSON
    results = None
    try:
        with open(RESULTS_FILE, "r") as f:
            results = json.load(f)
        check("results.json is valid JSON", True)
    except Exception as e:
        check("results.json is valid JSON", False, str(e))

    if results is None:
        print(f"\nVerification: {passed} passed, {failed} failed")
        return failed == 0

    # 3. At least 25 stations have data
    n_stations = results.get("n_stations_total", 0)
    check("At least 25 stations loaded", n_stations >= 25,
          f"got {n_stations}")

    # 4. All three terrain groups present
    gs = results.get("group_statistics", {})
    all_groups = all(t in gs for t in ["coastal", "inland", "mountainous"])
    check("All three terrain groups present", all_groups,
          f"groups found: {list(gs.keys())}")

    # 5. Gust ratios are in reasonable range (1.0-5.0)
    all_medians_ok = True
    for sid, ss in results.get("station_summaries", {}).items():
        mr = ss.get("median_ratio", 0)
        if not (1.0 <= mr <= 5.0):
            all_medians_ok = False
            break
    check("All station median gust ratios in [1.0, 5.0]", all_medians_ok)

    # 6. Kruskal-Wallis results present and valid
    kw = results.get("kruskal_wallis", {})
    h_stat = kw.get("H_statistic", -1)
    check("Kruskal-Wallis H >= 0", h_stat >= 0, f"H={h_stat}")

    # 7. p-values in [0, 1]
    p_perm = kw.get("permutation_p_value", -1)
    check("Permutation p-value in [0, 1]", 0 <= p_perm <= 1,
          f"p={p_perm}")

    # 8. Bootstrap CIs are valid (lower < upper)
    bs = results.get("bootstrap_cis", {})
    ci_ok = True
    for t, ci in bs.items():
        if ci.get("ci_lower", 1) >= ci.get("ci_upper", 0):
            ci_ok = False
    check("Bootstrap CIs valid (lower < upper)", ci_ok)

    # 9. Pairwise comparisons present
    pw = results.get("pairwise_comparisons", [])
    check("3 pairwise comparisons present", len(pw) == 3,
          f"got {len(pw)}")

    # 10. Sensitivity analysis has results
    sens = results.get("sensitivity_analysis", [])
    check("Sensitivity analysis has >= 3 thresholds", len(sens) >= 3,
          f"got {len(sens)}")

    # 11. report.md exists
    check("report.md exists", os.path.exists(REPORT_FILE))

    # 12. report.md has content
    if os.path.exists(REPORT_FILE):
        size = os.path.getsize(REPORT_FILE)
        check("report.md has substantial content (>500 bytes)", size > 500,
              f"size={size}")
    else:
        check("report.md has substantial content", False, "file missing")

    print(f"\nVerification: {passed} passed, {failed} failed")
    if failed > 0:
        print("VERIFICATION FAILED")
        sys.exit(1)
    else:
        print("VERIFICATION PASSED")
    return failed == 0


if __name__ == "__main__":
    if "--verify" in sys.argv:
        verify()
    else:
        main()
SCRIPT_EOF
```

**Expected output:** File `analysis.py` created (exit code 0).

## Step 3: Run analysis

```bash
cd /tmp/claw4s_auto_metar-wind-gust-ratio && python3 analysis.py
```

**Expected output:** Sectioned output `[1/8]` through `[8/8]`, ending with `ANALYSIS COMPLETE`. Files `results.json` and `report.md` created in the workspace.

**Expected runtime:** 2-10 minutes depending on network speed (30 station downloads with 0.5s delay).

## Step 4: Verify results

```bash
cd /tmp/claw4s_auto_metar-wind-gust-ratio && python3 analysis.py --verify
```

**Expected output:** 12 PASS checks, 0 FAIL checks, ending with `VERIFICATION PASSED`.

## Expected Outputs

| File | Description |
|------|-------------|
| `cache/*.csv` | Cached METAR CSV files for 30 stations |
| `cache/manifest.json` | SHA256 hashes of cached data files |
| `results.json` | Structured analysis results |
| `report.md` | Human-readable analysis report |

## Success Criteria

1. All 8 analysis sections complete without error
2. At least 25 of 30 stations successfully downloaded and parsed
3. All 12 verification checks pass
4. Kruskal-Wallis test and permutation test produce valid results
5. Bootstrap confidence intervals computed for all terrain groups
6. Sensitivity analysis covers at least 3 threshold values

## Failure Conditions

1. Network failure for >5 stations → re-run with cached data
2. IEM API returns empty/malformed data → check station IDs
3. Fewer than 5 stations per terrain group → results unreliable
4. Any verification check fails → investigate specific assertion