環境情報論　第7回 コンポジット解析

担当教員: 神山 翼 (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データセットの読み込みと関数の定義（前回までの復習）

今日もいきなり海面水温です。まず，先週計算した「デトレンド済みの海面偏差」のデータを読み込みます。

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は一番最後に書き換えないとダメ

次に，気象場を描画する関数を定義します。今回からは，描画する2次元場（fieldと呼ぶことにする）自体を引数にとることで，SST偏差に限らず色々な場を表示できるようにしてみましょう。

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

# 月別の時系列を2つ描画する関数
def plot_2_mon_time(time_series1, time_series2, lower = -3, upper = 3, init_year=1982, fin_year=2020):
    mon = np.arange(1982, 2020, 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()

コンポジット解析（合成図解析）

前回までは，大きなエルニーニョ現象が起こったときの例として，以下のようにある特定の月の海面水温偏差をお見せしてきました。

# 1982年12月の海面水温偏差を描画
draw_field(np.squeeze(ssta[:, :, (y==1982)*(m==12)]), 'SSTA 1982/12')

# 1997年12月の海面水温偏差を描画
draw_field(np.squeeze(ssta[:, :, (y==1997)*(m==12)]), 'SSTA 1997/12')

# 2015年12月の海面水温偏差を描画
draw_field(np.squeeze(ssta[:, :, (y==2015)*(m==12)]), 'SSTA 2015/12')

しかし，地球上のどの部分がエルニーニョ現象に関係のある偏差で，どの部分が関係のない偏差なのかは，特定の月の偏差を見るだけではわからないですよね。

そこで，大きなエルニーニョ現象が起こった月を選んで，そのときの海面水温偏差のみの平均をとって描画することによって，エルニーニョ現象に関係のある（=取り出した月に共通する特徴である）海面水温偏差のみを取り出すことができます。

衛星観測が始まった1979年以降，「ゴジラ級」のエルニーニョが起こったのは，上記の1982年12月，1997年12月，2015年12月の3回ですので，それらを平均した海面水温偏差を計算してみましょう。

# 「(1982年または1997年または2015年)かつ(12月)」がTrueになるデータだけを抜き出す
godzilla_nino_data = ssta[:, :, ((y==1982)+(y==1997)+(y==2015))*(m==12)]

# 抜き出したデータの時間方向（配列の第三次元）について平均をとる
godzilla_nino_composite = np.mean(godzilla_nino_data, 2)

# 描画
draw_field(godzilla_nino_composite, 'Composited SSTA for the Godzilla Niño events')

このように平均をとることによって，いま注目したい部分（今の場合は熱帯東太平洋）以外のノイズがかなり抑えられていることがわかります。これは，冒頭で示した個々の月の海面水温偏差の図と比べると一目瞭然です。

以上のように，ある特定の条件を満たすようなデータのみを取り出して偏差の平均をとることにより，その条件を満たすデータに共通する性質を抜き出す解析のことをコンポジット解析（あるいは合成図解析）といいます。コンポジットを行うことによって，上記のように「ゴジラ級のエルニーニョが発生した際の典型的な海面水温偏差」を描くことができます。

客観的なコンポジットの取り方

先ほどは，「ゴジラ級のエルニーニョが起こったのは3回」と言いましたが，僕があなたをデータで騙すために，意図的に嘘をついているかもしれません。

そのようなことがないように，コンポジットを見せるときは，必ずコンポジットをとる基準をはっきりする必要があります。

たとえば，上記のコンポジットは，「Niño3.4が2℃以上になった年の12月のコンポジット」という風に表現できますので，客観的なコンポジットだと言えます。

nino34 = aave(190, 240, -5, 5)
plot_2_mon_time(nino34, 2*np.ones(nino34.shape))
# np.ones(nino34.shape)は，nino34と同じ大きさを持った配列の中身に1を敷き詰めたもの

この図を見ると，確かにゴジラエルニーニョはお見せした3年分なので，僕は嘘をついていなかったことがわかります。

あるいは，ゴジラや12月にこだわらなければ，ちょっとエルニーニョの基準を緩めることもできます。たとえば「Niño3.4指数が1℃以上になった月」のことをエルニーニョ現象だと呼ぶことに決めれば，そのようなコンポジットを作ることもありえます。

# 「Niño3.4指数が1℃以上になった月」なら結構いっぱいある
plot_2_mon_time(nino34, 1*np.ones(nino34.shape))

# 「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)

こうすると，振幅は小さくなったものの，かなり美しい「典型的エルニーニョ像」が導き出されました。

一般に，平均するデータ数が多ければ多いほど，ノイズ（＝データのうち注目していない部分）の少ないコンポジットを得ることができます。

逆に，平均するデータの基準がゆるすぎると，シグナル（＝データのうち注目している部分）がぼやけてしまいます。

シグナル(Signal)とノイズ（Noise）の比をS/N比と呼びますが，S/N比の大きなコンポジットを得るために，コンポジットの客観的基準をしっかり考える必要があります。

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

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

A問題. エルニーニョ現象とは逆に，Niño3.4海域が冷たくなる現象である「ラニーニャ現象」の海面水温偏差コンポジット図を，上の例にならって適切な客観的基準を考えて定めることによって描画してください。ラニーニャ現象の際の海面水温分布は，エルニーニョ現象をただ寒暖反転させたものとはどのように異なっているでしょうか。

B問題. コンポジットをとる基準は，何もコンポジットを取られる側のデータに由来するものである必要はありません。たとえば，地上2 m気温偏差の様子を知るためにコンポジットを取るからといって，地上2 m気温偏差自身から作ったインデックスを用いる必要はありません。

エルニーニョ現象およびラニーニャ現象が発生した際の，地上2 m気温偏差 http://web.is.ocha.ac.jp/~tsubasa/env_info/t2ma_erai.npz のコンポジット図をそれぞれ作成してください（.npzファイルの中身の調べ方が分からない方は，第2回を復習してください）。これにより，エルニーニョ現象とラニーニャ現象が発生した際には，世界はどのような気温の分布になる傾向があるのかを自分の言葉で説明してください。

※両方のデータに共通する期間である1982年から2018年のみ抜き出して解析してくれれば大丈夫です（これを最初にやり忘れるとboolean index did not matchのエラーが出ます）。

ヒント：まずエルニーニョ現象の発生を定義するために，海面水温から計算したNiño3.4インデックスを用いて客観的基準を定めることによって，「エルニーニョ現象が発生した月」のみを地上2 m気温偏差のデータから抜き出します。次に，その抜き出したデータを平均し，気温の空間分布を求めればOKです。抜き出すデータが海面水温ではなくなったこと以外，この資料と手順は全く同じです。それが出来たら，同様のことをラニーニャ現象についてもやってみてください。

C問題 「1982年-2019年の東京の気温偏差（計算方法は第4回を復習）」そのものをインデックスとして，「東京が暑かった/寒かった月」を適切な客観的基準を定めることによって適当な数抜き出し，海面水温データをその抜き出した期間のみについて平均することで，東京が暑かった月や寒かった月の海面水温偏差のコンポジット図を描画することができます。これにより，過去に東京の気温が高く/低くなった際に，どういう海面水温分布が観測される傾向があったかを調べてみてください。

ヒント：今度はB問題とは逆で，インデックスが海面水温由来ではなくなったこと以外，やはりこの資料と手順は全く同じです。余裕があれば，季節を限定したり，カラーバーを変えたりしながら，グローバルに，あるいは日本近海のみに着目しても構いません。自由に考察してください。

D問題. Brier and Bradley (1964) https://journals.ametsoc.org/jas/article/21/4/386/17078 という論文のPDFを読んで，月の満ち欠けと雨量の関係を調べるために，コンポジット解析がどのように使われているのかを自分の言葉で説明しましょう。

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