気象情報解析特論　第5回　バンドパスフィルタと気象の時間スケール

担当教員: 神山 翼 (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
from matplotlib.colors import Normalize
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LatitudeFormatter,LongitudeFormatter
import matplotlib.ticker as mticker

下準備：東京の気温

今回の話は，まず東京の気温から始めます。

東京の気温データ：http://web.is.ocha.ac.jp/~tsubasa/env_info/Tokyo_temp.csv （出典は環境情報論第3回に記述）

In [2]:
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)]

# 気候値の計算
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 = signal.detrend(tokyoa)

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

バンドパスフィルタ

前回，「ローパスフィルタ」「ハイパスフィルタ」という概念を勉強しました。それぞれ，低周波領域と高周波領域のみを通すフィルタで，ゆったりとした変動を取り出したい，バタバタした変動を取り出したいときに使うものでした。

それぞれの応答関数は，次のようになりました。

In [3]:
# サンプリング間隔
delta_t = 1

# 青で1次ローパスバタワースフィルタの応答関数を描画
b, a = signal.butter(1, 0.08, btype='low', fs = 1/delta_t)
w, h = signal.freqz(b, a)
plt.plot(w*1/(2*np.pi*delta_t), np.real(h), 'b', label="1st-Order Low-pass Butterworth Filter") 

# 赤で1次ハイパスバタワースフィルタの応答関数を描画
b, a = signal.butter(1, 0.08, btype='high', fs = 1/delta_t)
w, h = signal.freqz(b, a)
plt.plot(w*1/(2*np.pi*delta_t), np.real(h), 'r', label="1st-Order High-pass Butterworth Filter") 

plt.legend()
plt.xlabel("Frequency (/month)") 
plt.ylabel("Response Function")
plt.show()

それでは，低周波でもなく，高周波でもなく，ある特定の周波数帯のみを通したい場合はどのようにすれば良いでしょうか。

そこで登場するのが，バンドパスフィルタです。バンドパスフィルタは，カットオフ周波数をローパス側とハイパス側で二つ持つようなフィルタです。

In [4]:
# 引数は，時系列，サンプリング間隔（月別データだと1ヶ月），カットオフ周波数（2次元配列で与える！），フィルタオーダー
def butterworth_bandpass(time_series, delta_t, f_cut, order):
    b, a = signal.butter(order, f_cut, btype='bandpass', fs = 1/delta_t) # フィルタオーダーは1
    time_series_bp = signal.filtfilt(b, a, time_series) # 与えられた時系列に対してフィルターを実行する
    # 発展的な注：filtfiltは順方向と逆方向の両方でフィルタすることで位相のずれないフィルタにする関数
    return time_series_bp

では，早速東京の気温偏差に，1次のバンドパスバタワースフィルタをかけてみましょう（カットオフ周波数は0.04 /month と0.16 /month）。

In [5]:
tokyoa_bw_bp = butterworth_bandpass(tokyoa, 1, [0.04, 0.16], 1)
plot_2_mon_time(tokyoa, tokyoa_bw_bp)

ローパスとハイパスの中間のような時系列になっているのがわかるでしょうか。これによって，たとえばある特定の周期を持つような現象が生じた際に，そのシグナルを抜き出すことができます。

それでは，応答関数を見てみましょう。

In [6]:
delta_t = 1

# 1次バンドパスバタワースフィルタの応答関数を描画
b, a = signal.butter(1, [0.04, 0.16], btype='bandpass', fs = 1/delta_t)
w, h = signal.freqz(b, a)
plt.plot(w*1/(2*np.pi*delta_t), np.real(h), 'r', label="1st-Order Band-pass Butterworth Filter") 

plt.legend()
plt.xlabel("Frequency (/month)") 
plt.ylabel("Response Function")
plt.show()

ちょうど0.04 /month と0.16 /month，つまり6ヶ月-2年周期くらいの周波数帯のみ通すようなフィルタになっています。

スペクトル解析のまとめ：日本の気象を決定する時間スケール

（参考文献：Dennis L. Hartmann著"Global Physical Climatology" Fig. 8.4）

第一回から，スペクトル解析の基礎を学んできました。その総括として，日本の気象を決定する現象の時間スケールを解析してみましょう。

用いるデータは，ERA-Interimの500 hPaジオポテンシャル高度の日別データ(1979-2011年)で，特に冬に注目します。ジオポテンシャル高度は，ここでは上空の気圧のようなものだと思ってくれれば良いです。

（たとえば飛行機で地上から上に向かって移動すると，気圧が地上付近の1000 hPaからどんどん下がっていくというのを知っていると思います。500 hPaジオポテンシャル高度は，「500 hPaになるような気圧はどの高さか」を表す変数です）

