07.5 周期ごとの関係：クロススペクトル・コヒーレンス・位相差

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import signal

# ============================================================
# 0. 設定
# ============================================================
soi_file = "SOI_timeseries_updated.txt"
sst_file = "sstoi.indices"

outdir = "fig_cross_spectrum_soi_nino34"
os.makedirs(outdir, exist_ok=True)

# 月平均データなので、1年あたり12点
fs = 12

# Welch法の窓長
# 256か月 = 約21.3年
nperseg = 256
noverlap = nperseg // 2

# 表示する周期範囲
period_min = 1.0
period_max = 10.0

# ENSO周期帯
enso_min = 2.0
enso_max = 7.0


def savefig(fig, filename_base):
    png_path = os.path.join(outdir, filename_base + ".png")
    pdf_path = os.path.join(outdir, filename_base + ".pdf")
    fig.tight_layout()
    fig.savefig(png_path, dpi=300, bbox_inches="tight")
    fig.savefig(pdf_path, bbox_inches="tight")
    print("Saved:", png_path)
    print("Saved:", pdf_path)
    plt.show()


# ============================================================
# 1. データ読み込み
# ============================================================
soi = pd.read_csv(soi_file, sep=r"\s+")

cols = [
    "YEAR", "MONTH",
    "NINO12", "ANOM12",
    "NINO3",  "ANOM3",
    "NINO4",  "ANOM4",
    "NINO34", "ANOM34"
]

sst = pd.read_csv(
    sst_file,
    sep=r"\s+",
    skiprows=1,
    names=cols
)

for c in cols:
    sst[c] = pd.to_numeric(sst[c], errors="coerce")

for c in ["YEAR", "MONTH", "SOI", "TIME"]:
    soi[c] = pd.to_numeric(soi[c], errors="coerce")


# ============================================================
# 2. 年月で結合
# ============================================================
df = pd.merge(
    soi,
    sst,
    on=["YEAR", "MONTH"],
    how="inner"
)

df = df.dropna(subset=["SOI", "ANOM34", "TIME"]).reset_index(drop=True)

print("Number of months =", len(df))
print("Start:", int(df["YEAR"].iloc[0]), int(df["MONTH"].iloc[0]))
print("End  :", int(df["YEAR"].iloc[-1]), int(df["MONTH"].iloc[-1]))

if len(df) < nperseg:
    raise ValueError("Data length is shorter than nperseg. Use smaller nperseg, e.g. 128.")


# ============================================================
# 3. 解析に使う時系列
# ============================================================
time = df["TIME"].to_numpy()
x_raw = df["SOI"].to_numpy()
y_raw = df["ANOM34"].to_numpy()

# 平均を引く
x = x_raw - np.nanmean(x_raw)
y = y_raw - np.nanmean(y_raw)

# 標準化
x = x / np.nanstd(x)
y = y / np.nanstd(y)

good = np.isfinite(x) & np.isfinite(y) & np.isfinite(time)
x = x[good]
y = y[good]
time = time[good]

r0 = np.corrcoef(x, y)[0, 1]
print("Correlation between SOI and Niño3.4 SST anomaly =", r0)


# ============================================================
# 4. パワースペクトル、クロススペクトル、コヒーレンス
# ============================================================
f, Pxx = signal.welch(
    x, fs=fs, window="hann",
    nperseg=nperseg, noverlap=noverlap,
    detrend="constant", scaling="density"
)

f, Pyy = signal.welch(
    y, fs=fs, window="hann",
    nperseg=nperseg, noverlap=noverlap,
    detrend="constant", scaling="density"
)

# クロススペクトル：Pxy は複素数
f, Pxy = signal.csd(
    x, y, fs=fs, window="hann",
    nperseg=nperseg, noverlap=noverlap,
    detrend="constant", scaling="density"
)

# コヒーレンス
f, coh = signal.coherence(
    x, y, fs=fs, window="hann",
    nperseg=nperseg, noverlap=noverlap,
    detrend="constant"
)

valid = f > 0
f2 = f[valid]
period = 1.0 / f2

