環境情報論　第10回 地図の描画と気象のテレコネクション

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

今回は，まず必要な地図ツールの使い方も学んだあと，回帰図や相関図の考え方を用いて遠くの地域の気象がどのように繋がっているかを理解します。

cartopyのインストール

Jupyter Notebookにてインストールするのは，通る場合はあっさり通るらしいのですが，インストールがうまくいかない場合も結構あるみたいです。

その場合は，今回だけGoogle Colaboratoryを用いて課題を提出してください（使用は簡単なので心配しないでください）。

Jupyter Notebookの場合

第0回で環境構築したときに無事にAnacondaが入っていれば，ターミナルに次のように打ち込むだけでインストールが完了します。

conda install -c scitools cartopy

（インストール中に[y]?n?と聞かれたら，yesの意味でyと打ち込んでください。）

Pythonのバージョンで怒られた方は，

conda install -c conda-forge cartopy

でうまくいく場合もあるようです（うまくいかないときは何回か試してJupyter Notebookを再起動するなどもしてみてください）。

もし独自の環境構築を行った方で，Anacondaを使っていなければ（上記のcondaコマンドが通らなければ），以下のサイト等を参考にしてインストールしてみてください。

https://yoheikoga.github.io/2017/09/19/jupyter-notebook-cartopy-map-test/

Google Colaboratoryの場合

Jupyter Notebookでインストールがうまくいかない場合も，Google Colaboratoryを使ってください。

神山研の学生さんが書いてくださったマニュアルがあるので，こちらを参考にしてください。

http://web.is.ocha.ac.jp/~tsubasa/env_info/colab_cartopy.html

Anacondaを入れていない人（pipでインストールしたい人）、かつWindowsユーザの場合

こちらの方法がうまくいくということを学生さんに教えていただきました。

http://py3dtech.com/?p=61

地図の描画

モジュールのインポート

In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import Normalize
# 以下の3つが新しい
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LatitudeFormatter,LongitudeFormatter
import matplotlib.ticker as mticker

ModuleNotFoundError: No module named 'cartopy'

というエラーが出てしまった人は，ターミナルからjupyter notebookと書いてJupyter Notebookを起動するのではなく，Anaconda NavigatorからGUIでJupyter Notebookを起動すると良いらしいです。

データセットの読み込み

データはここです（出典は第2回参照） http://web.is.ocha.ac.jp/~tsubasa/env_info/sst_OISST.npz

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

正距円筒図法

それでは，正距円筒図法で1997年12月の海面水温分布を地図上に描画してみましょう。このとき，（いままで書かずにお茶を濁していた）海岸線も書きます。

＜参考＞Cartopyについて困ったことはかなりこちらに書いてあって，大変勉強になりました。https://yyousuke.github.io/matplotlib/index.html

In [3]:
# 描画したい年・月・変数
draw_year = 1997
draw_month = 12
draw_var = np.squeeze(sst[:, :, (y==draw_year)*(m==draw_month)])

# 図の中心経度
c_lon = 180

# vminはカラーバーの下限，vmaxはカラーバーの上限
# vintはカラーバーの間隔
vmin = -5
vmax = 35
vint = 5
# 深い青から深い赤に向かうカラーバーを指定
cm = plt.get_cmap('seismic') 

# 描画する枠を作る
fig = plt.figure()

# 枠の中に絵を入れる（図を1枚しか書かないときは1, 1, 1でOK）
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree(central_longitude=c_lon))

# 色で塗られた等高線を描く（前回とちょっとだけ違う）
cs = ax.contourf(lon2, lat2, draw_var,\
                  cmap=cm, norm=Normalize(vmin=vmin, vmax=vmax),\
                  levels=np.arange(vmin,vmax+vint,vint), extend='both', \
                  transform=ccrs.PlateCarree())
                  #最後の一行が大事！！！（lon2, lat2は正距円筒の書き方なので，変換が必要）

