環境情報論　第4回 偏差の計算

担当教員: 神山 翼 (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

東京の気温を読み込む（先週の復習）

こちらからデータをダウンロードしてください（出典は第3回C問題に記載）：http://web.is.ocha.ac.jp/~tsubasa/env_info/Tokyo_temp.csv

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

簡単のために，今回は1990年から2019年の30年間のデータのみを用いて考えることにしましょう。

In [3]:
temp = temp[(1990 <= y)*(y <= 2019)]
m = m[(1990 <= y)*(y <= 2019)]
y = y[(1990 <= y)*(y <= 2019)] # y自身のサイズが変わってしまうので，yは一番最後に書き換えないとダメ

気候値を計算しておく（これも先週の復習）

In [4]:
month = np.arange(1, 13, 1)
temp_clim = np.zeros((12))
for mm in range(1, 13): 
    temp_clim[mm-1] = np.nanmean(temp[m==mm], 0)
plt.plot(month, temp_clim)
plt.xticks(month)
plt.show()

偏差の定義と計算

まず，360か月分の月平均気温をそのまま描画してみましょう。

In [5]:
mon = np.arange(1990, 2020, 1/12) # 横軸は1990から2020年の1/12年（=1ヶ月）刻み
plt.plot(mon, temp)
Out[5]:
[<matplotlib.lines.Line2D at 0x116c3c810>]

このままでは，季節変動が大きすぎて，どの年が暑くてどの年が寒かったのかがわかりづらいですよね。

そこで，気候値からの差として偏差(anomaly)を定義します。

1990年1月気温偏差　= 1990年1月気温 - 1月の気候値
1990年2月気温偏差　= 1990年2月気温 - 2月の気候値
...
1990年12月気温偏差　= 1990年12月気温 - 12月の気候値
1991年1月気温偏差　= 1991年1月気温 - 1月の気候値
1991年2月気温偏差　= 1991年2月気温 - 2月の気候値
...
2020年9月気温偏差　= 2020年9月気温 - 9月の気候値

偏差とは，「平年からのずれ」のことです。

たとえば東京の気温で，偏差が正の場合は平年よりも暖かく，偏差が負の場合は平年より寒かったということを表します。

それでは，早速偏差を計算していきます。

In [6]:
# 配列のサイズ
tmt = temp.shape 

# 0で埋められた行列を使って，欲しいサイズの行列を作っておく（初期化）
tempa = np.zeros((tmt))
# tempaはtemperature anomalyの略。天パではない。

# 偏差の計算
## 1990年1月から順番にfor文を回していく
for yy in range(1990, 2020):
    for mm in range(1, 13):
        tempa[(y==yy)*(m==mm)] = temp[(y==yy)*(m==mm)] - temp_clim[mm-1]

plt.plot(mon,tempa)
Out[6]:
[<matplotlib.lines.Line2D at 0x116d7b510>]

これが東京の偏差です。0よりも大きい数字が，平年（気候値）よりも暖かい月，0よりも偏差が小さい月が，平年（気候値）よりも寒かった月です。

偏差を計算することで，気温をそのまま描画するより，とても暑かった年やとても寒かった年がどの辺にあったかよくわかります。たとえば，1993年ごろにものすごく寒かった年があるのがわかりますが，これは1993年が記録的な冷夏だったときです。

あと，今の時点では蛇足に感じるかもしれませんが，偏差の平均はゼロ という超重要な性質があるので，覚えておいてください。

In [7]:
np.mean(tempa) # 数値計算の誤差はあるがほぼゼロ
Out[7]:
2.3931474086364484e-16

SST分布の偏差を描画してみる

さて，先ほどまで話してきた日本の異常気象をもたらす気候のゆらぎの一つとして，エルニーニョ現象があります。エルニーニョ現象は，東太平洋のペルー沖の海面水温が例年より暖かくなる現象です。

この現象の様子を調べるために，次は海面水温偏差を描画してみましょう。

前回までの海面水温データを持ってくるか，こちらからダウンロードしてください（出典は第2回に記述）：http://web.is.ocha.ac.jp/~tsubasa/env_info/sst_OISST.npz

先週計算しておいた海面水温気候値はこちらにもあります：http://web.is.ocha.ac.jp/~tsubasa/env_info/sstc_OISST.npz

データセットの読み込み（前回までの復習）

先週と同じように海面水温のデータを読み込みます。このとき，先週計算しておいた海面水温気候値も読み込んでおきます。

In [8]:
loadfile1 = 'sst_OISST.npz' # 海面水温データの入力ファイル名を定義
sst_dataset = np.load(loadfile1) # データセットはまずデータセットごと入力します 
sst = sst_dataset['sst'] # 海面水温(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に保存

loadfile2 = 'sstc_OISST.npz' # 先週計算した気候値の入力ファイル名を定義
sstc_dataset = np.load(loadfile2) # データセットはまずデータセットごと入力します 
sst_clim = sstc_dataset['sst_clim'] # 海面水温気候値を変数sst_climに保存

SST偏差の計算と保存

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

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

# 0で埋められた行列を使って，欲しいサイズの行列を作っておく（初期化）
ssta = np.zeros((imt, jmt, tmt)) 

# 偏差の計算
## やはり1982年1月から順番にfor文を回していく
## 1次元目（経度），2次元目（緯度）の方向には手をつけない
## →「配列の要素全部」という意味でコロンを書いておけばOK
## 3次元目（時間の方向）について，sstからsst気候値を引き算
for yy in range(1982, 2020):
    for mm in range(1, 13):
        ssta[:, :, (y==yy)*(m==mm)] = \
           sst[:, :, (y==yy)*(m==mm)] - (sst_clim[:, :, mm-1])[:, :, np.newaxis]
# np.newaxisは，360x180の配列を360x180x1として扱うためのコマンドです。
# （3次元の配列から2次元の配列を引き算することはできないため）

# 偏差は次回使うので保存しておく
savefile = 'ssta_OISST.npz'
np.savez(savefile, ssta=ssta, lon2=lon2, lat2=lat2, y=y, m=m)

ある月のSST偏差の描画

いつも通り描画します。だいぶ慣れてきた頃でしょうか。

In [10]:
# 描画したい年・月
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(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 = str(draw_year) + '/' + str(draw_month)
plt.title(title)
plt.xlim(0, 360)
plt.ylim(-90, 90)
plt.show()

第2回でちょっとだけ勉強しましたが，1997年12月は，大きなエルニーニョ現象があった年です。

東太平洋のペルー沖の海面水温が例年より5度ほど暖かくなっています。

この規模のエルニーニョ現象が発生したのは，1982年，1997年，2015年の3回のみです。「ゴジラエルニーニョ」などと言われました。

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

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

A問題. 1960年から1989年の東京の気温のデータで上記と同様の解析を行うことで，冷夏・猛暑・暖冬・寒冬などをいくつか見つけてください。wikipediaなどを引用し，自分の解析が正しいことを確かめてください。

B問題. 2010年はラニーニャ現象が起きた年です。2010年の12月の海面水温偏差を描画することで，ラニーニャ現象とはどのような現象なのかを自分の言葉で簡潔に説明してください（メカニズムの話とかを調べてもらう必要は要りません。描画した海面水温の図を見て，思ったことを書いてください）。

C問題. 偏差の平均がゼロになることを，数学的に示してください。このとき，Markdownで数式を書くためにはLaTeXのコマンドが使えますので，これを機に勉強してみてください（どうせ卒研でやらされます）。

D問題. 日本政府観光局のサイト https://www.jnto.go.jp/jpn/statistics/visitor_trends/ に，2003年-2020年の月別訪日外国人数のデータが公開されています。このデータを気候値や偏差と同様の考え方で解析することによって，一年のうちどの時期に訪日外国人数が多い傾向にあるか，また過去18年間の間に訪日外国人が特に多かったり少なかった時期があれば，理由を考えてみてください。

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