Pxx2 = Pxx[valid]
Pyy2 = Pyy[valid]
Pxy2 = Pxy[valid]
coh2 = coh[valid]

cross_power = np.abs(Pxy2)
phase_rad = np.angle(Pxy2)
phase_deg = np.degrees(phase_rad)

phase_unwrap_rad = np.unwrap(phase_rad)
phase_unwrap_deg = np.degrees(phase_unwrap_rad)

lag_year = phase_rad / (2.0 * np.pi * f2)
lag_month = lag_year * 12.0

m = (period >= period_min) & (period <= period_max)
period_plot = period[m]
Pxx_plot = Pxx2[m]
Pyy_plot = Pyy2[m]
cross_plot = cross_power[m]
coh_plot = coh2[m]
phase_plot = phase_deg[m]
phase_unwrap_plot = phase_unwrap_deg[m]
lag_plot = lag_month[m]

enso_band_all = (period >= enso_min) & (period <= enso_max)
mean_coh_enso = np.nanmean(coh2[enso_band_all])
max_coh_enso = np.nanmax(coh2[enso_band_all])
period_max_coh = period[enso_band_all][np.nanargmax(coh2[enso_band_all])]

print(f"Mean coherence in ENSO band ({enso_min:g}–{enso_max:g} years): {mean_coh_enso:.3f}")
print(f"Max coherence in ENSO band ({enso_min:g}–{enso_max:g} years): {max_coh_enso:.3f}")
print(f"Period of max coherence in ENSO band: {period_max_coh:.2f} years")


# ============================================================
# 5. 図1：時系列
# ============================================================
fig, ax = plt.subplots(figsize=(12, 4))
ax.plot(time, x, label="SOI (standardized)", linewidth=1.2)
ax.plot(time, y, label="Niño3.4 SST anomaly (standardized)", linewidth=1.2)
ax.axhline(0, color="k", linewidth=0.8)
ax.set_xlabel("Year")
ax.set_ylabel("Standardized anomaly")
ax.set_title("SOI and Niño3.4 SST anomaly")
ax.legend(loc="lower right")
ax.grid(True, alpha=0.3)
ax.text(0.02, 0.95, f"r = {r0:.2f}", transform=ax.transAxes,
        va="top", ha="left",
        bbox=dict(facecolor="white", edgecolor="0.7", alpha=0.8))
savefig(fig, "fig01_timeseries_soi_nino34")


# ============================================================
# 6. 図2：パワースペクトル
# ============================================================
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(period_plot, Pxx_plot, label="SOI", linewidth=1.8)
ax.plot(period_plot, Pyy_plot, label="Niño3.4 SST anomaly", linewidth=1.8)
ax.axvspan(enso_min, enso_max, alpha=0.15, label="ENSO band (2–7 years)")
ax.set_xlabel("Period (years)")
ax.set_ylabel("Power spectral density")
ax.set_title("Power spectra")
ax.set_xlim(period_max, period_min)
ax.legend(loc="upper left")
ax.grid(True, alpha=0.3)
savefig(fig, "fig02_power_spectra")


# ============================================================
# 7. 図3：クロススペクトル
# ============================================================
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(period_plot, cross_plot, linewidth=1.8)
ax.axvspan(enso_min, enso_max, alpha=0.15, label="ENSO band (2–7 years)")
ax.set_xlabel("Period (years)")
ax.set_ylabel("|Cross spectrum|")
ax.set_title("Cross spectrum: SOI vs Niño3.4 SST anomaly")
ax.set_xlim(period_max, period_min)
ax.legend(loc="upper left")
ax.grid(True, alpha=0.3)
savefig(fig, "fig03_cross_spectrum")


# ============================================================
# 8. 図4：コヒーレンス
# ============================================================
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(period_plot, coh_plot, linewidth=1.8)
ax.axvspan(enso_min, enso_max, alpha=0.15, label="ENSO band (2–7 years)")
ax.set_xlabel("Period (years)")
ax.set_ylabel("Coherence")
ax.set_title("Coherence: SOI vs Niño3.4 SST anomaly")
ax.set_xlim(period_max, period_min)
ax.set_ylim(0, 1.05)
ax.legend(loc="upper left")
ax.grid(True, alpha=0.3)
text = (f"Mean coherence ({enso_min:g}–{enso_max:g} yr) = {mean_coh_enso:.2f}\n"
        f"Max coherence = {max_coh_enso:.2f}\n"
        f"Peak period = {period_max_coh:.2f} yr")