データはこちら：http://web.is.ocha.ac.jp/~tsubasa/env_info/z500_erai_daily_1979_2011.npz

（出典：https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era-interim）

In [7]:
loadfile = 'z500_erai_daily_1979_2011.npz' # 入力ファイル名を定義
sst_dataset = np.load(loadfile) # データセットはまずデータセットごと入力します 
z500 = sst_dataset['z'] # 海面水温(sea surface temperature)を変数sstに保存
lon2 = sst_dataset['lon2'] # 経度(longitude)を変数lon2に保存（2は「2次元配列」の意味）
lat2 = sst_dataset['lat2'] # 緯度(latitude)を変数lat2に保存
y = sst_dataset['y'] # 年(year)を変数yに保存
m = sst_dataset['m'] # 月(month)を変数mに保存
d = sst_dataset['d'] # 日(day)を変数dに保存

今回の趣旨は，日本の冬の気圧変動がどのような現象に支配されているのかを時間スケールごとに分離して見てみたいということです。

このようなときには，まず環境情報論第9回で学習した一点回帰図を描いてみましょう。

In [8]:
# 冬のデータのみ取り出す
z500_DJF = z500[:, :, ((m == 12)+(m == 1)+(m == 2))] #12+1+2月のみ用いる
# 東京付近の点の気圧偏差を抜き出してインデックスとする
z500_Tokyo = z500_DJF[30, 28, :] #東経135度，北緯36度

# varのref_timeに対する回帰図を描く関数
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
    r_var = np.zeros((imt, jmt)) # 相関係数r

    # 回帰図と相関図の計算
    # 経度方向に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 % 20 == 0):
            print(ii) # ちゃんと計算が進んでるかチェックするために，20回に1回iiを出力する
    
    return a_var

a_z500_DJF = reg_map(z500_Tokyo, z500_DJF)

0
20
40
60
In [9]:
# 日本周辺の偏差を描画する関数
def draw_around_Japan(draw_var, vmin = 0.1, vmax = 1, vint = 0.1):

    # 描画する枠を作る
    fig = plt.figure()
    
    # 図の中心経度
    c_lon = 120

    # 枠の中に絵を入れる（図を1枚しか書かないときは1, 1, 1でOK）
    ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree(central_longitude=c_lon))
    
    # 深い青から深い赤に向かうカラーバーを指定
    cm = plt.get_cmap('seismic') 

    # 色で塗られた等高線を描く
    cs = ax.contour(lon2, lat2, draw_var, \
                    np.hstack([np.arange(-vmax, -vmin+0.001, vint), np.arange(vmin, vmax+0.001, vint)]), \
                    cmap=cm, norm=Normalize(vmin=-vmax, vmax=vmax),extend='both', \
                    transform=ccrs.PlateCarree())
                    #colors=['black'], transform=ccrs.PlateCarree())

    # 海岸線を書く（これが気象場を描画するときにめちゃくちゃ重要になる）
    ax.coastlines(lw=0.5,color='gray',resolution='50m')

    # 軸のラベルの間隔（写経でOK）
    dlon,dlat=60,30
    xticks=np.arange(60,360.1,dlon)
    yticks=np.arange(-60,60.1,dlat)
    ax.set_xticks(xticks,crs=ccrs.PlateCarree())
    ax.set_yticks(yticks,crs=ccrs.PlateCarree())

    # 軸のラベルのフォーマット（写経でOK）
    latfmt=LatitudeFormatter()
    lonfmt=LongitudeFormatter(zero_direction_label=True)
    ax.xaxis.set_major_formatter(lonfmt)
    ax.yaxis.set_major_formatter(latfmt)
    ax.axes.tick_params(labelsize=12)
    
    # カラーバーをつける
    plt.colorbar(cs, shrink=0.6) 

    # 描画範囲の指定
    # 正距円筒図法であることを明示するために2つ目の引数が必要
    # 経度の指定では，c_lon±180の値以外を書くと全領域が表示される
    # 0, 360と書くと360が0と同じ数字だと思われて地図が棒になる
    ax.set_extent([0, 240, -30, 90], crs=ccrs.PlateCarree())
    
    plt.show()
    
draw_around_Japan(a_z500_DJF)

これが，何もフィルタリング処理を行わなかった場合の，日本の気圧変動に対する 回帰図（=日本が高気圧のとき，世界の気圧がどうなっている傾向があるか）です。

日本自身が高気圧であることを一番よく説明していることはわかりますが，いまいちどのような現象が日本の気圧を決めているのかはわかりにくいです。

そこで，次に1週間から1ヶ月程度の時定数をもつ変動（7日-30日周期）のみを通すバンドパスフィルタをかけてみましょう。

