気象情報解析特論　第6回　主成分分析を用いた気象データの分析（EOF解析）

担当教員: 神山 翼 (tsubasa/at/is.ocha.ac.jp, @t_kohyama, 理学部3号館703号室)

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal
import math
from matplotlib.colors import Normalize
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LatitudeFormatter,LongitudeFormatter
import matplotlib.ticker as mticker

下準備：SSTデータセットの読み込みと関数の定義

今回から3回分の話は，行列を用いたデータ解析の応用です。また海面水温から話を始めます。

In [2]:
loadfile = 'detrended_ssta_OISST.npz' # デトレンド済みの海面水温偏差の入力ファイル名を定義
ssta_dataset = np.load(loadfile) # データセットはまずデータセットごと入力します 
ssta = ssta_dataset['ssta'] # 海面水温偏差を変数sstaに保存
lon2 = ssta_dataset['lon2'] # 経度(longitude)を変数lon2に保存（2は「2次元配列」の意味）
lat2 = ssta_dataset['lat2'] # 緯度(latitude)を変数lat2に保存
y = ssta_dataset['y'] # 年(year)を変数yに保存
m = ssta_dataset['m'] # 月(month)を変数mに保存

# 1月始まり12月終わりになるように，1982年から2019年のみのデータを用いる
ssta = ssta[:, :, (1982 <= y)*(y <= 2019)]
m = m[(1982 <= y)*(y <= 2019)]
y = y[(1982 <= y)*(y <= 2019)] # y自身のサイズが変わってしまうので，yは一番最後に書き換えないとダメ

# 領域平均をとる関数
def aave(west, east, south, north, var = ssta):
    var = var[(west<=lon2[:, 1])*(lon2[:, 1]<=east),:, :]
    var = var[:, (south<=lat2[1, :])*(lat2[1, :]<=north), :]
    aave_var = np.nanmean(np.nanmean(var, 0), 0)
    return aave_var

# 月別の時系列を2つ描画する関数
def plot_2_mon_time(time_series1, time_series2, lower = -3, upper = 3, init_year=1982, fin_year=2019):
    mon = np.arange(init_year, fin_year+1, 1/12)
    plt.plot(mon, time_series1)
    plt.plot(mon, time_series2, 'r')
    plt.plot(mon, 0*time_series1, 'k')
    plt.xlim(init_year, fin_year)
    plt.ylim(lower, upper)
    plt.show()
    
# 気象場を描画する関数
def draw_field(field, fig_title, vmin = -3, vmax = 3, vint = 1):        
    plt.figure()
    cm = plt.get_cmap('seismic')
    cs = plt.contourf(lon2, lat2, field,\
                  cmap=cm, norm=Normalize(vmin=vmin, vmax=vmax),\
                  levels=np.arange(vmin,vmax+vint,vint), extend='both')    
    plt.colorbar(cs)
    plt.xlabel('Longitude')
    plt.ylabel('Latitude')
    title = fig_title
    plt.title(title)
    plt.xlim(0, 360)
    plt.ylim(-90, 90)
    plt.show()

# 回帰図を計算する関数    
def reg_map(ref_time, var):
    # z500のサイズをそれぞれ変数imt, jmt, tmtに保存
    [imt, jmt, tmt] = var.shape 

    # 0で埋められた行列を使って，欲しいサイズの行列を作っておく（初期化）
    a_var = np.zeros((imt, jmt)) # 回帰係数a
    b_var = np.zeros((imt, jmt)) # 切片b
    
    # np.polyfitがエラーを吐かないようにするために，陸地の場所(nanが入っている)に一度ゼロを入れておきたい
    # （nanmeanのようなnanpolyfitという関数はない）
    # is_land_grids_3Dは，sstaの値がnanのところだけTrueが入っているような3次元配列（360x180x456）
    is_land_grids_3D = (np.isnan(var)==True) 
    var[is_land_grids_3D]=0

    # 回帰図の計算
    # 経度方向にimt回，緯度方向にjmt回forループを回す
    for ii in range(0, imt):
        for jj in range(0, jmt):
            [a_var[ii, jj], b_var[ii, jj]] = np.polyfit(ref_time, var[ii, jj, :], 1)

        if (ii % 60 == 0):
            print(ii) # ちゃんと計算が進んでるかチェックするために，60回に1回iiを出力する
    
    # さっきゼロにしておいた陸地の場所にもう一度nanを戻す
    # is_land_grids_2Dは，sstaの1982年1月の値がnanのところだけTrueが入っているような2次元配列（360x180）
    var[is_land_grids_3D]=np.nan
    is_land_grids_2D = np.squeeze(is_land_grids_3D[:, :, 0])
    a_var[is_land_grids_2D]=np.nan
    b_var[is_land_grids_2D]=np.nan
    
    return a_var

