気象情報解析特論　第3回　移動平均

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

import numpy as np
from scipy.special import factorial
import matplotlib.pyplot as plt
import math

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

今回も，Niño3.4指数を用いて説明します。

まず，使用データ https://web.is.ocha.ac.jp/~tsubasa/env_info/detrended_ssta_OISST.npz をダウンロードして，この.ipynbファイルと同じディレクトリに置いてください。

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

# 月別の時系列を描画する関数
def plot_mon_time(time_series, lower = -3, upper = 3, init_year=1982, fin_year=2019):
    mon = np.arange(init_year, fin_year+1, 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ño3.4指数の計算
nino34 = aave(190, 240, -5, 5)

移動平均

今日は，「移動平均」と呼ばれる操作を勉強します。まずはいつも通り，Niño3.4指数のグラフを見てみましょう。

plot_mon_time(nino34)

このグラフには，注目したいエルニーニョとラニーニャの間の振動の他に，細かいギザギザが入っているのがわかるでしょうか。今回は，もう少し長い時間スケールのゆっくりとした変動に興味があるから，細かいギザギザを取り除きたい，というお話（の第一回）です。

ギザギザというのは，言い換えると「数ヶ月以内のバタバタした細かい振動」です。それゆえ，数ヶ月以内におさまった変動は全部なかったことにしたいということで，全ての月のデータを，それぞれの月の前後2ヶ月と合わせた計5ヶ月のデータの平均で置き換える，ということをやります。これを5ヶ月移動平均 (five-month running mean)といいます。

つまり，数式で表すと，全ての月$i$について，

i月の5ヶ月移動平均データ = { (i-2月のデータ) + (i-1月のデータ) + (i月のデータ) + (i+1月のデータ) + (i+2月のデータ) } / 5

と計算するということです（もちろん，$i=1$とかの場合は，前年の11月や12月のデータを使います）。

ただし，この時系列データは1982年に始まって2019年に終わるので，たとえば1982年の1月のデータを計算したくても，1981年の11月のデータがありません。

移動平均を計算するときのデータの始端と終端の処理は，目的によって色々あり得ますが，ここでは

1982年1月の5ヶ月移動平均データ = (1982年1月のデータ + 1982年2月のデータ + 1982年3月のデータ) /3

というように，使える分のデータだけで平均をとることにしましょう。

それでは，Niño3.4指数の5ヶ月移動平均を計算してみましょう。

# 移動平均の関数を定義
def running_mean(time_series, wn):
    b = np.ones(wn)/wn
    time_series_r = np.convolve(time_series, b, mode="same")
    # 愚直にfor文で平均しても良いが，numpyにはconvolveという便利な関数が用意されている

    # 端点を補正する（参考 https://zenn.dev/bluepost/articles/1b7b580ab54e95）
    n_conv = math.ceil(wn/2) #2で割って切り上げ
    time_series_r[0] = time_series_r[0] * wn/n_conv
    for i in range(1, n_conv):
        time_series_r[i] = time_series_r[i] * wn/(i+n_conv)
        time_series_r[-i] = time_series_r[-i] * wn/(i + n_conv - (wn % 2))
         # size%2は奇数偶数での違いに対応するため

    return time_series_r

# Niño3.4の移動平均を計算
nino34_r5 = running_mean(nino34, wn = 5)
# wnはwindowの略。いまは5ヶ月移動平均なのでwn=5とする。
# 必ずしも奇数である必要はないが，偶数にすると真ん中が半月ずれる

これで移動平均の時系列が計算できました。比較のために，まずは元の時系列をもう一度描画しておきます。

plot_mon_time(nino34)

次に，これが移動平均後の時系列です。

plot_mon_time(nino34_r5)

なんだかさっきよりもツルツルになってますね。笑 これで，熱帯太平洋のゆっくりとした変動のみに注目できます。移動平均成功です！

次に，もっと大きなwindowで移動平均を取るとどうなるでしょうか？

# 15ヶ月移動平均
nino34_r15 = running_mean(nino34, wn = 15)
plot_mon_time(nino34_r15)

フニャフニャです。そして，1982年，1997年，2015年に3回あった大きなエルニーニョが目立たなくなっています。

2つ同時に描いてみましょう。

# 月別の時系列を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()
    
plot_2_mon_time(nino34, nino34_r15)

たしかに，細かいギザギザが消えて，ゆったりとした変動のみが取り出されていることがわかると思います。

一般に，windowが大きくなればなるほど細かいシグナルは消されて滑らかな曲線になりますが，極端現象には注目しづらくなります。

移動平均のもう一つの目的：周期的な変動を除去する

先ほど，移動平均の目的は「ゆっくりとした変動のみに注目する」ことだと述べました。

しかし，実はもう一つ移動平均が使える場面があります。それは，周期的な変動を除去できるということです。

また，東京の気温を例に説明してみましょう。

tokyo_temp = np.genfromtxt("Tokyo_temp.csv",  # ファイルのパスを書く
                  delimiter=",",    # 区切り文字
                  usecols=(0, 1, 2) # 読み込みたい列番号
                 )
y = tokyo_temp[:, 0]
m = tokyo_temp[:, 1]
tokyo = tokyo_temp[:, 2]

# 今回は1990年から2019年の30年分のデータを用いる
tokyo = tokyo[(1990 <= y)*(y <= 2019)]
m = m[(1990 <= y)*(y <= 2019)]
y = y[(1990 <= y)*(y <= 2019)]

まず，季節変動を引く前の気温を描画してみます。

plot_mon_time(tokyo, lower = 0, upper = 30, init_year=1990, fin_year=2019)

相変わらず，季節変動がバタバタしていますね。これでは，どの年が猛暑・暖冬でどの年が冷夏・寒冬だったのかわかりません。

そこで，いつもなら気候値を引いて偏差を出すところですが，今回はここで12ヶ月移動平均をとってみましょう。

（ただしこのとき，季節変動がとても大きいことによって端点では嫌な動きをすることが予想されるので，1991年から2018年のみ描画することにします。）

tokyo_r = running_mean(tokyo, wn = 12)
plot_mon_time(tokyo_r[12:-12], lower = 15, upper = 18, init_year=1991, fin_year=2018)

このように，12ヶ月移動平均だと，ちょうど季節一周分のデータが満遍なく入った期間で平均することになるので，季節変動を具合良く取り除くことができます。

試しに，気候値を引いて偏差を引いた場合の12ヶ月移動平均と比べてみましょう。

# 気候値の計算
tokyoc= np.zeros((12))
for mm in range(1, 13): 
    tokyoc[mm-1] = np.nanmean(tokyo[m==mm], 0)

# 偏差の計算
tokyoa = np.zeros((tokyo.shape))
for yy in range(1990, 2020):
    for mm in range(1, 13):
        tokyoa[(y==yy)*(m==mm)] = tokyo[(y==yy)*(m==mm)] - tokyoc[mm-1]
        
tokyoa_r = running_mean(tokyoa, wn = 12)
plot_mon_time(tokyoa_r[12:-12], lower = -1.5, upper = 1.5, init_year=1991, fin_year=2018)

さっきとほとんど同じ時系列になりました！ただ，偏差の場合は，平均が0になっていますね。このように，データ全体の平均が消えてしまうのが面倒なときは，移動平均で周期変動を除去するとシンプル ということを覚えておいてください。

たとえば，新型コロナウイルスの感染者数は，月曜日が一番少なくて木曜日が一番多いことが知られています。その場合は，使用できるデータ全体でまず曜日ごとに平均をとりそれを差し引いた偏差を用いても良いですが，それだとデータ全体の平均からの外れ具合がわかるだけで，肝心の感染者数自体の大きさの情報がわからなくなるので，データ全体の平均をもう一度足さなければいけなくなってしまいます。そんな面倒なことをするよりは曜日の変動があるデータは7日間移動平均を取りましょうというのがシンプルですね。

逆に，主成分分析を行う時のように，偏差がゼロになる方が都合が良いこともあります 。その場合は，今まで習ってきた通りに区間平均（気候値）を引きましょう。

最後に，せっかくなので東京の気温でも，移動平均前後の偏差を並べて描画してみましょう。

plot_2_mon_time(tokyoa, tokyoa_r, init_year=1990, fin_year=2019)

やはり移動平均をするとゆったりとした変動のみがうまく取り出されていて，「どの年が暑かった（寒かった）」などを議論するには都合が良さそうです。

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

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

A問題. 移動平均の計算方法と，その目的（2つ）について述べてください。

B問題. 太平洋十年規模振動（Pacific Decadal Oscillation, PDO）のインデックスを https://www.ncdc.noaa.gov/teleconnections/pdo/ のview dataから取得し，移動平均を取る前のデータと5年移動平均を取ったデータを比べてみてください。十年規模のゆっくりとした振動に注目しやすくなるでしょうか。

C問題. あなたの趣味などに関する好きな時系列データについて，移動平均を計算して考察してください。

D問題. 株式投資で用いられている移動平均の使用方法について，調べてまとめてください。

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