# 海岸線を書く（これが気象場を描画するときにめちゃくちゃ重要になる）
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, 359.9, -90, 90], crs=ccrs.PlateCarree())

だいぶプロフェッショナルな図に見えるようになったと思います。

ランベルト正角円錐図法

ほぼお遊びですが，こんな図もかけます。ランベルト正角円錐図法 (LambertConformal)と呼ばれる図法です。

（ちょっと描画に時間がかかります）

In [4]:
fig = plt.figure()
ax = plt.axes(projection=ccrs.LambertConformal())

# 色で塗られた等高線を描く
cs = ax.contourf(lon2, lat2, \
                  np.squeeze(sst[:, :, (y==draw_year)*(m==draw_month)]), \
                  cmap=cm, norm=Normalize(vmin=vmin, vmax=vmax),\
                  levels=np.arange(vmin,vmax+vint,vint), extend='both', \
                  transform=ccrs.PlateCarree())  
                  #最後の一行が大事！！！（lon2, lat2は正距円筒の書き方なので，変換が必要）

# 海岸線を書く
ax.coastlines(lw=0.5,color='gray', resolution='50m')

# 緯度経度線を書く
ax.gridlines(xlocs=mticker.MultipleLocator(10), \
               ylocs=mticker.MultipleLocator(10), \
               linestyle='-', color='gray')

plt.colorbar(cs, shrink = 0.8) # カラーバーをつける
Out[4]:
<matplotlib.colorbar.Colorbar at 0xb23ea89b0>

他にも，http://metpost.hatenablog.com/entry/2015/11/05/180006 あたりを参考に，色々遊んでみてください（A問題）。

気象のテレコネクション

大気を介して，遠くの気象が影響を及ぼすことをテレコネクション（遠隔影響）といいます。

今日は，東京付近の気圧変動はどのような場所の気象と関係があるかを調べてみましょう。

まず，海面更正気圧（=地面があるところも海抜0mだと思って補正した気圧）の偏差データを読み込みます。

データはここです http://web.is.ocha.ac.jp/~tsubasa/env_info/slpa_erai.npz

（データの出典：https://apps.ecmwf.int/datasets/data/interim-full-moda/levtype=sfc/ ）

In [5]:
loadfile = 'slpa_erai.npz' # 入力ファイル名を定義
slpa_dataset = np.load(loadfile) # データセットはまずデータセットごと入力します 
slpa = slpa_dataset['slpa'] # 海面更正気圧偏差(sea level pressure anomalies)を変数slpaに保存
lon2 = slpa_dataset['lon2'] # 経度(longitude)を変数lon2に保存（2は「2次元配列」の意味）
lat2 = slpa_dataset['lat2'] # 緯度(latitude)を変数lat2に保存
y = slpa_dataset['y'] # 年(year)を変数yに保存
m = slpa_dataset['m'] # 月(month)を変数mに保存

次に，東京付近の1点の気圧偏差をインデックスとして，気圧の回帰図を計算してみます（一点回帰と呼ばれる手法です）。

まず，lon2とlat2を用いて，東京付近の格子が配列の何番目なのかを特定します。

In [6]:
# なんだかんだlon2とlat2の中身を理解して数えるより
# 適当に数字を入れて探すのが早かったりする（ブルートフォース作戦）
print(lon2[46, 1]) 
print(lat2[1, 18]) 

138.0
36.0

経度方向は46番目，緯度方向は18番目に，東京付近（東経138度，北緯36度）の点があることがわかったので，その点の気圧偏差をslpa_Tokyoというインデックスとして，気圧偏差の回帰図を作成します。

In [7]:
# 東京付近の点の気圧偏差を抜き出してインデックスとする
slpa_Tokyo = slpa[46, 18, :] #東経138度，北緯36度

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

