02 南方振動指数の計算

1. 背景：ENSO と南方振動指数

エルニーニョとラニーニャは、海洋と大気が結びついた変動であり、まとめて ENSO（El Niño-Southern Oscillation）と呼ばれる。エルニーニョ時には赤道太平洋東部で暖水が広がり、大気循環も平年と異なる。

その大気側の変動を表す代表的な指標のひとつが 南方振動指数（SOI） である。SOI は、タヒチとダーウィンの海面気圧の変動を用いて計算される。

エルニーニョ時（上）は赤道太平洋東部で暖水が広がり、ラニーニャ時（下）は逆に冷水が強まる。これに対応して大気循環も変化する。

2. 使用データ

本演習では、タヒチおよびダーウィンの月平均海面気圧データを用いる。ファイルはこのページと同じディレクトリに置いてある。

3. 計算の流れ

1. タヒチ・ダーウィンの海面気圧データを読み込む
2. 1981–2010 年を基準期間として climatology（各月平均）を求める
3. anomaly（平均からのずれ）を計算する
4. 標準偏差を用いて標準化する
5. MSD（標準化した差のばらつき）を求める
6. SOI を計算する
7. 時系列として描画する

4. Pythonでよく出る機能の意味

SOI の計算では、計算そのものより前に、データの形を整える操作がたくさん出てくる。ここを飛ばすと、コードを読んでいて急に難しく感じる。

5. データの読み込み

まずは Python でテキストファイルを読み込む。ここでは pandas を用いる。

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

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

print(dar.head())
print(tah.head())
print(dar.columns)
print(dar.shape)
print(dar.tail())

ポイント

• sep は区切り文字の指定
• 今回は空白の数が一定でないので r"\s+" を使う
• na_values=[-99.9] により、欠損値を NaN として読む
• engine="python" は柔軟な読み込みを使う指定

6. climatology の計算

基準期間 1981–2010 年の各月平均を求める。これが climatology である。

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

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

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

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

print(mslpd)
print(mslpt)
print(len(mslpd))
print(mslpd["JAN"], mslpd["AUG"])

ポイント

• set_index("YEAR") により、年をインデックスとして使う
• .loc[行, 列] はラベルで範囲を指定して取り出す
• axis=0 は縦方向、つまり各月ごとの平均を意味する

7. anomaly の計算

観測値から climatology を引いて anomaly（平均からのずれ）を求める。

aDar = dar[months] - ______
aTah = tah[months] - ______

print(aDar.head())
print(aTah.head())
print(aDar.loc[1981:2010].mean(axis=0))
print(aTah.loc[1981:2010].mean(axis=0))

ポイント

• anomaly = 観測値 − climatology
• 単位はまだ hPa のまま
• 基準期間 1981–2010 の平均はおおむね 0 に近くなる

8. 標準化の計算

anomaly を標準偏差で割り、スケールをそろえる。

aDar_base = aDar.loc[1981:2010, months].to_numpy().ravel()
aTah_base = aTah.loc[1981:2010, 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) / ______)
sTah = np.sqrt(np.sum(aTah_base**2) / ______)

nDar = aDar / ______
nTah = aTah / ______

print(sDar, sTah)
print(nDar.head())
print(np.nanmean(nDar.loc[1981:2010, months].to_numpy()))
print(np.nanmean(nTah.loc[1981:2010, months].to_numpy()))

ポイント

• to_numpy().ravel() により、30年×12ヶ月の表を 1 次元の配列にする
• 全データをまとめて 1 つの標準偏差を出したいので 1 列に並べる
• ~np.isnan(...) は NaN でない値だけを残す処理
• anomaly はすでに平均を引いた値なので、標準偏差は sqrt(sum(x^2)/N) で求められる

9. MSD の計算

標準化したタヒチとダーウィンの差のばらつきを MSD として求める。

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

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

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

ポイント

• MSD は 標準化した差 のばらつき
• SOI は単なる差ではなく、差をさらに標準化した指標である

10. SOI の計算

SOI = (______ - ______) / ______
print(SOI.head())

SOI は、標準化したタヒチとダーウィンの差を MSD で割って求める。

11. 描画とデータ保存

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"])
plt.axhline(0, color="k")
plt.xlabel("Year")
plt.ylabel("SOI")
plt.title("Southern Oscillation Index")
plt.show()

# 月別テーブル型
SOI_out = SOI.copy()
SOI_out.insert(0, "YEAR", SOI_out.index)
SOI_out.to_csv("SOI_monthly_table_updated.txt",
               sep="\t", index=False, float_format="%.4f")

# 1次元時系列型
soi_ts.to_csv("SOI_timeseries_updated.txt",
              sep="\t", index=False, float_format="%.4f")

ポイント

• 2次元表の SOI を 1 次元の時系列に並べ直す
• 描画して、ENSO の時間変動を視覚的に確認する
• 計算した SOI を txt ファイルとして保存し、次回の課題で利用する
• to_csv(...) を使うと、表をテキストファイルとして保存できる
• sep="\t" はタブ区切りを意味する
• index=False は左端の行番号を保存しない指定である

12. 解答ページ

完全スクリプトは別ページにしてある。まず自分で考えてから必要なときだけ参照すること。

13. 提出方法（重要）

本課題は、Jupyter Lab で作成した Notebook を HTML として保存し、Teams に提出すること。

提出手順

1. Jupyter Lab 上で Notebook を開く
2. 上部メニューから File → Save and Export Notebook As → HTML を選択する
3. HTML ファイルを保存する
4. 保存した HTML ファイルを Teams に提出する

Jupyter Lab のメニューから HTML として保存する。