主成分分析の復習

環境情報論の第11回と第12回で，主成分分析について学びました。今回はそれをさらに発展させて，lat-lon格子データに主成分分析を用いる方法を学びます。

主成分分析の流れを覚えているでしょうか。前回通り，たとえば$\vec{x1}$を東京の気温偏差の時系列，$\vec{x2}$を宇都宮の気温偏差の時系列とします。

I. データ行列の定義

$$X := \left( \begin{array}{c} \vec{x_1}\\ \vec{x_2} \end{array} \right) = \left( \begin{array}{cccc} x_{1, 1} & x_{1, 2} & \cdots & x_{1, N}\\ x_{2, 1} & x_{2, 2} & \cdots & x_{2, N} \end{array} \right)$$

II. 共分散行列の計算

$$C := \left( \begin{array}{cc} \mathrm{cov}(\vec{x_1}, \vec{x_1}) & \mathrm{cov}(\vec{x_1}, \vec{x_2})\\ \mathrm{cov}(\vec{x_2}, \vec{x_1}) & \mathrm{cov}(\vec{x_2}, \vec{x_2}) \end{array} \right) = XX^T/N$$

III. 共分散行列の固有値固有ベクトルを求める

$$CE = E \Lambda $$

IV. 固有値の大きい固有ベクトルから順に，第1主成分，第2主成分，...の方向とする（このとき，第1主成分は分散最大方向となり，第2主成分以降は「固有ベクトルが直交する」という制約の中で分散の大きい順となる。）

$$ E = \left( \begin{array}{cc} \vec{e_1} & \vec{e_2} \end{array} \right)$$$$ \Lambda = \left( \begin{array}{cc} \lambda_1 & 0\\ 0 & \lambda_2 \end{array} \right), \ \ \lambda_1 \geq \lambda_2 > 0$$

V. 固有ベクトルを並べた行列に，元の行列を射影すると，PC1時系列，PC2時系列，...を得る（このとき，それぞれのPC時系列は無相関であり，固有値は各PC時系列の分散を与える。主成分分析前後で，分散の和は保存される。）

$$ Z = E^{T}X \\ (\iff \vec{PC1} = \vec{e_1}^T X, \vec{PC2} = \vec{e_2}^T X)$$

ただし

$$Z := \left( \begin{array}{c} \vec{PC1}\\ \vec{PC2} \end{array} \right)$$

ここまでで，頭の中にはてなマークがいっぱい浮かんでいる方は，まず環境情報論の第11回と第12回を復習してください（笑）。

http://web.is.ocha.ac.jp/~tsubasa/env_info/Lecture11.html

http://web.is.ocha.ac.jp/~tsubasa/env_info/Lecture12.html

海面水温の主成分分析

さて今回は，いよいよ主成分分析をlat-lon格子データに適用します。主成分分析を行うときは，必ず事前に平均値（気候値）とトレンドが除去された偏差データとなっていることを再確認してください。

データ行列の作成

早速データ行列づくりから入りましょう。データ行列は，普通は行方向が時間方向（「サンプリング次元」といいます），列方向が空間方向（「構造次元」といいます）であるように設計します。

ただし，lat-lon格子データの場合は，空間方向だけで2次元分ありますから，まずはこれを1次元になるように，順番に並び替えます。すなわち

SST: lon x lat x time → X: space x time

となるように行列を並び替えるということです。ちゃんと自分で戻せればどんな順でも構いませんが，あまり複雑なルールはやめた方が良いでしょう。

またこのとき，大気海洋のlat-lonデータ特有の注意が二つあります。

