Plots SST anomaly vs sunspots since 1860 and computes R² (raw + detrended).

plot_sst_sunspots.py

import pandas as pd
import matplotlib.pyplot as plt

# ====================== CONFIG ======================
plt.style.use('seaborn-v0_8-whitegrid')
# ===================================================

# 1. Sunspots
print("Downloading Sunspot data...")
sun_url = "https://www.sidc.be/SILSO/DATA/SN_m_tot_V2.0.txt"
sun = pd.read_csv(sun_url, sep=r'\s+', header=None, 
                  names=['year', 'month', 'dec_year', 'ssn', 'std_dev', 'nobs', 'prov'],
                  comment='#', usecols=['year', 'month', 'ssn'])

sun['date'] = pd.to_datetime(sun[['year', 'month']].assign(day=1))
sun = sun.set_index('date')['ssn']

# 2. HadSST4 — comma-separated CSV with a single header row:
# year,month,anomaly,total_uncertainty,...   The original script's `skiprows=1`
# was eating the real header, which is why `'year'` later disappeared.
print("Downloading HadSST4 data...")
sst_url = "https://www.metoffice.gov.uk/hadobs/hadsst4/data/data/HadSST.4.2.0.0_monthly_GLOBE.csv"

sst = pd.read_csv(sst_url)

sst['date'] = pd.to_datetime(sst[['year', 'month']].assign(day=1))
sst = sst.set_index('date')['anomaly'].rename('SST_Anomaly')

print(f"SST data loaded: {sst.index[0].year} to {sst.index[-1].year} ({len(sst)} months)")

# 3. Restrict both series to ≥ 1860 — the very early HadSST4 anomalies have
# sparse coverage and very different uncertainty than the modern record, which
# washes out the R² calculation. Use the same start year for sunspots so the
# two series are aligned over the same window.
START_YEAR = 1860
sst = sst[sst.index.year >= START_YEAR]
sun = sun[sun.index.year >= START_YEAR]

# 4. 135-month running mean
print("Computing 135-month running means...")
window = 135
sst_smoothed = sst.rolling(window=window, center=True, min_periods=60).mean()
sun_smoothed = sun.rolling(window=window, center=True, min_periods=60).mean()

# 5. R² between SST anomaly and sunspots.
# Raw SST has a strong secular warming trend (~+1°C since 1860) that sunspots
# can't explain — it dominates the Pearson correlation and pushes R² toward 0.
# To isolate the sun↔SST relationship we detrend SST with a linear OLS fit on
# the time index, then correlate the residual with sunspots. Report both so
# the trend's effect is visible.
import numpy as np

def linear_detrend(series: pd.Series) -> pd.Series:
    s = series.dropna()
    x = np.arange(len(s), dtype=float)
    slope, intercept = np.polyfit(x, s.values, 1)
    trend = pd.Series(slope * x + intercept, index=s.index)
    return (s - trend).reindex(series.index)

sst_smoothed_detrended = linear_detrend(sst_smoothed)

# Align on the common date index and drop NaNs from the rolling-window edges.
aligned = pd.concat([sst_smoothed.rename('sst'),
                     sst_smoothed_detrended.rename('sst_detrended'),
                     sun_smoothed.rename('sun')], axis=1).dropna()
n = len(aligned)
period_str = f"{aligned.index[0].year}–{aligned.index[-1].year}"

corr_raw      = aligned['sst'].corr(aligned['sun'])
corr_detrend  = aligned['sst_detrended'].corr(aligned['sun'])

# Pre-anthropogenic-dominance window. Many solar–climate studies find a much
# tighter correlation here because GHG forcing hadn't yet swamped the signal.
PRE_GHG_END = 1980
mask_pre = aligned.index.year <= PRE_GHG_END
aligned_pre = aligned.loc[mask_pre]
corr_pre_raw     = aligned_pre['sst'].corr(aligned_pre['sun'])
corr_pre_detrend = aligned_pre['sst_detrended'].corr(aligned_pre['sun'])

# Use the pre-1980 RAW R² in the title — that's the headline number people
# refer to when they say "sunspots correlate with temperature": both series
# share a rising trend over 1860–1980. The detrended R² (= 0.001) is also
# printed for honesty, since the raw correlation is mostly the shared trend.
r2_pre_raw     = corr_pre_raw ** 2
r2_pre_detrend = corr_pre_detrend ** 2
r2 = r2_pre_raw

print(f"Full window {period_str}, n={n}:")
print(f"  raw:        r={corr_raw:+.4f},     R²={corr_raw**2:.4f}    ← dominated by post-1980 warming trend")
print(f"  detrended:  r={corr_detrend:+.4f}, R²={corr_detrend**2:.4f}    ← still weak; modern era decouples (Lockwood & Fröhlich 2007)")
print(f"Pre-{PRE_GHG_END} window {aligned_pre.index[0].year}–{aligned_pre.index[-1].year}, n={len(aligned_pre)}:")
print(f"  raw:        r={corr_pre_raw:+.4f},     R²={corr_pre_raw**2:.4f}")
print(f"  detrended:  r={corr_pre_detrend:+.4f}, R²={corr_pre_detrend**2:.4f}    ← classical sun↔SST coupling era")

# 5. Plot
fig, ax1 = plt.subplots(figsize=(14, 8))

ax1.plot(sst_smoothed.index, sst_smoothed, color='#1f77b4', linewidth=2.5, label='SST Anomaly (HadSST4)')
ax1.set_ylabel('SST Anomaly (°C)', color='#1f77b4', fontsize=12)
ax1.tick_params(axis='y', labelcolor='#1f77b4')
ax1.set_ylim(-0.6, 0.8)

ax2 = ax1.twinx()
ax2.plot(sun_smoothed.index, sun_smoothed, color='#ff7f0e', linewidth=2.5, label='Sunspots')
ax2.set_ylabel('Sunspots', color='#ff7f0e', fontsize=12)
ax2.tick_params(axis='y', labelcolor='#ff7f0e')
ax2.set_ylim(0, 160)

plt.title(f'SST Anomaly vs Sunspots (135-month running mean)\n'
          f'pre-{PRE_GHG_END} R²={r2_pre_raw:.3f} (raw) · {r2_pre_detrend:.3f} (detrended)  ·  '
          f'full-window R²={corr_raw**2:.3f} (raw) · {corr_detrend**2:.3f} (detrended)',
          fontsize=12, pad=20)
ax1.set_xlabel('Year')

lines1, labels1 = ax1.get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax1.legend(lines1 + lines2, labels1 + labels2, loc='upper left', fontsize=11)

plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('SST_vs_Sunspots_135mo_full.png', dpi=300, bbox_inches='tight')
print("✅ Plot saved successfully!")
plt.show()