ax.text(0.98, 0.05, text, transform=ax.transAxes,
        va="bottom", ha="right",
        bbox=dict(facecolor="white", edgecolor="0.7", alpha=0.85))
savefig(fig, "fig04_coherence")


# ============================================================
# 9. 図5：位相差
# ============================================================
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(period_plot, phase_plot, linewidth=1.6)
ax.axhline(0, color="k", linewidth=0.8)
ax.axhline(180, color="gray", linestyle="--", linewidth=0.8)
ax.axhline(-180, color="gray", linestyle="--", linewidth=0.8)
ax.axvspan(enso_min, enso_max, alpha=0.15, label="ENSO band (2–7 years)")
ax.set_xlabel("Period (years)")
ax.set_ylabel("Phase difference (degree)")
ax.set_title("Phase difference: SOI vs Niño3.4 SST anomaly")
ax.set_xlim(period_max, period_min)
ax.set_ylim(-190, 190)
ax.legend(loc="upper left")
ax.grid(True, alpha=0.3)
ax.text(0.98, 0.05,
        "Note: SOI and Niño3.4 SST anomaly\nare mostly opposite in sign.\nPhase near ±180° indicates anti-phase.",
        transform=ax.transAxes, va="bottom", ha="right",
        bbox=dict(facecolor="white", edgecolor="0.7", alpha=0.85))
savefig(fig, "fig05_phase_raw")


# ============================================================
# 10. 図6：まとめ図
# ============================================================
fig, axes = plt.subplots(3, 1, figsize=(10, 12), sharex=False)

axes[0].plot(period_plot, Pxx_plot, label="SOI", linewidth=1.6)
axes[0].plot(period_plot, Pyy_plot, label="Niño3.4 SST anomaly", linewidth=1.6)
axes[0].axvspan(enso_min, enso_max, alpha=0.15)
axes[0].set_ylabel("Power spectral density")
axes[0].set_title("Power spectra")
axes[0].set_xlim(period_max, period_min)
axes[0].legend(loc="upper left")
axes[0].grid(True, alpha=0.3)

axes[1].plot(period_plot, cross_plot, linewidth=1.6)
axes[1].axvspan(enso_min, enso_max, alpha=0.15)
axes[1].set_ylabel("|Cross spectrum|")
axes[1].set_title("Cross spectrum")
axes[1].set_xlim(period_max, period_min)
axes[1].grid(True, alpha=0.3)

axes[2].plot(period_plot, coh_plot, linewidth=1.8)
axes[2].axvspan(enso_min, enso_max, alpha=0.15, label="ENSO band (2–7 years)")
axes[2].set_xlabel("Period (years)")
axes[2].set_ylabel("Coherence")
axes[2].set_title("Coherence")
axes[2].set_xlim(period_max, period_min)
axes[2].set_ylim(0, 1.05)
axes[2].legend(loc="upper left")
axes[2].grid(True, alpha=0.3)
axes[2].text(0.98, 0.05,
             f"Mean coherence = {mean_coh_enso:.2f}\nMax coherence = {max_coh_enso:.2f}\nPeak = {period_max_coh:.2f} yr",
             transform=axes[2].transAxes,
             va="bottom", ha="right",
             bbox=dict(facecolor="white", edgecolor="0.7", alpha=0.85))
savefig(fig, "fig08_summary_power_cross_coherence")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pywt

# =========================================================
# 0. 設定
# =========================================================
soi_file = "SOI_timeseries_updated.txt"
nino_file = "sstoi.indices"

fig_nino_name = "wavelet_nino34_improved.png"
fig_soi_name  = "wavelet_soi_improved.png"

# 月平均データ：dt は年単位
# 1か月 = 1/12年
dt_year = 1.0 / 12.0

# モール波
wavelet_name = "morl"