注意1: lat-lon格子では，各格子の面積は等しくありません。北極や南極に近いほど，1°x1°の領域の面積は狭くなります。そのため，共分散行列を面積比（緯度のコサイン）で重み付けする必要があります。共分散行列はX同士の掛け算なので，各データに事前に緯度のcosの平方根を乗じておかなければいけません。

注意2: NaNのデータは除外します（たとえばSSTの場合，陸の部分）。いまはSSTデータを扱っていますので，あとで戻せるようにどこがNaNだったかをlsmask（Land Sea mask）という変数に保存しておきましょう。

In [3]:
# latlon格子のデータをデータ行列Xに変換する関数
## var: 主成分分析したい3次元配列(lon, lat, time)
## lat2: 各格子の緯度が入った2次元行列
## mask: NaNではないところに1,NaNのところは0が入っている2次元行列
def latlon_to_X(var, lat2, mask):
    
    # 行列のサイズを得る
    imt, jmt, tmt = var.shape
    
    # 出力するデータ行列を初期化しておく（このとき，最大サイズ（すべて海）を仮定）
    X = np.zeros((imt*jmt, tmt))
    
    # 緯度のコサインの平方根をかけておく
    for tt in range(0, tmt):
        var[:, :, tt] = var[:, :, tt] * np.sqrt(np.cos(lat2*math.pi/180))
    
    # データ行列の作成
    cc = 0 # 何番目のNaNではない格子なのか数えるカウンター
    for ii in range(0, imt):
        for jj in range(0, jmt):
            if mask[ii, jj] == 1 :  # もしii, jj番地の格子がNaNではなければ
                X[cc, :] = np.squeeze(var[ii, jj, :]) # cc行目に，ii, jj番地の経度緯度のSST時系列を入れて
                cc = cc + 1 # カウンターをインクリ
                
    X = np.squeeze(X[0:cc, :]) # すべて海を仮定したが，実際は陸があるので結局cc-1行目までしかデータは入っていない
    
    return X

# NaNではないところに1,NaNのところは0が入っている2次元行列
lsmask = np.isfinite(np.squeeze(ssta[:, :, 0])) 

# SST偏差のデータについて，データ行列Xの計算
X = latlon_to_X(ssta, lat2, lsmask)

出来上がった行列のサイズを見てみましょう。

In [4]:
X.shape
Out[4]:
(44219, 456)

44219行，456列の行列ができました。つまり，海の格子は44218箇所あり，その各場所において456ヶ月分のデータがあるということですね。

主成分分析

「あとはいつも通り共分散行列を作って，固有値固有ベクトルを求めるだけ！」と言いたいところなのですが，大気海洋データの場合は実はもう一つ追加の注意が要ります。

注意3: 大気海洋のlat-lonデータは，大抵は空間方向の方が次元が大きいです。今回のケースだと，44219行，456列のデータから，前回と同じように共分散行列を作ると，44219行44219列の正方行列（！）になります。こんなに大きな正方行列の対角化は，さすがにコンピュータでもかなり時間がかかります。

そこでテクニックとして，データ行列の行と列を入れ替えて（つまり，空間の次元と時間の次元を入れ替えて）主成分分析をします。すると，固有ベクトルがPC時系列を与え，Zが空間構造（普通のやり方で固有ベクトルとして得られたもの。大気海洋分野ではEmpirical Orthogonal Functionの略でEOFと呼びます）を与えます。 こうすることで，456行456列の正方行列を対角化すれば済みます（それでも結構大変そうに見えますが，numpyならこのくらいは余裕です）。

主成分分析では，サンプリング次元（時間方向）と構造次元（空間方向）を入れ替えても数学的に等価な結果が得られるのです。これは，二重直交性とも整合します。

それではさっそく，固有値展開をしていきましょう。

