02 南方振動指数の計算 解答

完全スクリプト

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

base_start = 1981
base_end   = 2010
months = ["JAN","FEB","MAR","APR","MAY","JUN",
          "JUL","AUG","SEP","OCT","NOV","DEC"]

dar = pd.read_csv("Darwin.txt", sep=r"\s+", engine="python", na_values=[-99.9])
tah = pd.read_csv("Tahiti.txt", sep=r"\s+", engine="python", na_values=[-99.9])

dar = dar.set_index("YEAR")
tah = tah.set_index("YEAR")

dar_base = dar.loc[base_start:base_end, months]
tah_base = tah.loc[base_start:base_end, months]

mslpd = dar_base.mean(axis=0)
mslpt = tah_base.mean(axis=0)

aDar = dar[months] - mslpd
aTah = tah[months] - mslpt

aDar_base = aDar.loc[base_start:base_end, months].to_numpy().ravel()
aTah_base = aTah.loc[base_start:base_end, months].to_numpy().ravel()

aDar_base = aDar_base[~np.isnan(aDar_base)]
aTah_base = aTah_base[~np.isnan(aTah_base)]

sDar = np.sqrt(np.sum(aDar_base**2) / aDar_base.size)
sTah = np.sqrt(np.sum(aTah_base**2) / aTah_base.size)

nDar = aDar / sDar
nTah = aTah / sTah

diff_base = (nTah.loc[base_start:base_end, months]
             - nDar.loc[base_start:base_end, months]).to_numpy().ravel()

diff_base = diff_base[~np.isnan(diff_base)]

MSD = np.sqrt(np.sum(diff_base**2) / diff_base.size)

SOI = (nTah - nDar) / MSD

records = []
for year in SOI.index:
    for imonth, month in enumerate(months, start=1):
        records.append({
            "YEAR": year,
            "MONTH": imonth,
            "SOI": SOI.loc[year, month]
        })

soi_ts = pd.DataFrame(records)
soi_ts["TIME"] = soi_ts["YEAR"] + (soi_ts["MONTH"] - 1) / 12.0
soi_ts = soi_ts.dropna(subset=["SOI"])

plt.figure(figsize=(12, 4))
plt.plot(soi_ts["TIME"], soi_ts["SOI"], linewidth=1.0)
plt.axhline(0, color="k", linewidth=0.8)
plt.xlabel("Year")
plt.ylabel("SOI")
plt.title("Southern Oscillation Index (Tahiti - Darwin)")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("SOI_timeseries.png", dpi=300, bbox_inches="tight")
plt.show()

plt.figure(figsize=(12, 4))
pos = soi_ts["SOI"] >= 0
neg = soi_ts["SOI"] < 0
plt.bar(soi_ts.loc[pos, "TIME"], soi_ts.loc[pos, "SOI"], width=0.07)
plt.bar(soi_ts.loc[neg, "TIME"], soi_ts.loc[neg, "SOI"], width=0.07)
plt.axhline(0, color="k", linewidth=0.8)
plt.xlabel("Year")
plt.ylabel("SOI")
plt.title("Southern Oscillation Index (bar plot)")
plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("SOI_bar.png", dpi=300, bbox_inches="tight")
plt.show()

Python で計算した SOI の時系列図。

正負を色分けした SOI の棒グラフ。