15.5 ERSSTで見るPDO-like EOF解析

1. 何をしたいのか

ここでは、ERSSTの月平均海面水温を使って、PDOに似た変動パターンを自分で計算してみる。

このPC1をPDO-like indexとして扱い、全球の各グリッドのSST偏差とどのような関係を持つかを調べる。

2. PDOとは何か

PDO（Pacific Decadal Oscillation）は、北太平洋の海面水温偏差に見られる低周波の変動パターンである。正位相では、北太平洋中央部が低温偏差、北米西岸沿いが高温偏差になりやすい。負位相ではおおむね逆の分布になる。

ただし、正位相と負位相は完全な赤青反転とは限らない。特に、南半球や熱帯ではENSO、IPO、全球的な昇温、解析期間の偏りなどが混ざるため、コンポジットでは非対称性が出る。

3. 解析の流れ

1. ERSSTを読み込む。
2. 月別気候値を引いて、SST anomalyを作る。
3. 必要に応じて、全球平均SST anomalyを各時刻から引く。
4. 北太平洋 20–70°N, 120–260°E を切り出す。
5. 高緯度ほど格子面積が小さいため、sqrt(cos(lat)) の面積重みをかける。
6. 時間 × 空間グリッドの2次元行列に変換する。
7. SVDでEOF/PCAを行い、第1モードを取り出す。
8. PC1を標準化してPDO-like indexとする。
9. PC1と全球SST anomalyの相関係数を計算する。
10. PC1が強く正・負の月を集め、SST anomalyのコンポジットを作る。

4. PDO-like PC1時系列

第1主成分 PC1 を標準化したものが、ここでのPDO-like indexである。正の値を赤、負の値を青で塗り、黒線は平滑化した時系列である。

図1：ERSSTから計算したPDO-like PC1時系列。正の月は赤、負の月は青で表示した。

正の期間と負の期間が数年〜十年程度で入れ替わっているように見える。これは、PDOが単なる月々のノイズではなく、比較的ゆっくり変動する成分を含むことを示している。

5. EOF第1モードの空間パターン

次に、北太平洋だけで計算したEOF第1モードの空間パターンを見る。中央北太平洋と北米西岸側で符号が反対になっており、PDO-likeな構造が確認できる。

図2：北太平洋SST anomalyに対するEOF第1モード。寄与率は約22%。

6. PC1と全球SSTの相関係数マップ

PDO-like PC1と、全球の各グリッドのSST anomalyとの相関係数を計算した。これは「PDO-like indexが正のとき、どの海域のSSTも同時に高くなりやすいか、低くなりやすいか」を見る図である。

図3：PDO-like PC1と全球SST anomalyの相関係数マップ。

北太平洋では、EOF第1モードと対応する相関パターンが見える。一方、熱帯太平洋や南半球にも相関が出ている。これはPDO-like PC1が北太平洋だけで完結するわけではなく、太平洋全体の変動やENSO的な変動とも関係している可能性を示している。

7. PDO正位相・負位相コンポジット

相関係数マップとは別に、PDO-like PC1が強く正の月、強く負の月を集めて、それぞれのSST anomaly平均を作った。これは、PDO正位相・負位相の典型的なSST偏差場を見るための図である。

図4：PDO正位相・負位相のSST anomalyコンポジット。半球表示により、NOAA PSLのPDO説明図に近い見せ方にしている。

図5：同じコンポジットを平面図で表示したもの。学生が緯度経度で確認しやすい。

8. 使ったPythonコードの骨格

完全コードは長いが、中心となる処理は次の流れである。

import numpy as np
import xarray as xr

ds = xr.open_dataset("ersst.v5.nc")
sst = ds["sst"].sel(time=slice("1950-01-01", "2021-12-31"))

# 1. 月別偏差
clim = sst.groupby("time.month").mean("time", skipna=True)
ssta = sst.groupby("time.month") - clim

# 2. 全球平均SST anomalyを除去
weights = np.cos(np.deg2rad(ssta["lat"]))
global_mean = ssta.weighted(weights).mean(("lat", "lon"), skipna=True)
ssta = ssta - global_mean

# 3. 北太平洋を切り出す
# ERSSTは緯度が北から南へ降順なので注意
ssta_np = ssta.sel(lat=slice(70, 20), lon=slice(120, 260))

# 4. 面積重み
lat = ssta_np["lat"].values
w = np.sqrt(np.cos(np.deg2rad(lat)))
data_w = ssta_np.values * w[None, :, None]

# 5. time × space に変換
nt, ny, nx = data_w.shape
X = data_w.reshape(nt, ny * nx)
valid = np.all(np.isfinite(X), axis=0)
X = X[:, valid]
X = X - X.mean(axis=0)

# 6. SVDでEOF/PCA
U, S, Vt = np.linalg.svd(X, full_matrices=False)
pc1 = U[:, 0] * S[0]
pdo_index = (pc1 - pc1.mean()) / pc1.std()

9. 解釈の注意点