In [5]:
## 主成分分析（EOFとPC時系列の計算）
def calc_EOF_PC(X):

    ## 共分散行列の計算
    # データ行列の定義
    N = X.shape[0] #データ行列の空間方向の大きさ

    # 注意点3に従って，データ行列を転置してから共分散行列を計算
    X = X.T
    C = X@X.T/N # 最初からX.T@X/Nって書いてもいいよ

    ## 固有値展開！！！
    # 共分散行列の固有値固有ベクトルを求める
    [Lam, E] = np.linalg.eig(C)
    # なぜか中身は実数なのに複素数型にキャストされる（未解決...）ので，実数型に直す
    Lam = np.real(Lam) # 固有値の入ったベクトル
    E = np.real(E) # 固有ベクトルが入った行列
    
    ## 後処理
    # 固有値の大きい順にソート
    index = Lam.argsort()[::-1]   
    Lam = Lam[index]
    E = E[:, index] # PC時系列を並べた行列
    # 対角成分に固有値を並べる
    L = np.diag(Lam) # 対角成分に大きい順に固有値の入った行列
    # 元のデータ行列を固有ベクトルに射影
    Z = E.T@X # 空間構造(EOF)を並べた行列
    
    return E, L, Z

E, L, Z = calc_EOF_PC(X)

EOFの描画

それでは，各EOFを描画するために，Z行列をlat-lon格子に戻しましょう。つまり，

すなわち

Z: mode x space → EOFs: lon x lat x mode

となるように行列を並び替えます。このとき，先ほどデータ行列Xを作ったときのルールを逆にたどります。

In [6]:
# 展開行列Zをlatlon格子のデータに変換する関数
## Z: 主成分の空間構造の入った展開行列Z(mode, space)
## mask: NaNではないところに1,NaNのところは0が入っている2次元行列
def Z_to_latlon(Z, mask):
        
    # 行列のサイズを得る
    mmt, smt = Z.shape
    imt, jmt = mask.shape
    
    # 出力するデータ行列を初期化しておく（このとき，最大サイズ（すべて海）を仮定）
    EOFs = np.zeros((imt, jmt, mmt))
    
    # lat-lon行列の作成
    cc = 0 # 何番目のNaNではない格子なのか数えるカウンター
    for ii in range(0, imt):
        for jj in range(0, jmt):
            if mask[ii, jj] == 1 :  # もしii, jj番地の格子がNaNではなければ
                EOFs[ii, jj, :] = Z[:, cc] # ii, jj番地に，cc列目の格子のデータを入れる
                cc = cc + 1 # カウンターをインクリ
            else:  # もしii, jj番地の格子がNaNであれば
                EOFs[ii, jj, :] = np.nan # EOFsの方のii, jj番地にもNaNを入れる
                 
    return EOFs

sst_EOFs = Z_to_latlon(Z, lsmask)

それでは，第1主成分のEOF（=EOF1といいます）を描画してみましょう！

In [15]:
eof1 = sst_EOFs[:, :, 0];
eof1 = eof1/np.sqrt(np.nansum(eof1*eof1)) # EOFは単位ベクトルに直すのが通例
draw_field(eof1, 'EOF1 of global SST', -0.04, 0.04, 0.004)

SSTの第一主成分は，エルニーニョ南方振動(ENSO)でした！

Niño3.4指数などは，太平洋の変動を見たいために主観的に決めたインデックスでした。しかし，なんと今回は，SSTの分散最大方向を探すことによって，客観的にENSOを取り出すことができました。この辺になると，ようやく主成分分析のご利益が分かってくる感じがしますね。

PC時系列

PC1時系列も見てみましょう。これは，Niño3.4指数と似たようなものになっていて欲しいですね。

ちなみに，PC時系列は標準偏差が1になるようにするのが通例です（これを，標準偏差による規格化(normalization) ，または単に 正規化(standardization) といいます）。

In [8]:
# PC1時系列は，最大固有値に対応する固有ベクトルとして得られる
##（符号は逆転させても良いが，その場合はEOF1の符号も逆転させる）
pc1 = E[:, 0] # 固有ベクトルの抜き出し
pc1 = pc1/np.std(pc1) # 正規化
nino34 = aave(190, 240, -5, 5) # 参考までに，nino34と比べてみる
nino34 = nino34/np.std(nino34) # 振幅を合わせるためにこっちも正規化

plot_2_mon_time(nino34, pc1) # nino34を青で，PC1時系列を赤で描画する

ご覧の通り，ほぼ一致する時系列が得られました！！！

Niño3.4インデックスは，全世界のSSTの第一主成分と一致するからこそ，客観的で良いインデックスの定義だと言えるわけです。

PC時系列への回帰図（振幅の情報を持つEOF）

