環境情報論 第9回 回帰図と相関図¶
担当教員: 神山 翼 (tsubasa/at/is.ocha.ac.jp, @t_kohyama, 理学部3号館703号室)
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
下準備：SSTデータセットの読み込みと関数の定義（前回までの復習）¶
今回は，前回学んだ回帰係数と相関係数を，SSTのような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 draw_field(field, fig_title, vmin = -6, vmax = 6, 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 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
# 月別の時系列を描画する関数
def plot_mon_time(time_series, lower = -3, upper = 3, init_year=1982, fin_year=2020):
mon = np.arange(1982, 2020, 1/12)
plt.plot(mon, time_series)
plt.plot(mon, 0*time_series, 'k')
plt.xlim(init_year, fin_year)
plt.ylim(lower, upper)
回帰図¶
まず，前々回の授業で作ったエルニーニョ現象のコンポジット図を見てみましょう。
# Niño3.4指数を計算
nino34 = aave(190, 240, -5, 5)
# 「Niño3.4指数が1℃以上になった月」がTrueになるデータだけを抜き出す
godzilla_nino_data = ssta[:, :, (nino34>1)]
# 抜き出したデータの時間方向（配列の第三次元）について平均をとる
godzilla_nino_composite = np.mean(godzilla_nino_data, 2)
# 描画
draw_field(godzilla_nino_composite, \
'Composited SSTA for months when Niño3.4 > 1 ℃', -3, 3, 0.5)
前々回の復習ですが，この図はNiño3.4が1℃以上の月を取り出して平均を取ることによって作成しました。
plot_mon_time(nino34)
mon = np.arange(1982, 2020, 1/12)
plt.plot(mon,np.ones(nino34.shape), 'r:')
コンポジット解析は，上記のように極端現象の特徴を解析するためには非常に優れた方法である一方で，基準を満たさない大部分の残りのデータは捨てることになります。
エルニーニョ現象のように卓越した現象ならば，比較的S/Nを高くしやすいのでコンポジット解析で良いのですが，全てのデータを使って少しでもS/N比をあげたい場合には，必ずしも良い方法とは言えません。
加えて，インデックスが負になったときの特徴を記述するには，もう一枚図が必要になりますので，最もシンプルに物事を理解したいときには向きません。
このようなコンポジット解析の短所を補うために，コンポジット図と相補的に用いられる方法として，回帰図というものが存在します。
回帰図の作り方¶
回帰図とは，1つのインデックスと各地点の気象場の間で計算される各地点の回帰係数を描画したものです。たとえば，
$x =$ Niño3.4指数 と $y=$「(1°N, 1°E)のSST」との回帰係数を(1°N, 1°E)に描く
$x =$ Niño3.4指数 と $y=$「(1°N, 2°E)のSST」との回帰係数を(1°N, 2°E)に描く
...
$x =$ Niño3.4指数 と $y=$「(90°S, 1°W)のSST」との回帰係数を(90°S, 1°W)に描く
という具合に，地球全体の回帰係数をそれぞれ求めて，全部の回帰係数を一つの地図上に描画します。
数式を使って書くと，$N$か月分のインデックス$x_i$と各点の気象場$y_i(\lambda, \phi)$の組
$(x_i, y_i(\lambda, \phi)) \ \ \ (i = 1, 2, ..., N)$
について，各点の回帰係数$a(\lambda, \phi)$を求めて描画したものといえます。ただしここで，$i$は時刻方向の添字，$N$はデータ数，$\lambda$は経度，$\phi$は緯度です。
# sstのサイズをそれぞれ変数imt, jmt, tmtに保存
[imt, jmt, tmt] = ssta.shape
# 0で埋められた行列を使って，欲しいサイズの行列を作っておく（初期化）
a_ssta = np.zeros((imt, jmt)) # 回帰係数a
b_ssta = 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(ssta)==True)
ssta[is_land_grids_3D]=0
# 回帰図の計算
# 経度方向にimt（=360）回，緯度方向にjmt（=180）回forループを回す
for ii in range(0, imt):
for jj in range(0, jmt):
[a_ssta[ii, jj], b_ssta[ii, jj]] = np.polyfit(nino34, ssta[ii, jj, :], 1)
if (ii % 30 == 0):
print(ii) # ちゃんと計算が進んでるかチェックするために，30回に1回iiを出力する
# さっきゼロにしておいた陸地の場所にもう一度nanを戻す
# is_land_grids_2Dは，sstaの1982年1月の値がnanのところだけTrueが入っているような2次元配列（360x180）
ssta[is_land_grids_3D]=np.nan
is_land_grids_2D = np.squeeze(is_land_grids_3D[:, :, 0])
a_ssta[is_land_grids_2D]=np.nan
b_ssta[is_land_grids_2D]=np.nan
draw_field(a_ssta, 'Regression map of SSTA on Niño 3.4', -1.6, 1.6, 0.2)
この回帰図は，「Niño3.4が1℃上がったとき，各地点の海面水温が何℃くらい上がる傾向があったかを表した図」ということになります。
言い換えると，Niño3.4が1℃下がったときには，赤青が反転した図に近い状態が実現される傾向にあることにもなります。これは，コンポジット図にはなかった特徴です。
発展的な注1：回帰図にはコンポジット図と比べて，
・全てのデータを用いるので，微小な変動成分でもS/N比を上げやすく，検出しやすくなる
・1℃上がったときと1℃下がったときの空間パターンがだいたい鏡写しのときに，図が1枚で済む
というメリットがあります。逆にコンポジットでは，1℃上がったときと1℃下がったときの空間パターンが全く異なるときに，その非対称性を記述できるという点が優れています（たとえばエルニーニョとラニーニャの空間パターンは少し違う）。
発展的な注2：トレンドの描画は，月を表す単調増加のインデックス（mon）に対する回帰図を描いていたということになります。言い換えると，トレンドの描画と同じことを，monだけでなく様々なインデックスに一般化したものが，回帰図の描画であるとも言えます。そして，トレンドに興味がないときにそれを除去すること，つまりmonをインデックスとして計算した線型変動を除去するのがデトレンドでした。これと同様に，mon以外のインデックスに対する回帰から求められる線型変動を除去することもあります。たとえば，ENSOの変動に興味がない解析では，デトレンドと全く同じようにしてENSOの変動を抜き去ることができます。これをENSOの回帰除去（リグレスアウト）といいます。
相関図¶
次に，上記の回帰図の回帰係数について，相関係数で全く同じことをやると，相関図を書くことができます。
# 0で埋められた行列を使って，欲しいサイズの行列を作っておく（初期化）
r_ssta = np.zeros((imt, jmt)) # 相関係数r
ssta[is_land_grids_3D]=0
# 相関図の計算
# 経度方向にimt（=360）回，緯度方向にjmt（=180）回forループを回す
for ii in range(0, imt):
for jj in range(0, jmt):
r_ssta[ii, jj] = np.corrcoef(nino34, ssta[ii, jj, :])[1, 0]
if (ii % 30 == 0):
print(ii) # ちゃんと計算が進んでるかチェックするために，30回に1回iiを出力する
# nanが入っているから変なエラーが出るけど無視でOK
draw_field(r_ssta, 'Correlation map of SSTA on Niño 3.4', -1, 1, 0.2)
この相関図は，「Niño3.4と各地点の海面水温の間にどれだけ比例に近い関係が成り立つかを表した図」ということになります。ここでは，その地点の海面水温がNiño3.4と正の相関を持っていれば赤，負の相関を持っていれば青に塗られています。
回帰図と比べた時に，図がよりのっぺり（？）しているのは，エルニーニョやラニーニャにともなって変動する海面水温の振幅（＝「エルニーニョへの感度」ともいいます）は必ずしも大きくなくても，Niño3.4と高い相関関係がある海域は世界各地に広がっているということを表しています。
課題9（締切：次回の授業開始前まで）¶
この講義資料のようなNotebook形式で次の問に答え，html形式またはpdf形式に保存して提出して下さい。NotebookにはMarkdownセルを用いて，それぞれのセルについて（この講義資料のように）説明を加えて下さい。
A問題. 回帰図・相関図とは何か，自分の言葉で簡単に説明してください。答案は短くて良いです。
B問題. エルニーニョ現象が起こると，南極付近の「ロス海」と呼ばれる海域の海氷が減ることが知られています。ロス海の海氷面積偏差 http://web.is.ocha.ac.jp/~tsubasa/env_info/siaa_RossSea_NSIDC.npz （出典は第5回に記載）をインデックスとして，海面水温への回帰図および相関図を描画してください。熱帯の海にシグナルは見えるでしょうか。
ヒント：海氷面積偏差は1990年から2019年の30年のデータが入っています。回帰図を計算する際，まず海面水温偏差のデータから1990年から2019年を抜き出したのち，その期間について回帰図と相関図を作成してください。海氷面積の数字が大きいので，回帰図のカラーバーはとても小さい数字にしないといけないので注意です。
C問題. 地上気圧偏差 http://web.is.ocha.ac.jp/~tsubasa/env_info/slpa_erai.npz のNiño 3.4に対する回帰図を計算して，海面水温のNiño 3.4に対する回帰図と重ねて表示してください。このとき，地上気圧の回帰図の方はcontour，海面水温の回帰図の方はcontourfを使うようにすると，重ね描きがうまくいきます。熱帯の現象であるエルニーニョ現象が，大気を介して中緯度の気温を変化させるメカニズムが見えるでしょうか。
ヒント：contourの方は，全て黒の線にして，contourfの方はカラーバーをcoolwarmあたりにすると，重ねたときに色が見づらくならないのではないかと思います。contourの等高線間隔も，見やすいようにうまく調整してください。
D問題. 1997年12月および2010年12月の海面水温偏差の空間パターンについて，ENSOをリグレスアウトする前とした後の海面水温偏差データで比較してください。リグレスアウトによって，ENSOに関わる変動成分（=ENSOという現象が存在することによって発生することによって生じる海面水温変動全般のこと）はどの程度除去されているでしょうか。完全に除去されていないとしたら，理由を考えてください。
提出方法：課題提出専用メールアドレスkohyama.coursework/at/gmail.com（普段のメールと違います！！/at/を@に変えてください）に対して，件名を「環境情報論課題9，〇〇〇〇」として（〇〇〇〇にはご自身の氏名を入れてください），課題ファイルをメールに添付して送信してください。