環境情報論　第5回線型トレンドとその除去¶

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

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize

東京の気温を読み込んで気候値と偏差の計算（前回と同じだが，140年分）¶

まず前回と同様のことを，30年ではなく140年についてやってみましょう。

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

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

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

# 偏差の計算
tempa = np.zeros((temp.shape))
for yy in range(1880, 2020):
    for mm in range(1, 13):
        tempa[(y==yy)*(m==mm)] = temp[(y==yy)*(m==mm)] - temp_clim[mm-1]

mon = np.arange(1880, 2020, 1/12)
plt.plot(mon,tempa)
plt.title('Temperature Anomalies in Tokyo')

Text(0.5, 1.0, 'Temperature Anomalies in Tokyo')

前回は使ったデータの期間が短かったのであまり気づかなかったかもしれませんが，実は東京の気温は140年でこんなに上がっています。

このように上がっているのは，地球温暖化と都市化（ヒートアイランド）の両者によると考えられていますが，ここではこれらをまとめて東京の「温暖化」と呼ぶことにしましょう。

気温のデータのように，様々な揺らぎに隠されながら，じわじわと温暖化がどのくらいのペースで上がっているのかを定量的に計算する方法を勉強しましょう。

線型トレンドの定義と計算¶

温暖化のペースを定量化するためにもっとも基礎的な方法は，上記の気温偏差を一次関数で近似することです。これには，（何度か他の授業でも出てきているかもしれませんが）最小二乗法を用います。まず，データの構造がわかりやすいように，先ほどの図を折れ線グラフではなく散布図で書いてみます。

plt.scatter(mon, tempa)

<matplotlib.collections.PathCollection at 0x10d278160>

このように点（$x_i$, $y_i$）($i = 1, 2, 3, ...$)という点列が与えられたとき，その点全体の散らばりを最もよく近似するような直線$y=ax+b$における傾き$a$と切片$b$は，[a, b] = np.polyfit(x, y, 1)という関数で求めることができます。

# 第一引数に横軸，第二引数に縦軸のデータを書く。
# 第三引数の1は，1次関数で近似することを表す。
[a, b] = np.polyfit(mon, tempa, 1)

直線$y=ax+b$を重ねて，赤色で描画してみましょう。

plt.scatter(mon,tempa)
plt.plot(mon, a*mon + b, 'r')

[<matplotlib.lines.Line2D at 0x10d2befd0>]

何をやっているかわかったと思うので，折れ線グラフに戻してみます。

plt.plot(mon,tempa)
plt.plot(mon, a*mon + b, 'r')
plt.title('Temperature Anomalies in Tokyo')

Text(0.5, 1.0, 'Temperature Anomalies in Tokyo')

この気温の上昇をもっともよく説明するような直線$y=ax+b$を線型トレンドといい，特に$a$を回帰係数と呼びます。ここでは横軸が年，縦軸が℃ですので，回帰係数の単位は「℃/年」です。

aの値は

a

0.02505409050955032

ですから，東京の気温は過去140年間で，平均的には1年につき0.025℃ずつ上がってきたことになります（100年で2.5℃）。これは，全地球平均を上回るスピードであるといえます。

線型トレンドの除去（デトレンド）¶

線型トレンドがあるようなデータでは，偏差の大きさが必ずしも「異常さ（=いかに近くの年と違うか）」を反映していません。

たとえば，2003年に明らかに寒かった年の下向きピークが見られますが，偏差の大きさ自体は100年前の「普通の年」とあまり変わりません。

一方，1885年頃の冷夏のピークは，偏差こそ-4℃という負に大きな値ですが，温暖化する前であったことを考えるとその頃にはそこまで「異常」ではなかったかもしれません。

そのような場合，線型トレンドを偏差のデータから差っ引いた方が，その年の「異常さ」を検出するには都合が良いといえます。

このような操作を「線型トレンドの除去」または「デトレンド(detrend)」といいます。

detrended_tempa = tempa - (a*mon + b)
plt.plot(mon,detrended_tempa)

[<matplotlib.lines.Line2D at 0x11e1c6828>]

こうすることで，温暖化シグナルを差っ引いて考えると，2003年の冷夏（デトレンド後の偏差-3.5℃）はやはり「異常」であったことがわかります。

また，1885年ごろの冷夏（デトレンド後の偏差-2.5℃）よりも極端な現象であったことがわかります。

以上のように，大きなトレンドを含むデータを見る時に，トレンド自体には興味がなく，むしろそれよりも短い時間スケールの変動に興味がある場合は，事前にデータをデトレンドしておくと，さまざまな解析がしやすくなります。

pythonでは，わざわざnp.polyfitを経由しなくても，デトレンド自体が関数として与えられています。

from scipy import signal #モジュールのインポート
detrended_tempa_scipy = signal.detrend(tempa) # 「気温偏差をデトレンドせよ」という関数
plt.plot(mon,detrended_tempa_scipy) #上記と同じ結果になる

[<matplotlib.lines.Line2D at 0x1c207db668>]

SST分布のトレンドを描画してみる¶

データセットの読み込み¶

いつも通り海面水温のデータを読み込みます。今日以降しばらく，先週計算しておいた海面水温偏差しか使わないので，海面水温の生データと気候値は読み込む必要がありません。

先週の海面水温偏差データはこちらからもダウンロードできます：http://web.is.ocha.ac.jp/~tsubasa/env_info/ssta_OISST.npz

loadfile = '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に保存