なお，先ほどEOFは単位ベクトルに直すのが通例と言いましたが，実際に読者に見せる際には，振幅の情報が落ちてしまうのであまり便利なものではありません。そこで論文や学会発表では，EOFを見せずに，PC時系列の回帰図を示すことが多いです。

In [9]:
# SSTのPC1への回帰図を計算
a_pc1 = reg_map(pc1, ssta)
# 描画
draw_field(a_pc1, 'Regression map of SST onto PC1', vmin = -1.5, vmax = 1.5, vint = 0.1)

0
60
120
180
240
300

EOFと全く同じ図が出てきました（回帰もEOFもどちらも線型解析なので，結果はほぼ一致します）！

こうすると，PC1時系列が自身の1標準偏差と等しくなったときに，赤道東太平洋上では約1.5℃程度の振幅があることがわかります。ものすごく大雑把な議論ですが，（仮にPC1に1年で1自由度あり，かつ正規分布に従うと仮定すると，）20年に1度くらいは3℃程度の大きなENSOイベントが発生することになります。

つまりPC1への回帰図を見せることで，EOF1の構造を保ったまま，振幅の情報を加えることができるのです。

寄与率の計算

最後に，寄与率を計算しておきましょう。

In [10]:
ev1 = L[0, 0]/np.trace(L) # トレースの計算はnp.traceを用いる

ev1
Out[10]:
0.16438340262378257

つまり，SSTのデータが持つ全分散のうち，16%がENSOによって説明されることになります。「ENSOは地球上で最も卓越するSSTの経年変動モードである」などという言い方をすることもありますが，具体的には，世界のSST経年変動（=数年程度の時間スケールの変動）のうち約1/5弱をENSOだけで説明できてしまうことになります。

EOF解析のまとめ

以上のように，EOF解析を行うことによって得られる第1主成分に関する情報は

・PC1時系列

・EOF1の空間パターン（実際はPC1時系列への回帰図として示す）

・EOF1の寄与率

の3点です。第2主成分以降も必要な場合は，同様に計算していってください。

課題6（締切：次回の授業開始前まで）

この講義資料のようなNotebook形式で次の問に答え，html形式またはpdf形式に保存して提出して下さい。NotebookにはMarkdownセルを用いて，それぞれのセルについて（この講義資料のように）説明を加えて下さい。

A問題. EOF1，PC1時系列，およびその寄与率とは何か，簡潔に説明してください。

B問題. SSTの第2主成分と第3主成分について，それぞれEOFとPCを計算し，どのような現象を示しているか簡単に説明してください（難しく考えず，見たまま書いてくれればそれで良いです）。また，それぞれの寄与率はどの程度でしょうか？

発展的な注：参考文献として，

Hartmann, D. L. (2015), Pacific sea surface temperature and the winter of 2014. Geophys. Res. Lett., 42, 1894– 1902.

があげられます。

C問題. 太平洋十年規模振動(Pacific Decadal Oscillation)，いわゆるPDOの指数は，北太平洋の20°N以北における海面水温偏差のEOF第1モードとして計算されます。PDO指数を計算し，時系列と空間パターンを示してください。PDOとはどのような現象でしょうか？

ヒント：ちなみに正解は https://www.ncdc.noaa.gov/teleconnections/pdo/ のview dataから取得できますので，比べてみてください（使用データや考慮した期間が違うと思うので，厳密に数字がピッタリ合うことはないです。同じ現象が捉えられていそうなら，正解として良いです）。

D問題. 北極振動(Arctic Oscillation），いわゆるAOの指数は，北緯20度以北の海面気圧偏差におけるEOF第1モードとして定義されます。AO指数を計算し，時系列と空間パターンを示してください。北極振動とはどのような現象でしょうか？

ヒント：海面気圧偏差のデータはこちらに置いてあります。 http://web.is.ocha.ac.jp/~tsubasa/env_info/slpa_erai.npz

提出方法：課題提出専用メールアドレスkohyama.coursework/at/gmail.com（普段のメールと違います！！/at/を@に変えてください）に対して，件名を「気象情報解析特論6，〇〇〇〇」として（〇〇〇〇にはご自身の氏名を入れてください），課題ファイルをメールに添付して送信してください。