環境情報論　第6回 インデックスの定義

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

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize

SSTデータセットの読み込み（前回と同じ）

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

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

エルニーニョ発生時（1997年12月）とラニーニャ発生時（2010年12月）のSSTの描画

そろそろ描画に慣れてきた頃だと思うので，ついでに描画の一連の操作を関数にする練習もしてみましょう。

In [3]:
# 関数の定義
## デフォルトの値を引数の横に書いておけば，実行時に省略可能
def draw_ssta(draw_year, draw_month, vmin = -6, vmax = 6, vint = 1, fig_title = 'Detrended SST anomalies '):
        
    plt.figure()
    cm = plt.get_cmap('seismic')
    cs = plt.contourf(lon2, lat2, \
                  np.squeeze(ssta[:, :, (y==draw_year)*(m==draw_month)]), \
                  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 + str(draw_year) + '/' + str(draw_month)
    plt.title(title)
    plt.xlim(0, 360)
    plt.ylim(-90, 90)
    plt.show()
In [4]:
# エルニーニョが発生した年
draw_ssta(1997, 12)
# デフォルトとは違うカラーバーで書きたいときは，たとえば次のように書く
# draw_fig(1997, 12, -3, 3)
In [5]:
# ラニーニャ現象が発生した年
draw_ssta(2010, 12)
In [6]:
# ちょっと弱いエルニーニョが発生した年
draw_ssta(2006, 12)

このようにエルニーニョ状態とラニーニャ状態を行ったり来たりする現象を，エルニーニョ南方振動現象（El Niño Southern Oscillation）あるいはそれを省略してエンソ（ENSO）といいます。ENSOの強弱の時間変化を記述する指数（インデックス）の一つに，Niño3.4指数があるので，今日はまずそれを計算してみましょう。

Niño 3.4指数の計算

上で描画した3つの図を比較するに，熱帯東太平洋のペルー沖の海面水温を見れば，ENSOの変動を記述することができそうです。

そこで，南緯5度 - 北緯5度, 西経170度 - 西経120度で囲まれる領域（Niño3.4領域という）内の海面水温を平均したものを各月ごとに計算して，ENSOを代表する時系列としてプロットします。これがNiño3.4指数です。

In [7]:
#　Niño3.4領域内のデータの取り出し
## まず西経170度 - 西経120度 のデータを取り出す
nino34ssta_data = ssta[(190<=lon2[:, 1])*(lon2[:, 1]<=240),:, :]
## そこからさらに南緯5度 - 北緯5度のデータを取り出す
nino34ssta_data = nino34ssta_data[:, (-5<=lat2[1, :])*(lat2[1, :]<=5), :]
### Pythonでは上記2つを同時にやろうとするとうまくいかない

# データの領域内で平均をとる
## まず東西方向に平均し，次に南北方向に平均する（順番はどちらでも良い）
nino34 = np.nanmean(np.nanmean(nino34ssta_data, 0), 0)
## 今回のケースはnp.meanでも良いが，別のインデックスを定義する場合に
## 領域内に陸地が入ってくることもあるので，最初からnanmeanにしておいた方が安全
## 赤道付近の小さな領域なのでこれで良いが，本当はグリッドの面積による重み付けが必要（D問題参照）

# Niño 3.4指数の描画
mon = np.arange(1982, 2020, 1/12) # 横軸は1982から2020年の1/12年（=1ヶ月）刻み
plt.plot(mon, nino34)
plt.plot(mon, 0*nino34, 'k') # 0のところに黒線を引く
plt.xlim(1982, 2020)
plt.ylim(-3, 3)
plt.show()

Niño3.4指数が正の値のときエルニーニョ状態，負の値のときがラニーニャ状態を表しています。

この指数の絶対値が特に大きくなったとき，それぞれ「エルニーニョ現象が発生した」「ラニーニャ現象が発生した」などといいます。

Niño3.4指数を見ると，たしかに1997年は大きなエルニーニョ現象が発生しています。逆に，たとえば2010年12月には，大きなラニーニャ現象が発生しています。

第4回で言及した「ゴジラエルニーニョ」などと呼ばれる巨大エルニーニョは，1982年，1997年，2015年の3回です。

海面水温の領域平均で定義される色々なインデックス

気候学では，海面水温の領域平均を用いて定義される色々なインデックスがありますので，それを見ていきましょう。

領域平均・描画を行う関数の準備

「領域平均をとって描画する」という操作をいっぱいするので，まず関数に定義しておきます。

In [8]:
# 領域平均をとる関数
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.figure;
    plt.plot(mon, time_series)
    plt.plot(mon, 0*time_series, 'k')
    plt.xlim(init_year, fin_year)
    plt.ylim(lower, upper)
    plt.show()

Niño 3.4指数（上と同じ）

関数を定義しておけば，先ほど計算したNiño3.4インデックスは2行で書けます！

In [9]:
nino34 = aave(190, 240, -5, 5)
plot_mon_time(nino34)

Niño 3指数

ENSOを定義するインデックスには，Niño3指数というインデックスもあります。南緯5度 - 北緯5度, 西経150度 - 西経90度で囲まれる領域（Niño3領域という）内の海面水温を平均したものです。

日本の気象庁では，ENSOの記述にはこちらの指標を用いています（国ごとにエルニーニョの定義に流派のようなものがあります笑）。Niño3の方が定義領域が東に寄っていて，エルニーニョとラニーニャの非対称性（エルニーニョ現象の方がラニーニャ現象よりも強くなりやすい）を表現できるインデックスです。

In [10]:
nino3 = aave(210, 270, -5, 5)
plot_mon_time(nino3, -4, 4)

エルニーニョもどき指数

エルニーニョ現象に似て非なる現象として，エルニーニョもどき (El Niño Modoki)と呼ばれる現象が知られています。エルニーニョ現象では熱帯東太平洋の海面水温が暖まるのに対して，エルニーニョもどき現象では熱帯中央太平洋の海面水温が暖まります。

エルニーニョもどきを定義するエルニーニョもどきインデックス（EMI）は，次の式で計算されます。

EMI = SSTA$_{A}$ - 0.5 SSTA$_{B}$ - 0.5 SSTA$_{C}$

ただし，右辺の３つの項にあるSSTA${X}$は左から順に，領域A(165°E-140°W, 10°S-10°N)，B(110°W-70°W, 15°S-5°N)，C(125°E-145°E, 10°S-20°N)でそれぞれの領域平均した海面水温偏差を表します。ここで，°Eは東経，°Wは西経，°Sは南緯，°Nは北緯を表します。

In [11]:
emi = aave(165, 220, -10, 10) - 0.5*aave(250, 290, -15, 5) - aave(125, 145, -10, 20)
plot_mon_time(emi, -2, 2)

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

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

A問題. 株価の変動を表すインデックス（日経225, TOPIXなど）を一つ選んで，どのように計算されているかを調べて簡潔に記述してください（答案の長さは問わない）。

B問題. ダイポールモード指数について調べて，次の問に答えてください。

  1. 1982年-2019年のダイポールモード指数を計算して時系列を描画してください。

  2. ダイポールモード指数が正に大きかった月と負に大きかった月を1ヶ月ずつ取り出し，上記の最初の例と同様に（インド洋周辺の）海面水温偏差を描画することによって，インド洋ダイポールモード現象を特徴づける海面水温の分布について分かることを簡潔に記述してください。

C問題. あなたの趣味などの好きなデータについて，何か面白いインデックスを定義し，考察を自由に語ってください。

D問題. 上記の領域平均のプログラムには，少しだけ近似があります。赤道から離れたときに，より正確に領域平均を計算するには，各点が代表する面積による「重み付き平均」を計算する必要があります（緯度が大きい方が，経度1度×緯度1度で囲まれる面積が小さい）。重み付けを含んだ領域平均を計算するためのaave_weighted関数を実装してください。aave関数とaave_weighted関数では，結果にどの程度の違いが出るのか，適当な指数を計算して確かめてみましょう。

ヒント：緯度を$\phi$とすると，単位球面のヤコビアンは$\cos\phi$ですので，$\cos\phi$で重み付けをすれば良いことになります。

$\cos\phi$をマクローリン展開すると

$\cos \phi = 1 - \phi^2/2 + O(\phi^4) $

ですので，緯度$\phi$が十分小さい赤道付近の十分小さい領域で平均をとる分には，重み付けをしなくても結果はあまり変わらないということになります。

発展的な注：普通の極座標だとヤコビアンは$R^2\sin\theta$ですが，今回は平均の相対的な重み付けをするだけなので，地球半径を考える必要がないことに加え，緯度の$\phi$は通常の極座標の$\theta$で書くと$\phi = \pi/2 - \theta$ですので，ヤコビアンは$\cos\phi$でいいことになります。

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