# 0で埋められた行列を使って，欲しいサイズの行列を作っておく（初期化）
a_slpa = np.zeros((imt, jmt)) # 回帰係数a
b_slpa = np.zeros((imt, jmt)) # 切片b
r_slpa = np.zeros((imt, jmt)) # 相関係数r

# 回帰図と相関図の計算
# 経度方向にimt（=360）回，緯度方向にjmt（=180）回forループを回す
for ii in range(0, imt):
    for jj in range(0, jmt):
        [a_slpa[ii, jj], b_slpa[ii, jj]] = np.polyfit(slpa_Tokyo, slpa[ii, jj, :], 1)
        r_slpa[ii, jj] = np.corrcoef(slpa_Tokyo, slpa[ii, jj, :])[1, 0] 

    if (ii % 30 == 0):
        print(ii) # ちゃんと計算が進んでるかチェックするために，30回に1回iiを出力する

0
30
60
90

一点回帰図を描画してみます。

In [8]:
draw_var = a_slpa

# 一番右が東経357度までしか入っていないので，白い線が入ってしまう。
# よって，360度の行を追加する必要がある。
# lon2には360を入れ，lat2とdraw_varには0度の行と同じものを入れる
lon2 = np.vstack([lon2, (lon2[0, :]+360.0)[np.newaxis, :]])
lat2 = np.vstack([lat2, (lat2[0, :])[np.newaxis, :]])
draw_var = np.vstack([draw_var, (draw_var[0, :])[np.newaxis, :]])

fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection=ccrs.PlateCarree(central_longitude=180))

# 黒い等高線を描く
cs = ax.contour(lon2, lat2, draw_var, 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)

# 描画範囲の指定
# 正距円筒図法であることを明示するために2つ目の引数が必要
# 経度の指定では，c_lon±180の値までしか書いてはいけないっぽい
# 0, 360と書くと360が0と同じ数字だと思われて地図が棒になる
ax.set_extent([0, 359.9, -90, 90], crs=ccrs.PlateCarree())

実線で描かれているのが正の回帰係数，破線で描かれているのが負の回帰係数です。

北極やカナダ，ヨーロッパの方にまで負の回帰係数が見えます。分かりやすいように，北極からみた図を書いてみます。

In [9]:
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection=ccrs.NorthPolarStereo())

# 色で塗られた等高線を描く（前回とちょっとだけ違う）
cs = ax.contour(lon2, lat2, draw_var, colors=['black'],transform=ccrs.PlateCarree())

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

ax.set_extent([-180, 180, 20, 90], ccrs.PlateCarree())

東京付近の月平均気圧は，北極を中心とした気圧変動と強い逆位相の関係を持つことがわかります（北極振動とか北半球環状モードと呼ばれています）。

この図を見ると，大気のゆらぎを介して日本の天気は遠く北アメリカやヨーロッパと繋がっていることがわかります。テレコネクションの典型例です。

テレコネクションを調べる際には，回帰図・相関図や，コンポジット図が非常によく用いられます。

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

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

A問題. cartopyの機能をググって，上記の例以外の好きな投影法で，1つ海面水温の図を描いてみてください。

B問題. 東京付近以外の点を一つ選び，気圧偏差の一点回帰を行うことで，面白いテレコネクションを見つけてください。

C問題. 第2回のC問題では，日本地図を描かずに日本付近の天気を描画しました。cartopyで日本地図を重ねて，第2回のC問題をやり直してみてください（天気予報で見るような図に少し近づいたでしょうか...！）。

D問題. 渡部雅浩『テレコネクションー日本の天候を左右するものー』 http://www.metsoc-hokkaido.jp/saihyo/pdf/saihyo48/saihyo48-080.pdf を読んで，

  1. Wallace and Gutzler (1981) は，テレコネクション・パターンを検出するためにどのようなデータ解析手法を用いたか

  2. 日本の天候はテレコネクションにどのような影響を受けるか

の2点について，自分の言葉で簡単に説明してください。

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