In [10]:
z500_DJF_bp = butterworth_bandpass(z500_DJF, 1, [1/30, 1/7], 1) 
# 入力データは，時間の時系列が最後の次元になっていればOK
# 発展的注：filtfiltのデフォルトで，多次元配列の場合は最後の次元にfilteringがかかる
z500_Tokyo_bp = z500_DJF_bp[30, 28, :] #東経135度，北緯36度

a_z500_DJF_bp = reg_map(z500_Tokyo_bp, z500_DJF_bp)

0
20
40
60
In [11]:
draw_around_Japan(a_z500_DJF_bp)

今回は，日本と同じくらいの緯度帯に卓越する波のような構造が大事そうだということが見えてきました。

実は，定常ロスビー波という西進する波が，東向きに吹く偏西風の流れに逆らうことで定在することで，この時間スケールの気象を決定しています。

発展的な注：定常ロスビー波は偏西風の強いところに向かって屈折しやすく、ジェット気流が「導波管」となってエネルギーを遠方まで伝播しやすい性質があります。（参考 https://www.jma.go.jp/jma/kishou/books/kisetutext/30/chapter3.pdf ）

もっと短い時間スケールや長い時間スケールの気象を眺めるために，今度はそれぞれハイパスフィルタとローパスフィルタを構成してみましょう。

In [12]:
def butterworth_highpass(time_series, delta_t, f_cut, order):
    b, a = signal.butter(order, f_cut, btype='high', fs = 1/delta_t)
    time_series_lp = signal.filtfilt(b, a, time_series) 
    return time_series_lp

def butterworth_lowpass(time_series, delta_t, f_cut, order):
    b, a = signal.butter(order, f_cut, btype='low', fs = 1/delta_t)
    time_series_lp = signal.filtfilt(b, a, time_series) 
    return time_series_lp

# 7日以内の変動を通すハイパスフィルタ
z500_DJF_hp = butterworth_highpass(z500_DJF, 1, 1/7, 1) 
z500_Tokyo_hp = z500_DJF_hp[30, 28, :] #東経135度，北緯36度

# 30日以上の変動を通すローパスフィルタ
z500_DJF_lp = butterworth_lowpass(z500_DJF, 1, 1/30, 1) 
z500_Tokyo_lp = z500_DJF_lp[30, 28, :] #東経135度，北緯36度

# 回帰図の計算
a_z500_DJF_hp = reg_map(z500_Tokyo_hp, z500_DJF_hp)
a_z500_DJF_lp = reg_map(z500_Tokyo_lp, z500_DJF_lp)

# 回帰図の描画
draw_around_Japan(a_z500_DJF_hp) # 1週間以内の変動
draw_around_Japan(a_z500_DJF_lp) # 1か月以上の変動

0
20
40
60
0
20
40
60

1枚目が1週間以内，2枚目が1ヶ月以上の変動に絞ったときの回帰図です。ご覧の通り，注目する時間スケールによって，気象の決まり方が全然違うことが見えてきました。

まず1週間以内のスケールの日本の冬の気象は，偏西風ジェットによって東向きに流れされてくる短波長のロスビー波の位相によって決定 されます。傾圧不安定というメカニズムで，ロスビー波とともに成長した高気圧と低気圧（移動性擾乱といいます）が西の方から交互にやってきて，日々の天気を決定します。（参考：http://wwwoa.ees.hokudai.ac.jp/people/yamazaki/Lecture/tenki-7.pdf）

それに対して，1ヶ月以上のスケールの日本の冬の気象は，偏西風ジェット自体がどこを通るかで決定 されます。すなわち，偏西風ジェットが南北に移動することで，定在ロスビー波の導波管や，移動性擾乱の「通り道」（これをストームトラックといいます）自体の位置が気圧変動にともなってゆらぐことで，その月や年の天気の「傾向」が決定されるのです。

（偏西風ジェットがどこを通るかは，大気自身の内部変動（北極振動など）のほかに，エルニーニョ現象や境界流同期（Kohyama et al. 2021; 宣伝です）などの海洋現象との結びつきなども勘案しながら決定されます。）

以上のように，フーリエ解析を基にしたスペクトル解析やフィルタリングを行うことによって，混沌の中で何も見えなかったところから，それぞれの時間スケールに切り分けて現象を眺めることができるようになるのです！

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

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

A問題. バンドパスフィルタについて，簡潔に説明してください。

B問題. 日本以外の場所を一つ選んで，上記の資料と同様に一点回帰の解析を行い，時間スケールごとの気象の様相の違いを確かめてください。

C問題. あなたの趣味などに関する好きなデータにフィルタリングを適用し，考察してみましょう。

D問題. 好きな気象データにフィルタリングを施すことで，何か面白いことを考えてください。

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