# スケール範囲
# 大きくすると長周期まで見ることができる
scales = np.arange(1, 256)

# ENSO帯（年）
enso_min_year = 2.0
enso_max_year = 7.0


# =========================================================
# 1. データ読み込み
# =========================================================
colnames = [
    "YR", "MON",
    "NINO12", "ANOM12",
    "NINO3",  "ANOM3",
    "NINO4",  "ANOM4",
    "NINO34", "ANOM34"
]

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

nino_time = nino_df["YR"].values + (nino_df["MON"].values - 1) / 12.0
nino_x = nino_df["ANOM34"].values.astype(float)

soi_df = pd.read_csv(soi_file, sep=r"\s+")
soi_time = soi_df["TIME"].values
soi_x = soi_df["SOI"].values.astype(float)


# =========================================================
# 2. 前処理関数
# =========================================================
def standardize(x):
    x = np.asarray(x, dtype=float)
    if np.any(~np.isfinite(x)):
        x = pd.Series(x).interpolate(limit_direction="both").values
    x = x - np.mean(x)
    x = x / np.std(x)
    return x


# =========================================================
# 3. Wavelet図作成関数
# =========================================================
def wavelet_plot_with_coi(x, time, title, ylabel, save_name,
                          wavelet_name="morl", scales=None,
                          enso_min_year=2.0, enso_max_year=7.0):
    x = standardize(x)

    coef, freqs = pywt.cwt(
        x,
        scales,
        wavelet_name,
        sampling_period=dt_year
    )

    power = np.abs(coef) ** 2

    # 周期（年）
    period_year = 1.0 / freqs

    # COIの簡易的な目安
    # 端に近いほど長周期成分は信頼しにくい
    n = len(x)
    edge_dist = np.minimum(np.arange(n), np.arange(n)[::-1])
    coi = np.sqrt(2) * edge_dist * dt_year

    fig, ax = plt.subplots(figsize=(14, 7))

    levels = np.linspace(0, np.nanpercentile(power, 99), 20)
    cf = ax.contourf(time, period_year, power, levels=levels, cmap="RdBu_r", extend="max")

    # ENSO帯を線で示す
    ax.axhline(enso_min_year, color="k", linestyle="--", linewidth=1.0)
    ax.axhline(enso_max_year, color="k", linestyle="--", linewidth=1.0)
    ax.text(time[0] + 0.02 * (time[-1] - time[0]),
            (enso_min_year + enso_max_year) / 2,
            "ENSO band (2–7 years)",
            va="center", ha="left",
            bbox=dict(facecolor="white", edgecolor="none", alpha=0.8))

    # COIを描く
    ax.plot(time, coi, color="gray", linewidth=1.5, label="COI")
    ax.fill_between(time, coi, period_year.max(), color="gray", alpha=0.25)

    ax.set_yscale("log")
    ax.set_ylim(period_year.max(), period_year.min())
    ax.set_ylabel(ylabel)
    ax.set_xlabel("Time")
    ax.set_title(title, fontsize=16)
    ax.grid(True, alpha=0.25)
    ax.legend(loc="upper right")

    cbar = fig.colorbar(cf, ax=ax, pad=0.02)
    cbar.set_label("Power")

    fig.savefig(save_name, dpi=300, bbox_inches="tight", facecolor="white")
    plt.show()
    print(f"Saved: {save_name}")


# =========================================================
# 4. 実行
# =========================================================
wavelet_plot_with_coi(
    x=nino_x,
    time=nino_time,
    title="Wavelet Power Spectrum (Niño3.4)",
    ylabel="Period (years)",
    save_name=fig_nino_name,
    wavelet_name=wavelet_name,
    scales=scales,
    enso_min_year=enso_min_year,
    enso_max_year=enso_max_year
)

wavelet_plot_with_coi(
    x=soi_x,
    time=soi_time,
    title="Wavelet Power Spectrum (SOI)",
    ylabel="Period (years)",
    save_name=fig_soi_name,
    wavelet_name=wavelet_name,
    scales=scales,
    enso_min_year=enso_min_year,
    enso_max_year=enso_max_year
)