SSTトレンドとデトレンドした偏差の計算¶

早速トレンドを計算していきます。いつも通り，海面水温は2次元場なので，東京の気温の例と同じことを地球上の各点全てについて行います。

ちなみに2次元場の場合は，切片$b$を無視して回帰係数自体をトレンドと呼ぶことが多いです。

# sstのサイズをそれぞれ変数imt, jmt, tmtに保存
[imt, jmt, tmt] = ssta.shape 

# 0で埋められた行列を使って，欲しいサイズの行列を作っておく（初期化）
a_ssta = np.zeros((imt, jmt)) # 回帰係数a
b_ssta = np.zeros((imt, jmt)) # 切片b
detrended_ssta = np.zeros((imt, jmt, tmt)) # デトレンドした海面水温偏差

# 時間軸
mon = np.arange(1982, 2020, 1/12)

# 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(mon, ssta[ii, jj, :], 1) 
        # デトレンド
        detrended_ssta[ii, jj, :] = \
                       ssta[ii, jj, :] - (a_ssta[ii, jj]*mon + b_ssta[ii, jj]) 
    if (ii % 30 == 0):
        print(ii) # ちゃんと計算が進んでるかチェックするために，30回に1回iiを出力する
        
# さっきゼロにしておいた陸地の場所にもう一度nanを戻す
# is_land_grids_2Dは，sstaの1982年1月の値がnanのところだけTrueが入っているような2次元配列（360x180）
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
detrended_ssta[is_land_grids_2D]=np.nan

0
30
60
90
120
150
180
210
240
270
300
330

# デトレンドした偏差は次回使うので保存しておく
savefile = 'detrended_ssta_OISST.npz'
np.savez(savefile, ssta=detrended_ssta, lon2=lon2, lat2=lat2, y=y, m=m)
# detrended_sstaはちょっと長いので，単にsstaという変数名で保存する

SSTトレンドの描画¶

# vminはカラーバーの下限，vmaxはカラーバーの上限
# vintはカラーバーの間隔
vmin = -0.03
vmax = 0.03
vint = 0.005
        
plt.figure()
cm = plt.get_cmap('seismic')
cs = plt.contourf(lon2, lat2, a_ssta, \
                  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 = 'SST trend'
plt.title(title)
plt.xlim(0, 360)
plt.ylim(-90, 90)
plt.show()

過去40年の海面水温については，日本付近を含む大陸東岸（海洋の西岸。黒潮やメキシコ湾流を含む。）で温暖化が速いことがわかります。

逆に，東太平洋や南極海では温暖化が遅かった（むしろ冷たくなった）ことが知られています。

（余談：東太平洋の温暖化が遅い理由は，温室効果ガスのせいなのか自然変動なのかが未解決問題となっています。多くの研究者が自然変動のせいであろうと信じていた中で，神山の博士論文ではこれが温室効果ガスの強制でありうるという仮説を提唱しました。）

デトレンドしたSST偏差の描画¶

# 描画したい年・月
draw_year = 1997
draw_month = 12

# vminはカラーバーの下限，vmaxはカラーバーの上限
# vintはカラーバーの間隔
vmin = -6
vmax = 6
vint = 1
        
plt.figure()
cm = plt.get_cmap('seismic')
cs = plt.contourf(lon2, lat2, \
                  np.squeeze(detrended_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 = 'Detrended SST anomalies ' + str(draw_year) + '/' + str(draw_month)
plt.title(title)
plt.xlim(0, 360)
plt.ylim(-90, 90)
plt.show()

パッと見ただけではあまり違いは見られないかもしれませんが，次回以降にこの操作の効果が見られます。

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

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

A問題. 課題1と同様，あなたの趣味などに関する好きな時系列データ（横軸が時間軸になるようなデータ）について簡単なグラフを描き，ある期間に着目してトレンドを計算してください。単位も付記すること。

B問題. 地上2m気温偏差のデータ http://web.is.ocha.ac.jp/~tsubasa/env_info/t2ma_erai.npz のトレンドを，上記のSSTの例と同じように計算しましょう（.npzファイルの中身の調べ方が分からない方は，第2回を復習してください）。最近の約40年において，北半球と南半球の温暖化はどちらが速かったでしょうか。また，陸上の温暖化と海上の温暖化はどちらが速かったでしょうか。

C問題. C-1. 日経平均株価の月次データ https://indexes.nikkei.co.jp/nkave/historical/nikkei_stock_average_monthly_jp.csv を用いて，2000年から2019年についての日経平均株価の終値についてトレンドを求めてください。単位も付記すること。

ヒント：普通に読むと文字化けする可能性があります（たぶんShift-JIS）。左からデータ日付，終値，始値，高値，安値の順で入っていますので，ヘッダーが読めない場合は削除するなどしてください（あるいは，これを機にShift-JISで読む方法を勉強してみるのも良いかも。pandas等を使っても良いです）。

C-2. C-1.で用いたデータをデトレンドしたデータをグラフに描画し，特に株価が暴落したり高騰したイベントを見つけ，その理由を簡単に考察してください。

D問題. 海氷面積偏差のデータ https://nsidc.org/arcticseaicenews/sea-ice-tools/ を自由に用いて，北極と南極の海氷面積（sea ice area）のトレンドを調べてください。"3. Sea ice extent and area organized by year" というのが良いと思います。北極と南極の海氷はどちらが急激に変化しているでしょうか。

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

環境情報論 第5回 線型トレンドとその除去¶