#!/usr/bin/env python3 """Radial scaling of ICME magnetic-obstacle mean field strength with heliocentric distance. Pre-registered analysis of the HELIO4CAST Interplanetary Coronal Mass Ejection Catalog v2.3 (ICMECAT). Estimates the power-law exponent alpha in B ~ r^-alpha for the magnetic-obstacle mean field strength (mo_bmean) versus spacecraft heliocentric distance (mo_sc_heliodistance), using nonparametric and robust methods, and guards against the spacecraft / heliodistance sampling confound. Re-running this script regenerates every figure, table, and statistic deterministically. Data are fetched from the canonical figshare URL and md5-verified; if the network is unavailable an md5-gated local copy is used. Source data: DOI 10.6084/m9.figshare.6356420.v24 (MIT licence). """ import hashlib import json import os import sys import warnings from pathlib import Path import numpy as np import pandas as pd from scipy import stats # ---------------------------------------------------------------------------- # 0. Reproducibility # ---------------------------------------------------------------------------- SEED = 42 N_BOOT = 10000 np.random.seed(SEED) RNG = np.random.default_rng(SEED) # ---------------------------------------------------------------------------- # Pinned data provenance # ---------------------------------------------------------------------------- DATA_FILENAME = "HELIO4CAST_ICMECAT_v23.csv" DATA_MD5 = "d6ea054d1f6baf2311d14181f7e34e5b" DATA_URL = "https://ndownloader.figshare.com/files/58757791" LOCAL_FALLBACK = "/Users/markhahnel/datasetpapers-handfeed/helio4cast/" + DATA_FILENAME HERE = Path(__file__).resolve().parent FIG_DIR = HERE / "figures" TAB_DIR = HERE / "tables" for d in (FIG_DIR, TAB_DIR): d.mkdir(exist_ok=True) def _md5(path): h = hashlib.md5() with open(path, "rb") as fh: for chunk in iter(lambda: fh.read(1 << 20), b""): h.update(chunk) return h.hexdigest() def load_data(): """Fetch canonical figshare file (md5-verified); fall back to md5-gated local copy.""" local_cache = HERE / DATA_FILENAME # 1) Try the canonical download once, wrapped so any network failure is non-fatal. if not local_cache.exists(): try: import urllib.request print(f"[data] attempting canonical download: {DATA_URL}") urllib.request.urlretrieve(DATA_URL, local_cache) got = _md5(local_cache) if got != DATA_MD5: local_cache.unlink(missing_ok=True) raise ValueError(f"downloaded md5 {got} != expected {DATA_MD5}") print("[data] canonical download verified (md5 OK)") except Exception as exc: # noqa: BLE001 - network optional by design print(f"[data] canonical download unavailable ({exc}); using local fallback") if not os.path.exists(LOCAL_FALLBACK): raise FileNotFoundError( "Canonical download failed and no local fallback present at " f"{LOCAL_FALLBACK}" ) from exc got = _md5(LOCAL_FALLBACK) if got != DATA_MD5: # hard-fail on mismatch raise ValueError( f"Local fallback md5 {got} != expected {DATA_MD5}; refusing to proceed." ) from exc import shutil shutil.copyfile(LOCAL_FALLBACK, local_cache) print("[data] local fallback verified (md5 OK)") else: if _md5(local_cache) != DATA_MD5: raise ValueError("Cached file md5 mismatch; delete and re-run.") print("[data] using verified local cache") return pd.read_csv(local_cache) # ---------------------------------------------------------------------------- # Plot styling (self-contained, colourblind-safe Okabe-Ito palette) # ---------------------------------------------------------------------------- import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt plt.rcParams.update({ "figure.dpi": 150, "savefig.dpi": 200, "font.size": 9, "axes.titlesize": 10, "axes.labelsize": 9, "xtick.labelsize": 8, "ytick.labelsize": 8, "legend.fontsize": 7, "axes.spines.top": False, "axes.spines.right": False, "xtick.direction": "out", "ytick.direction": "out", "figure.constrained_layout.use": True, }) # Okabe-Ito colourblind-safe palette, one hue per spacecraft OKABE_ITO = [ "#0072B2", "#E69F00", "#009E73", "#CC79A7", "#56B4E9", "#D55E00", "#F0E442", "#999999", "#000000", "#8C564B", "#117733", ] def bh_correct(pvals): """Benjamini-Hochberg FDR correction. Returns corrected p-values in input order.""" p = np.asarray(pvals, float) n = len(p) order = np.argsort(p) ranked = p[order] * n / (np.arange(1, n + 1)) ranked = np.minimum.accumulate(ranked[::-1])[::-1] out = np.empty(n) out[order] = np.clip(ranked, 0, 1) return out def loglog_ols(r, b): """OLS of log10(b) on log10(r). Returns dict with slope, alpha, intercept, r2, n.""" lr, lb = np.log10(r), np.log10(b) res = stats.linregress(lr, lb) return { "slope": float(res.slope), "alpha": float(-res.slope), "intercept": float(res.intercept), "r2": float(res.rvalue ** 2), "p_slope": float(res.pvalue), "stderr": float(res.stderr), "n": int(len(r)), } def bootstrap_slope_ci(r, b, n_boot=N_BOOT, rng=RNG): """Percentile bootstrap 95% CI on the log-log OLS slope.""" lr, lb = np.log10(np.asarray(r)), np.log10(np.asarray(b)) n = len(lr) slopes = np.empty(n_boot) for i in range(n_boot): idx = rng.integers(0, n, n) slopes[i] = stats.linregress(lr[idx], lb[idx]).slope lo, hi = np.percentile(slopes, [2.5, 97.5]) return float(lo), float(hi) TABLE_DESCRIPTIONS = { "tbl-1-powerlaw-fit.csv": "Headline power-law fit: OLS and Theil-Sen estimates of the radial exponent alpha (B ~ r^-alpha) for mo_bmean vs heliocentric distance, with bootstrap 95% CIs.", "tbl-2-spacecraft-summary.csv": "Per-spacecraft sample sizes, heliocentric-distance range, and median magnetic-obstacle mean field.", "tbl-3-perspacecraft-slopes.csv": "Per-spacecraft log-log OLS slope/alpha and within-spacecraft Spearman correlation (spacecraft with n>=30), with BH-corrected p-values.", "tbl-4-loso-alpha.csv": "Leave-one-spacecraft-out global alpha refits; each row excludes one spacecraft.", "tbl-5-within-spacecraft.csv": "Within-spacecraft radial correlation for the three widest-baseline instruments (PSP, Juno, Ulysses).", } # Human-readable descriptions for known fields (Frictionless field.description). FIELD_DESCRIPTIONS = { "sc_insitu": "In-situ spacecraft name.", "n": "Number of ICMEs.", "r_min": "Minimum heliocentric distance in group (AU).", "r_max": "Maximum heliocentric distance in group (AU).", "r_median": "Median heliocentric distance in group (AU).", "r_span_AU": "Heliocentric-distance baseline, r_max - r_min (AU).", "b_median": "Median magnetic-obstacle mean field (nT).", "slope": "Log-log OLS slope (= -alpha).", "alpha": "Power-law exponent alpha in B ~ r^-alpha.", "r2": "Coefficient of determination of the log-log fit.", "spearman_rho": "Spearman rank correlation of B vs r.", "spearman_p": "Two-sided Spearman p-value (uncorrected).", "spearman_p_bh": "Benjamini-Hochberg FDR-corrected Spearman p-value.", "ols_slope": "Log-log OLS slope for this spacecraft (= -alpha).", "ols_alpha": "Log-log OLS power-law exponent for this spacecraft.", "method": "Estimation method.", "alpha_ci95_lo": "Lower bound of 95% CI on alpha.", "alpha_ci95_hi": "Upper bound of 95% CI on alpha.", "excluded_sc": "Spacecraft removed before refitting the global alpha.", } _FRICTIONLESS_TYPE = {"i": "integer", "f": "number", "O": "string", "b": "boolean"} def write_datapackage(): """Build a Frictionless tabular-data-package descriptor from the written CSVs.""" resources = [] for fname in sorted(TABLE_DESCRIPTIONS): path = TAB_DIR / fname if not path.exists(): continue d = pd.read_csv(path) fields = [] for col in d.columns: ftype = _FRICTIONLESS_TYPE.get(d[col].dtype.kind, "string") fields.append({ "name": col, "type": ftype, "description": FIELD_DESCRIPTIONS.get(col, ""), }) resources.append({ "name": fname.replace(".csv", "").replace(".", "-"), "path": fname, "format": "csv", "mediatype": "text/csv", "encoding": "utf-8", "title": fname, "description": TABLE_DESCRIPTIONS[fname], "schema": {"fields": fields}, }) package = { "profile": "tabular-data-package", "name": "icme-radial-field-scaling-tables", "title": "Radial scaling of ICME magnetic-obstacle field strength - derived tables", "description": ( "Derived tables for the pre-registered analysis of how ICME magnetic-obstacle " "mean field strength (mo_bmean) scales with heliocentric distance in the " "HELIO4CAST ICMECAT v2.3 catalogue." ), "licenses": [{"name": "MIT", "title": "MIT License"}], "sources": [{ "title": "HELIO4CAST Interplanetary Coronal Mass Ejection Catalog v2.3 (ICMECAT)", "path": "https://doi.org/10.6084/m9.figshare.6356420.v24", }], "resources": resources, } with open(TAB_DIR / "datapackage.json", "w") as fh: json.dump(package, fh, indent=2) def main(): results = {} # ---- Pre-registration block ---- prereg = { "question": ( "How does ICME magnetic-obstacle mean field strength (mo_bmean, nT) " "scale with spacecraft heliocentric distance (mo_sc_heliodistance, AU)? " "Estimate alpha in B ~ r^-alpha." ), "primary_test_T1": "Spearman rank correlation of mo_bmean vs mo_sc_heliodistance (two-sided).", "powerlaw_test_T2": "OLS of log10(mo_bmean) on log10(r); slope = -alpha; bootstrap 95% CI (n=10000).", "robust_test_T3": "Theil-Sen slope on log-log data as a sensitivity check to outliers/non-normality.", "confound_guards": [ "Per-spacecraft log-log OLS slopes (spacecraft with n>=30).", "Leave-one-spacecraft-out global alpha refit (range across 11 refits).", "Within-spacecraft Spearman for the three widest-baseline instruments PSP, Juno, and Ulysses.", ], "multiple_comparison_correction": "Benjamini-Hochberg FDR across per-spacecraft + within-spacecraft test family.", "expected": "Negative rho ~ -0.7; alpha ~ 1.5-1.9 consistent with published MO radial scaling.", "seed": SEED, "n_bootstrap": N_BOOT, "variables": { "response": "mo_bmean (nT)", "predictor": "mo_sc_heliodistance (AU)", "grouping": "sc_insitu", }, "test_direction": "two-sided", } results["preregistration"] = prereg # ---- Load & filter ---- df = load_data() n_total = len(df) sub = df[["mo_bmean", "mo_sc_heliodistance", "sc_insitu"]].copy() sub = sub[(sub["mo_bmean"] > 0) & (sub["mo_sc_heliodistance"] > 0)] sub = sub.dropna(subset=["mo_bmean", "mo_sc_heliodistance", "sc_insitu"]) n_kept = len(sub) results["data"] = { "source_doi": "10.6084/m9.figshare.6356420.v24", "source_file": DATA_FILENAME, "source_md5": DATA_MD5, "n_total_rows": int(n_total), "n_analysis_rows": int(n_kept), "n_dropped": int(n_total - n_kept), "heliodistance_min_AU": float(sub["mo_sc_heliodistance"].min()), "heliodistance_max_AU": float(sub["mo_sc_heliodistance"].max()), "n_spacecraft": int(sub["sc_insitu"].nunique()), } r = sub["mo_sc_heliodistance"].to_numpy() b = sub["mo_bmean"].to_numpy() # ---- T1: Spearman (primary) ---- sp = stats.spearmanr(r, b) results["T1_spearman"] = { "test": "Spearman rank correlation (mo_bmean vs mo_sc_heliodistance)", "statistic_rho": float(sp.statistic), "p_value": float(sp.pvalue), "effect_size": "rho", "effect_size_value": float(sp.statistic), "n": int(n_kept), } # ---- T2: OLS power-law ---- ols = loglog_ols(r, b) lo, hi = bootstrap_slope_ci(r, b) ols["slope_ci95"] = [lo, hi] ols["alpha_ci95"] = [float(-hi), float(-lo)] results["T2_powerlaw_ols"] = { "test": "OLS log10(mo_bmean) ~ log10(r)", "statistic_slope": ols["slope"], "alpha": ols["alpha"], "slope_ci95": ols["slope_ci95"], "alpha_ci95": ols["alpha_ci95"], "intercept": ols["intercept"], "r2": ols["r2"], "p_value": ols["p_slope"], "effect_size": "R2", "effect_size_value": ols["r2"], "n": ols["n"], } # ---- T3: Theil-Sen robust ---- lr, lb = np.log10(r), np.log10(b) ts = stats.theilslopes(lb, lr) # slope, intercept, lo_slope, hi_slope results["T3_theilsen"] = { "test": "Theil-Sen robust slope of log10(mo_bmean) ~ log10(r)", "statistic_slope": float(ts[0]), "alpha": float(-ts[0]), "intercept": float(ts[1]), "slope_ci95": [float(ts[2]), float(ts[3])], "alpha_ci95": [float(-ts[3]), float(-ts[2])], "effect_size": "robust_slope", "effect_size_value": float(ts[0]), "n": int(n_kept), } # ---- Per-spacecraft summary (tbl-2) ---- grp = sub.groupby("sc_insitu") sc_summary = grp.agg( n=("mo_bmean", "size"), r_min=("mo_sc_heliodistance", "min"), r_median=("mo_sc_heliodistance", "median"), r_max=("mo_sc_heliodistance", "max"), b_median=("mo_bmean", "median"), ).reset_index().sort_values("r_median") sc_summary.round(4).to_csv(TAB_DIR / "tbl-2-spacecraft-summary.csv", index=False) # ---- Guard (a): per-spacecraft log-log slopes, n>=30 ---- persc_rows = [] for name, g in grp: if len(g) >= 30: gr = g["mo_sc_heliodistance"].to_numpy() gb = g["mo_bmean"].to_numpy() fit = loglog_ols(gr, gb) # within-spacecraft Spearman too (for all n>=30 spacecraft) gsp = stats.spearmanr(gr, gb) persc_rows.append({ "sc_insitu": name, "n": fit["n"], "r_min": float(gr.min()), "r_max": float(gr.max()), "r_span_AU": float(gr.max() - gr.min()), "slope": fit["slope"], "alpha": fit["alpha"], "r2": fit["r2"], "spearman_rho": float(gsp.statistic), "spearman_p": float(gsp.pvalue), }) persc = pd.DataFrame(persc_rows).sort_values("r_span_AU", ascending=False) # ---- Guard (c): within-spacecraft Spearman for Ulysses & PSP ---- # Three widest-baseline instruments by heliocentric-distance span: # PSP (0.07-0.91 AU, widest in log-ratio), Juno (1.04-5.42 AU), Ulysses (1.34-5.41 AU). within_targets = ["PSP", "Juno", "ULYSSES"] within_rows = [] within_pvals = [] for name in within_targets: g = sub[sub["sc_insitu"] == name] if len(g) >= 10: gr = g["mo_sc_heliodistance"].to_numpy() gb = g["mo_bmean"].to_numpy() gsp = stats.spearmanr(gr, gb) fit = loglog_ols(gr, gb) within_rows.append({ "sc_insitu": name, "n": int(len(g)), "r_min": float(gr.min()), "r_max": float(gr.max()), "spearman_rho": float(gsp.statistic), "spearman_p": float(gsp.pvalue), "ols_slope": fit["slope"], "ols_alpha": fit["alpha"], "r2": fit["r2"], }) within_pvals.append(float(gsp.pvalue)) within = pd.DataFrame(within_rows) # ---- Multiple-comparison correction across the guard test family ---- # Family = per-spacecraft Spearman tests (n>=30) + within-spacecraft (Ulysses, PSP). fam_names, fam_p = [], [] for _, row in persc.iterrows(): fam_names.append(f"perSC_spearman_{row['sc_insitu']}") fam_p.append(row["spearman_p"]) for _, row in within.iterrows(): fam_names.append(f"within_spearman_{row['sc_insitu']}") fam_p.append(row["spearman_p"]) fam_p_corr = bh_correct(fam_p) family = pd.DataFrame({"test": fam_names, "p_raw": fam_p, "p_bh": fam_p_corr}) # attach corrected p back to per-spacecraft table bh_map = dict(zip(fam_names, fam_p_corr)) persc["spearman_p_bh"] = [bh_map[f"perSC_spearman_{s}"] for s in persc["sc_insitu"]] within["spearman_p_bh"] = [bh_map[f"within_spearman_{s}"] for s in within["sc_insitu"]] persc.round(6).to_csv(TAB_DIR / "tbl-3-perspacecraft-slopes.csv", index=False) within.round(6).to_csv(TAB_DIR / "tbl-5-within-spacecraft.csv", index=False) results["guard_perspacecraft"] = { "test": "Per-spacecraft log-log OLS slope + within-spacecraft Spearman (n>=30)", "spacecraft": persc.to_dict(orient="records"), "n_spacecraft_tested": int(len(persc)), "multiple_comparison_correction": "Benjamini-Hochberg FDR across full guard family", } results["guard_within_spacecraft"] = { "test": "Within-spacecraft Spearman (widest r-range instruments)", "targets": within.to_dict(orient="records"), } results["guard_family_bh"] = family.to_dict(orient="records") # ---- Guard (b): leave-one-spacecraft-out global alpha ---- loso_rows = [] for name in sorted(sub["sc_insitu"].unique()): keep = sub[sub["sc_insitu"] != name] fit = loglog_ols(keep["mo_sc_heliodistance"].to_numpy(), keep["mo_bmean"].to_numpy()) loso_rows.append({"excluded_sc": name, "alpha": fit["alpha"], "slope": fit["slope"], "n": fit["n"]}) loso = pd.DataFrame(loso_rows) loso.round(6).to_csv(TAB_DIR / "tbl-4-loso-alpha.csv", index=False) results["guard_loso"] = { "test": "Leave-one-spacecraft-out global alpha refit", "alpha_min": float(loso["alpha"].min()), "alpha_max": float(loso["alpha"].max()), "alpha_range": float(loso["alpha"].max() - loso["alpha"].min()), "full_alpha": ols["alpha"], "per_exclusion": loso.to_dict(orient="records"), } # ---- tbl-1: headline power-law fit ---- tbl1 = pd.DataFrame([{ "method": "OLS (log-log)", "slope": ols["slope"], "alpha": ols["alpha"], "alpha_ci95_lo": ols["alpha_ci95"][0], "alpha_ci95_hi": ols["alpha_ci95"][1], "r2": ols["r2"], "n": ols["n"], }, { "method": "Theil-Sen (robust)", "slope": float(ts[0]), "alpha": float(-ts[0]), "alpha_ci95_lo": float(-ts[3]), "alpha_ci95_hi": float(-ts[2]), "r2": np.nan, "n": int(n_kept), }]) tbl1.round(6).to_csv(TAB_DIR / "tbl-1-powerlaw-fit.csv", index=False) # ------------------------------------------------------------------ # Figures # ------------------------------------------------------------------ sc_order = sc_summary["sc_insitu"].tolist() color_map = {s: OKABE_ITO[i % len(OKABE_ITO)] for i, s in enumerate(sc_order)} # ---- fig-1: log-log scatter + fits ---- fig1, ax = plt.subplots(figsize=(6.4, 4.6)) for s in sc_order: g = sub[sub["sc_insitu"] == s] ax.scatter(g["mo_sc_heliodistance"], g["mo_bmean"], s=14, alpha=0.55, color=color_map[s], edgecolors="none", label=f"{s} (n={len(g)})") xs = np.logspace(np.log10(r.min()), np.log10(r.max()), 100) ax.plot(xs, 10 ** (ols["intercept"]) * xs ** ols["slope"], color="black", lw=2, label=f"OLS α={ols['alpha']:.2f}") ax.plot(xs, 10 ** (ts[1]) * xs ** ts[0], color="black", lw=1.5, ls="--", label=f"Theil–Sen α={-ts[0]:.2f}") ax.set_xscale("log"); ax.set_yscale("log") ax.set_xlabel("Heliocentric distance r (AU)") ax.set_ylabel("Magnetic-obstacle mean field B (nT)") ax.set_title("ICME magnetic-obstacle field strength declines with heliocentric distance") ax.legend(frameon=False, fontsize=6, ncol=2, loc="upper right") fig1.savefig(FIG_DIR / "fig-1-radial-bmean-powerlaw.png") plt.close(fig1) # ---- fig-2: per-spacecraft alpha vs radial baseline + LOSO alpha range ---- # A per-spacecraft slope is only meaningful when that spacecraft samples a # wide range of r. Plotting alpha against the spacecraft's radial baseline # (r-span) shows that every reliable-baseline spacecraft (r-span >= 0.30 AU) # yields a negative slope near the global value (alpha ~ 1.0-1.6), while # narrow-baseline spacecraft (near-fixed r) give ill-constrained slopes. fig2, (axa, axb) = plt.subplots(1, 2, figsize=(9.2, 4.4)) RELIABLE_SPAN = 0.30 # AU; below this a within-spacecraft slope is ill-constrained ps = persc.copy() reliable = ps["r_span_AU"] >= RELIABLE_SPAN axa.axhline(ols["alpha"], color="black", lw=1.5, zorder=1, label=f"global α={ols['alpha']:.2f}") axa.axhspan(-ols["slope_ci95"][1], -ols["slope_ci95"][0], color="grey", alpha=0.15, zorder=0, label="global α 95% CI") axa.axvline(RELIABLE_SPAN, color="grey", lw=1, ls=":", zorder=1) axa.scatter(ps.loc[reliable, "r_span_AU"], ps.loc[reliable, "alpha"], color=[color_map[s] for s in ps.loc[reliable, "sc_insitu"]], s=70, zorder=3, edgecolors="black", linewidths=0.6) axa.scatter(ps.loc[~reliable, "r_span_AU"], ps.loc[~reliable, "alpha"], facecolors="none", edgecolors="grey", s=55, zorder=2, label=f"r-span < {RELIABLE_SPAN} AU (ill-constrained)") # stagger labels vertically (in x-order) so the mid-range cluster doesn't collide rel_sorted = ps[reliable].sort_values("r_span_AU").reset_index(drop=True) stagger = [10, 22, -16, -28, 10, 10] for i, row in rel_sorted.iterrows(): dy = stagger[i % len(stagger)] axa.annotate(row["sc_insitu"], (row["r_span_AU"], row["alpha"]), textcoords="offset points", xytext=(6, dy), fontsize=7, arrowprops=dict(arrowstyle="-", lw=0.5, color="grey")) axa.set_xlabel("Spacecraft radial baseline (r-span, AU)") axa.set_ylabel("Per-spacecraft power-law exponent α") axa.set_title("Per-spacecraft exponent vs radial baseline:\nreliable-baseline slopes span α ≈ 1.0–1.6") axa.legend(frameon=False, fontsize=7, loc="lower right") lo_a = loso.sort_values("alpha") yl = np.arange(len(lo_a)) axb.scatter(lo_a["alpha"], yl, color="#0072B2", s=40, zorder=3, edgecolors="black", linewidths=0.6) axb.axvline(ols["alpha"], color="black", lw=1.5, label=f"all-data α={ols['alpha']:.2f}") axb.set_yticks(yl); axb.set_yticklabels([f"–{s}" for s in lo_a["excluded_sc"]]) axb.set_xlabel("Global α with one spacecraft removed") axb.set_title("Leave-one-spacecraft-out stability") axb.legend(frameon=False, fontsize=7, loc="lower right") fig2.savefig(FIG_DIR / "fig-2-confound-guards.png") plt.close(fig2) # ---- fig-3: within-spacecraft (Ulysses & PSP) ---- fig3, axes = plt.subplots(1, len(within_targets), figsize=(4.0 * len(within_targets), 4.0), squeeze=False) for j, name in enumerate(within_targets): axx = axes[0][j] g = sub[sub["sc_insitu"] == name] gr = g["mo_sc_heliodistance"].to_numpy(); gb = g["mo_bmean"].to_numpy() axx.scatter(gr, gb, s=16, alpha=0.6, color=color_map.get(name, "#0072B2"), edgecolors="none") fit = loglog_ols(gr, gb) gsp = stats.spearmanr(gr, gb) xr = np.logspace(np.log10(gr.min()), np.log10(gr.max()), 50) axx.plot(xr, 10 ** fit["intercept"] * xr ** fit["slope"], color="black", lw=1.8) axx.set_xscale("log"); axx.set_yscale("log") axx.set_xlabel("Heliocentric distance r (AU)") if j == 0: axx.set_ylabel("Magnetic-obstacle mean field B (nT)") axx.set_title(f"{name} (n={len(g)})\nα={fit['alpha']:.2f}, ρ={gsp.statistic:.2f}") fig3.suptitle("Radial decline holds within a single spacecraft", fontsize=10) fig3.savefig(FIG_DIR / "fig-3-within-spacecraft.png") plt.close(fig3) # ------------------------------------------------------------------ # tables/datapackage.json (Frictionless) # ------------------------------------------------------------------ write_datapackage() # results.json # ------------------------------------------------------------------ with open(HERE / "results.json", "w") as fh: json.dump(results, fh, indent=2) print("[done] results.json + figures + tables + datapackage written") print(f" T1 Spearman rho = {sp.statistic:.3f}, p = {sp.pvalue:.2e}, n = {n_kept}") print(f" T2 OLS alpha = {ols['alpha']:.3f} 95% CI [{ols['alpha_ci95'][0]:.3f}, {ols['alpha_ci95'][1]:.3f}]") print(f" T3 Theil-Sen alpha = {-ts[0]:.3f}") print(f" LOSO alpha range = [{loso['alpha'].min():.3f}, {loso['alpha'].max():.3f}]") return results if __name__ == "__main__": with warnings.catch_warnings(): warnings.simplefilter("ignore") main()