08 ENSOコンポジット解析

穴埋めの要点

使う列: ANOM34
閾値: ±0.5℃
平滑化: 5か月移動平均
継続条件: 5か月以上
対応付け: YYYY-MM 単位
差分: comp_el - comp_la

完全スクリプト

import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt

sst_file = "sst.mnmean.nc"
nino_file = "sstoi.indices"
threshold = 0.5
min_len = 5
fig_ts_name   = "enso_event_detection.png"
fig_comp_name = "enso_composite_sst.png"
fig_anom_name = "enso_composite_sst_anomaly.png"

colnames = [
    "YR", "MON",
    "NINO12", "ANOM12",
    "NINO3",  "ANOM3",
    "NINO4",  "ANOM4",
    "NINO34", "ANOM34"
]

df = pd.read_csv(
    nino_file,
    sep=r"\s+",
    skiprows=1,
    names=colnames
)

df["time"] = pd.to_datetime(dict(year=df["YR"], month=df["MON"], day=15))
df["ym"] = df["time"].dt.to_period("M")
nino34_anom = pd.Series(df["ANOM34"].values, index=df["time"])
nino34_smooth = nino34_anom.rolling(window=5, center=True).mean()

def detect_runs(series, threshold=0.5, mode="el", min_len=5):
    if mode == "el":
        mask = series >= threshold
    elif mode == "la":
        mask = series <= -threshold
    else:
        raise ValueError("mode must be 'el' or 'la'")

    events = np.zeros(len(series), dtype=bool)
    count = 0

    for i in range(len(series)):
        if pd.notna(mask.iloc[i]) and mask.iloc[i]:
            count += 1
        else:
            if count >= min_len:
                events[i-count:i] = True
            count = 0

    if count >= min_len:
        events[len(series)-count:len(series)] = True

    return pd.Series(events, index=series.index)

el_event = detect_runs(nino34_smooth, threshold=threshold, mode="el", min_len=min_len)
la_event = detect_runs(nino34_smooth, threshold=threshold, mode="la", min_len=min_len)

el_years = np.unique(el_event.index[el_event].year)
la_years = np.unique(la_event.index[la_event].year)
print("El Niño years :", el_years)
print("La Niña years :", la_years)

plt.figure(figsize=(14, 5))
plt.plot(nino34_anom.index, nino34_anom.values, color="0.65", lw=1.0, label="Nino3.4 anomaly (monthly)")
plt.plot(nino34_smooth.index, nino34_smooth.values, color="k", lw=2.0, label="5-month running mean")
plt.axhline(threshold, color="r", ls="--", lw=1)
plt.axhline(-threshold, color="b", ls="--", lw=1)
for t in el_event.index[el_event]:
    plt.axvspan(t - pd.Timedelta(days=15), t + pd.Timedelta(days=15), color="red", alpha=0.08)
for t in la_event.index[la_event]:
    plt.axvspan(t - pd.Timedelta(days=15), t + pd.Timedelta(days=15), color="blue", alpha=0.08)
plt.title("Niño3.4 anomaly and ENSO event detection")
plt.ylabel("Niño3.4 anomaly (°C)")
plt.xlabel("Time")
plt.legend(loc="upper right")
plt.tight_layout()
plt.savefig(fig_ts_name, dpi=300, bbox_inches="tight", facecolor="white")
plt.show()

ds = xr.open_dataset(sst_file)
sst = ds["sst"].astype(float)
sst = sst.where(sst > -100)
time_sst = pd.to_datetime(ds["time"].values)
ym_sst = pd.Series(time_sst).dt.to_period("M")

event_df = pd.DataFrame({
    "ym": df["ym"],
    "el": el_event.values,
    "la": la_event.values
})

el_map = dict(zip(event_df["ym"], event_df["el"]))
la_map = dict(zip(event_df["ym"], event_df["la"]))
el_mask_sst = np.array([el_map.get(ym, False) for ym in ym_sst], dtype=bool)
la_mask_sst = np.array([la_map.get(ym, False) for ym in ym_sst], dtype=bool)

sst_el = sst.isel(time=el_mask_sst)
sst_la = sst.isel(time=la_mask_sst)
comp_el = sst_el.mean(dim="time")
comp_la = sst_la.mean(dim="time")
comp_diff = comp_el - comp_la
clim_all = sst.mean(dim="time")
comp_el_anom = comp_el - clim_all
comp_la_anom = comp_la - clim_all

fig, axes = plt.subplots(1, 3, figsize=(18, 5), constrained_layout=True)
pcm1 = axes[0].pcolormesh(comp_el["lon"], comp_el["lat"], comp_el, shading="auto", cmap="turbo")
axes[0].set_title("El Niño Composite SST")
axes[0].set_xlabel("Longitude")
axes[0].set_ylabel("Latitude")
fig.colorbar(pcm1, ax=axes[0], shrink=0.9, label="SST (°C)")
pcm2 = axes[1].pcolormesh(comp_la["lon"], comp_la["lat"], comp_la, shading="auto", cmap="turbo")
axes[1].set_title("La Niña Composite SST")
axes[1].set_xlabel("Longitude")
axes[1].set_ylabel("Latitude")
fig.colorbar(pcm2, ax=axes[1], shrink=0.9, label="SST (°C)")
pcm3 = axes[2].pcolormesh(comp_diff["lon"], comp_diff["lat"], comp_diff, shading="auto", cmap="RdBu_r", vmin=-2, vmax=2)
axes[2].set_title("El Niño - La Niña")
axes[2].set_xlabel("Longitude")
axes[2].set_ylabel("Latitude")
fig.colorbar(pcm3, ax=axes[2], shrink=0.9, label="SST difference (°C)")
plt.savefig(fig_comp_name, dpi=300, bbox_inches="tight", facecolor="white")
plt.show()

fig, axes = plt.subplots(1, 3, figsize=(18, 5), constrained_layout=True)
pcm1 = axes[0].pcolormesh(comp_el_anom["lon"], comp_el_anom["lat"], comp_el_anom, shading="auto", cmap="RdBu_r", vmin=-2, vmax=2)
axes[0].set_title("El Niño Composite SST Anomaly")
axes[0].set_xlabel("Longitude")
axes[0].set_ylabel("Latitude")
fig.colorbar(pcm1, ax=axes[0], shrink=0.9, label="SST anomaly (°C)")
pcm2 = axes[1].pcolormesh(comp_la_anom["lon"], comp_la_anom["lat"], comp_la_anom, shading="auto", cmap="RdBu_r", vmin=-2, vmax=2)
axes[1].set_title("La Niña Composite SST Anomaly")
axes[1].set_xlabel("Longitude")
axes[1].set_ylabel("Latitude")
fig.colorbar(pcm2, ax=axes[1], shrink=0.9, label="SST anomaly (°C)")
pcm3 = axes[2].pcolormesh(comp_diff["lon"], comp_diff["lat"], comp_diff, shading="auto", cmap="RdBu_r", vmin=-2, vmax=2)
axes[2].set_title("El Niño - La Niña")
axes[2].set_xlabel("Longitude")
axes[2].set_ylabel("Latitude")
fig.colorbar(pcm3, ax=axes[2], shrink=0.9, label="SST difference (°C)")
plt.savefig(fig_anom_name, dpi=300, bbox_inches="tight", facecolor="white")
plt.show()