はじめに¶
前回は、KEOの時系列を例にNumPyの1次元配列の基本的な取り扱いを学びました。 一方、シミュレーション結果のようなデータの場合は空間2, 3次元に時間1次元が加わって、最大4次元のデータを扱うことになります。 NumPyの配列は多次元データの扱いにも長けており、データの切り出しや統計に便利な機能が備わっています。
15年にわたる北西太平洋域の海面水温の時系列データを例に、多次元データの取り扱いを学びましょう。 まずは、NumPyとmatplotlibをインポートします。
import numpy as np
import matplotlib.pyplot as plt
a = np.arange(20).reshape((4,5)) # 4 x 5の2次元配列
print(a)
b = np.arange(27).reshape((3,3,3)) # 3 x 3 x 3の3次元配列
print(b)
形状を表す属性¶
配列xの形状を知るとき、以下の3つの属性をよく使います。
x.ndim次元数x.shape形状（例えばN×M×Lのとき、(N, M, L)というタプル）x.size合計要素数
print(a.ndim)
print(a.shape)
print(a.size)
print(b.ndim)
print(b.shape)
print(b.size)
# 例：2次元の場合
print(a) # 上で定義した4x5配列。後ろの次元ほど内側のカッコで表示される
print(a[0,0]) # 第1次元が（0から数えて）0番目、第2次元が0番目の要素
print(a[0,1]) # 第1次元が（0から数えて）0番目、第2次元が1番目の要素
print(a[1,0]) # 第1次元が（0から数えて）1番目、第2次元が0番目の要素
print(a[1,4]) # 第1次元が（0から数えて）1番目、第2次元が4番目の要素
print(a[1,-1]) # 第1次元が（0から数えて）1番目、第2次元が最後から1番目の要素
複数の要素を取り出す¶
複数要素を取り出すスライスを、次元ごとに指定することもできます。
print(a[1,:]) # 第1次元が1番目の、すべての要素（1次元配列として得られる）
print(a[1,1:4]) # 第1次元が1番目、第2次元が1番目から4番目の手前
print(a[:,3]) # 第2次元が3番目のすべての要素（1次元配列として得られる）
練習4-1¶
はじめに、次の通りに2次元配列xを作りなさい。
x = np.arange(48).reshape((8,6))
print(x)
xから次の要素を得るには、どのように指定すればよいでしょうか。
32040[24, 25, 26, 27, 28, 29][0, 6, 12, 18, 24, 30, 36, 42][[18, 19, 20, 21, 22, 23], [24, 25, 26, 27, 28, 29]][[38, 39, 40, 41], [44, 45, 46, 47]]
p = np.arange(30)
print(p) # 0~29の数列を30要素の1次元配列に格納したもの
print(p.shape)
q = p.reshape((6,5)) # (6,5)というタプル（リストみたいなもの）で形状を指定している
print(q) # pを6×5の2次元配列に格納したもの
print(q.shape)
r = p.reshape((2,3,5)) # (2,3,5)というタプルで形状を指定
print(r) # pを2×3×5の3次元配列に格納したもの
print(r.shape)
ravel：形状を1次元にする¶
上の例とは逆に、多次元配列を1列に並べなおして1次元配列にするときにはravelメソッドを用います。1次元に並べなおすときにはa[0,0], a[0,1], ..., a[0,n], a[1,0], a[1,1], ..., a[1,n], ...のように、後ろの次元を小さい単位（一の位）とし、前の次元を大きい単位（十の位、百の位・・・）のように並べられます。
print(r)
print(r.ravel()) # rを1次元に並べたもの。r[0,0,0~4], r[0,1,0~4], r[0,2,0~4], r[1,0,0~4]...
転置¶
配列の転置はa.Tとすることで得られます。例えば2次元配列でb=a.Tのとき、b[i,j]はa[j,i]に等しくなり、bの形状はaの形状の1,2次元目の順序を入れ替えたものとなります。
()が無いことからわかる通り、転置配列は厳密にはメソッドではなく属性です。
print(q)
print(q.shape)
s = q.T
print(s)
print(s.shape)
最後に注意として、reshape, ravel, Tのメソッドや属性は配列そのものを上書き更新するわけではなく、操作を施したものを返していることに気をつけましょう。すぐ上の例でも、q.Tをsに代入することでsはqの転置となっていますが、q自体はもともとのままです。
print(a)
print(a.sum()) # 全体の合計をとり、定数で返す
print(a.sum(0)) # 1次元めの方向に和をとり、1次元配列で返す
print(a.sum(1)) # 2次元めの方向に和をとる、1次元配列で返す
print(b) # 3x3x3のルービックキューブのような配列
print(b.mean()) # 全体の平均をとり、定数で返す
print(b.mean(0)) # 1次元めの方向に平均をとり、2次元配列で返す
print(b.mean((1,2))) # 2次元めと3次元めの方向に平均をとり、1次元配列で返す。複数の次元を指定するときは、タプルで指定する
練習4-2¶
まず、次の通りに配列yを定義しなさい。（参考：第3回練習問題）
これは、あるマンションの部屋ごとの居住者数を表しているとします。マンションは7階建てで、各階に8部屋あります。
yの第1次元は階（1階は0, 2階は1, ..., 7階が6）を、第2次元は部屋番号（1号室が0, 2号室が1, ..., 8号室が7）を表しています。
rng = np.random.default_rng(seed=54321)
y = rng.integers(low=1, high=7, size=56).reshape((7,8))
以下を求めなさい。
- 2階3号室に住んでいる人数（5）
- マンション全体の居住者数の合計（201）
- マンション全体の一部屋の居住者数の平均（3.589...）
- 各階の居住者数の合計（31 31 25 25 36 26 27）
- 各階の居住者数の最大（6 6 5 6 6 5 6）
- 1号室に住む人数の合計（22）
- 4階以上に住む1~4号室の人数の合計（59）
- 4階以上に住む5~8号室の人数の合計（55）
練習4-3¶
n個のサイコロを振って出た目の平均を$\bar{X}$とすると、$\bar{X}$の期待値は3.5ですが、その標準偏差は$\sqrt{n}$に比例して小さくなります。そのことを確かめましょう。
n個のサイコロを振って$\bar{X}$を求める操作を100回行い、その100個の$\bar{X}$の平均値・標準偏差を出力します。
まず、練習4-2と同様に1から6の数列を返す乱数生成器rngを次のように定義します。（今回はseedなしで行きましょう）
rng = np.random.default_rng()
整数nを引数とし次のような操作を行う関数mean_and_stdを定義しなさい。
rng.integers(low=1, high=7, size=n*100)で、1~6の乱数が100n個連なった配列を得る- 上の配列をn×100の2次元配列に変形する
- 大きさnの次元に平均をとり、大きさ100の1次元配列を得る
- 上の配列の平均と標準偏差を
printする
mean_and_stdを$n=10, 100, 1000, 10000, 100000$について呼び、確かに標準偏差が$\sqrt{n}$に比例して小さくなっていることを確かめましょう。
実例：海面水温の時空間データ¶
データの読み込み¶
それでは、実際の多次元データを使ってデータ解析をしてみましょう。 今回は、ヨーロッパ中期予報センター（ECMWF）が制作した再解析プロダクトERA5の海面水温データを用います。
今回使うデータ¶
再解析データとは、観測と数値シミュレーションのデータの持つ以下の長所・短所を補い合うため、これらのデータを組み合わせて得られた推測値です。
- 観測：現実を反映しているが、時間的・空間的に飛び飛びにしかない
- 数値シミュレーション：時間的・空間的に抜けがないが、現実からかけ離れるおそれがある
数値シミュレーションと同様、空間2or3次元+時間1次元の多次元データとして得られますので、今回紹介した多次元配列の処理が必要となります。
今回用いるデータは以下のような物です。
- 東経120～180度、北緯20～60度の北西太平洋域、0.25度ごと
- 2007年1月から2021年12月まで、毎月の平均海面水温データ
こちらからダウンロードし、保存してください。
https://yasushifujiwara.github.io/dataanalysistutorial/data_ERA5/sst_nwpac.nc
NetCDF形式とは¶
今回のデータは拡張子.ncをもつ、NetCDF形式で保存されたバイナリデータです。NetCDFの特徴はすぐ下で述べますが、大気や海洋分野で非常によく使われるデータ形式です。
以下の方法で読み込んでみましょう。
import netCDF4
data = netCDF4.Dataset(r"C:\Users\yfuji\Desktop\lesson\data\data_ERA5\sst_nwpac.nc")
NetCDF形式の特徴として、自己記述的であることが挙げられます。
例えば（今回の場合）時間1次元・空間2次元の格子状に格納されたSSTデータだけでは、それぞれの点がいつ・どの経度・どの緯度の情報なのかわかりません。
NetCDFデータは目的のデータに加え、緯度・経度・時間といった軸データも同時に含んでいます。
読み込んだdata変数を出力することで表示される、以下の情報のなかのdimensionsやvariablesに注目してみましょう。
print(data)
dimensionsにはこのファイルの中で登場する次元と、それぞれの次元が持つ点の数が表記されています。今回の場合、経度longitudeが241点、緯度latitudeが161点、時間timeが180点あることがわかります。
variablesには、このファイルに含まれているデータと、そのデータがどの次元を持つのかが表記されています。
軸データであるlongitude, latitude, timeは自分自身の次元を持つ1次元配列で、海面水温sstは時間・緯度・経度の順に次元を持っていることがわかります。
データの中身へのアクセス¶
さて、それではそれぞれのデータにアクセスしてみましょう。以下のようにすることで、それぞれのデータをndarrayとほぼ同じ形（厳密にはndarrayを拡張したmaskedarray）で取得できます。
lon = data["longitude"][:]
lat = data["latitude"][:]
time = data["time"][:]
sst = data["sst"][:,:,:]
# masked arrayを普通のndarrayに戻す。陸地の値をnanにする
sst = sst.filled(np.nan)
lon # 経度
sst.shape # SST配列の形状
数値の意味を調べる¶
ある地点の海面水温を覗いてみましょう。適当に、時間0番目, 緯度80番目、経度120番目の点を見てみます（チェックする箇所はどこでも問題ありません。）
その時間・緯度・経度がどこを指しているのかも、time, lat, lonの配列の同番号から調べてみます。
print(sst[0,80,120]) # SST
print(time[0]) # 時間
print(lat[80]) # 緯度
print(lon[120]) # 経度
緯度・経度は北緯40度・東経150度だと想像できますが、SSTと時間が何の数字なのかよくわかりませんね。
自己記述的であるNetCDFファイルには、それぞれの数値がどんな単位で格納されているのかの説明も含まれています。
元データdata変数に立ち戻り、以下のように調べてみましょう。
print(data["sst"])
この中で、unitsに注目しましょう。今回の場合、単位（unit）はKであることがわかります。
上空で気温が非常に低くなる大気分野では℃にこだわらずKを使うのが便利なことが多いのですが、最低水温が結氷温度（-2℃くらい）で制限される海洋分野では℃標記のほうが用いられることが多いです。ここではSSTを℃表記で扱うことにしましょう。
sst = sst - 273.15 # KからCelciusへ
print(sst[0,80,120])
さて、次は時間937944が何を表すのかを調べましょう。
print(data["time"])
unitsを見てみると、「1900年1月1日0時0分0秒から数えた時間（hour）」と書いてあります。つまり、時間937944は「1900年1月1日0時0分0秒から937944時間後」を表します。
datetimeモジュールを使って日付を扱う¶
これが具体的にどの日を表すのか計算するのは、月ごとの日数が違ったりうるう年を考慮したりでなかなか面倒なのですが、現代は便利なものでネットですぐ調べられます（答え：2007年1月1日0時）。
また、以前存在だけ紹介したdatetimeモジュールを使ってもすぐ調べられます。参考までに方法を紹介します。
from datetime import datetime, timedelta
# datetimeはある日時、timedeltaは日時の差を表す
# datetime + timedelta -> datetime
# timedelta + timedelta -> timedelta
# datetime + datetime -> エラー
print(937944 * timedelta(hours=1) + datetime(1900, 1, 1))
等値線図（空間2次元のプロット）¶
基本の描き方¶
それでは、いよいよ海面水温データの中身を図示しましょう。
sstは時間1次元・空間2次元の3次元配列ですが、例えばsst[0,:,:]などのようにある時間を切り出して、空間2次元の平面図を描くのが一番わかりやすいと思います。
2次元的なデータ分布を平面図に描くには、同じ値の点を線で結んだ（あるいは同じ色で塗った）等値線図を用いるのが有力な方法です。matplotlibの機能で、等値線図を描いてみましょう。
等値線図を描く基本の関数はplt.contour（線を引く）またはplt.contourf（色を塗る）です。基本的な使い方はplt.contourf(X, Y, Z)で、Zが描こうとしているデータの入った2次元配列、XとYがZと同じ形状の配列で、Zの格子のx座標, y座標を表します。
長方形状に格子が配置されている場合、2次元配列X, Yはnp.meshgrid関数を用いることで1次元配列から生成できます。
# 2次元の位置座標配列を生成。
Lon, Lat = np.meshgrid(lon, lat) # SSTの配列の次元がlat, lonの順の場合。
# Lon, Lat = np.meshgrid(lon, lat, indexing="ij") # SSTの配列の次元がlon, latの順の場合。
print(lat.shape)
print(lat)
print(Lat.shape)
print(Lat)
# 時間0について図を描く
plt.contourf(Lon, Lat, sst[0,:,:], # 時間0について、色塗り等値線図を描く
levels=31, # 色分けは30段階（色の区切りが31個）で行う
cmap="plasma" # 色調に"plasma"を用いる
)
plt.colorbar() # カラーバーを図の右側に表示する
南から北に水温が下がっていたり、黒潮が暖水を運んでいる様子が見て取れますね。
陸地が白く塗りつぶされていますが、これは陸地の点のsst配列にnanが入っているからです。より地図っぽくするには海岸線を線で引いてやるとよいのですが、ここではそこまで踏み入らないことにします。
また、levels引数に色分けの個数（あるいは線の本数）を表す整数を与えています。こうすると色塗りの温度範囲を自動で判別してくれますが、自分で色塗り範囲を指定したいときにはlevelsに色を区切る温度を表す1次元配列を与えることもできます（違うデータを比較するときは色範囲を統一するべきなのでよく使います）。
cmap引数には、色調列の名前を与えられます。matplotlibで使える色調列はこちらなどから見られます。
余談ですが、よく使われる虹色の"jet"は前線域の誤認につながったり色覚異常の方に見えづらいことから使うのを避けるべきとされています。（参考：https://tos.org/oceanography/assets/docs/29-3_thyng.pdf ）
誰にでも見やすく、データの意味をよく伝える色調列を選ぶようにしましょう。
練習4-4¶
上を参考に、適当な時間について海面水温の分布図を描きなさい。
- 時間0は2007年1月、時間1は2007年2月...となっています。自分が描いているのは何年何月の図ですか？
levelsに与える数を変えると、見え方はどのように変わるでしょうか。どの数を選ぶのが見やすいでしょうか。
関数を使って楽をする¶
違う時刻の図を見たいときにいちいち上のコードを書き直すのは面倒なので、図を描く関数を作ってしまうと楽です。
datetimeモジュールも駆使して、日時入りの図をつくる関数を定義してみます。
def draw_sst(n):
lev = np.linspace(-2, 33, 36) # -2度から33度まで、1度区切りの配列
plt.figure() # 新しいキャンバスを立ち上げる; 複数枚図を描くとき必要
plt.contourf(Lon, Lat, sst[n,:,:], # 時間nについて、色塗り等値線図を描く
levels=lev, # 色分けはlevに従って行う
cmap="plasma" # 色調に"plasma"を用いる
)
plt.colorbar() # カラーバーを図の右側に表示する
plt.title(time[n] * timedelta(hours=1) + datetime(1900, 1, 1)) # 日時を図タイトルに表示する
plt.show() # 図を表示する; 複数枚図を描くとき必要
draw_sst(169)
draw_sst(175)
sst_mean = sst.mean(0)
lev = np.linspace(-2, 33, 36)
plt.contourf(Lon, Lat, sst_mean, # sst_meanについて、色塗り等値線図を描く
levels=lev, # 色分けは上で定義したlevに従い行う
cmap="plasma" # 色調に"jet"を用いる
)
plt.colorbar() # カラーバーを図の右側に表示する
渦が取り除かれ、平均的でなだらかな水温分布が見えてきました。
基本的には南から北へと水温が下がっていきますが、黒潮や対馬海流に乗って周囲より暖かい水が北上している様子が見えます。つまり、これらの海流は常在している特徴と考えられるでしょう。
平均からのズレ¶
さて、平均的な描像がわかったところで今度は偏差 (anomaly)、すなわち平均からのズレを図にしてみましょう。
例えば、2017年以降のデータなら黒潮大蛇行の特徴が捉えられるかもしれません。
def draw_sst_a1(n):
lev = np.linspace(-17,17,35)
plt.figure() # 新しいキャンバスを立ち上げる; 複数枚図を描くとき必要
plt.contourf(Lon, Lat, sst[n,:,:]-sst_mean, # 時間nについて、色塗り等値線図を描く
levels=lev, # 色分けは1度ごとに行う
cmap="seismic" # 色調に"seismic"を用いる
)
plt.colorbar() # カラーバーを図の右側に表示する
plt.title(time[n] * timedelta(hours=1) + datetime(1900, 1, 1)) # 日時を図タイトルに表示する
plt.show() # 図を表示する; 複数枚図を描くとき必要
draw_sst_a1(133)
draw_sst_a1(139)
・・・さて、これは「通常からのズレ」をよく表しているのでしょうか？そもそも通常とは？
気候値¶
平均からの偏差の図を見ると、夏は単に赤っぽく、冬は単に青っぽくなってしまいました。これは年間平均からの偏差を調べたからで、当然夏は平均より高く、冬は平均より低く出てしまいます。
ですが、多くの場合見たいのはそういった季節の違いによる差ではなく、夏に水温が高いのも冬に水温が低いのも「通常」だと考えるのが自然です。
このような背景から、大気・海洋の分野では基準的な状態として「同じ月の長期平均（例えば過去30年間の1月の平均）」がよく用いられます。これを気候値と呼びます。
ここでは15年分のデータを使っていますので、「15年分の1月の平均」「15年分の2月の平均」...「15年分の12月の平均」を気候値として求めましょう。
その結果の配列は(時間12点, 緯度161点, 経度241点)の配列となるはずです。
さて、この計算はどのようになされるでしょうか。
練習4-5¶
sst (時間180点, 緯度161点, 経度241点)を月ごとに平均し、(時間12点, 緯度161点, 経度241点)の形状を持つ気候値の配列sst_climを作りなさい。
以下で3つの解答例を示しますが、これまでの知識で十分にできることですのでまずは自分で考えてみましょう。
気候値描画の関数¶
解答例の前に、気候値ssh_climを図示する関数を定義してしまいます。
def draw_sst_clim(month):
# monthは何月か
# sst_climは0: 1月、1: 2月、...、11: 12月なので1個ずつずれてる
month_name = ["Jan", "Feb", "Mar", "Apr", "May", "Jun",
"Jul", "Aug", "Sep", "Oct", "Nov", "Dec"]
lev = np.linspace(-2, 33, 36) # -2度から33度まで、1度区切りの配列
plt.figure() # 新しいキャンバスを立ち上げる; 複数枚図を描くとき必要
plt.contourf(Lon, Lat, sst_clim[month-1,:,:], # 時間nについて、色塗り等値線図を描く
levels=lev, # 色分けはlevに従って行う
cmap="plasma" # 色調に"plasma"を用いる
)
plt.colorbar() # カラーバーを図の右側に表示する
plt.title(month_name[month-1] + " climatology") # 月名を図タイトルに表示する。"climatology"は気候値のこと
plt.show() # 図を表示する; 複数枚図を描くとき必要
練習4-5の解答例¶
# 解答例1: forループを用いて、時間軸12個とばしの平均をとる
sst_clim = np.zeros((12, 161, 241)) # 入れ物をつくる
for n in range(12):
sst_clim[n,:,:] = sst[n::12,:,:].mean(0)
# n=0なら 0, 12, 24, ...
# n=1なら 1, 13, 25, ...
# 解答例2: forループとマスクを用いて、月の通し番号を12で割った余りごとに平均をとる
sst_clim = np.zeros((12, 161, 241)) # 入れ物をつくる
nmonth = np.arange(180) # 月の通し番号, 0からスタート
for n in range(12):
mask = (nmonth % 12 == n) # 月番号を12で割った余りがn
sst_clim[n,:,:] = sst[mask,:,:].mean(0)
# 解答例3: reshapeで時間軸を総月数180から年数15x月数12に変形して（折りたたんで）、年数方向に平均をとる
sst_clim = sst.reshape((15, 12, 161, 241)).mean(0)
draw_sst_clim(2) # 2月
draw_sst_clim(8) # 8月
気候値からの偏差を描画する¶
「通常の状態」として気候値を求めたところで、改めて「偏差」を描く関数を定義しましょう。
通し番号nのデータについて、同月の気候値からの偏差を求めます。
def draw_sst_a2(n):
lev = np.linspace(-5,5,21)
plt.figure() # 新しいキャンバスを立ち上げる; 複数枚図を描くとき必要
plt.contourf(Lon, Lat, sst[n,:,:]-sst_clim[n%12,:,:], # 時間nについて、色塗り等値線図を描く
levels=lev, # 色分けはlevに従って行う
cmap="seismic" # 色調に"seismic"を用いる
)
plt.colorbar() # カラーバーを図の右側に表示する
plt.title(time[n] * timedelta(hours=1) + datetime(1900, 1, 1)) # 日時を図タイトルに表示する
plt.show() # 図を表示する; 複数枚図を描くとき必要
では、先ほどと同様に2018年2月と8月の平年からの偏差を見てみましょう。
draw_sst_a2(133)
draw_sst_a2(139)
先ほどは季節変化で見えにくくなってしまっていた、様々な構造が見えています。
2月には東海・近畿南岸に大きな冷水塊が見えています。これは黒潮が蛇行し、その北側に普段はない冷水が位置しているから見えているものです。
また、2月と8月を比べると2月は数百km（数度）規模の構造が目立つのに対し、8月は数千km（数十度）規模の正負の構造が目立っています。
前者は海洋の中規模渦に起因する水温変動なのに対し、後者は大気の気圧配置に関連しています。夏は強い季節躍層が発達し、海洋内部の状態よりも大気の状態に左右されやすくなるのではと推測できます。
等値線図（時間1次元・空間1次元のプロット）¶
空間分布の時間変化を見るために¶
ここまでは、プロットする時間を1つ選んで、空間2次元の等値線図を見てきました。
空間分布の時間変化を示したい場合には時間1次元・空間1次元の等値線図を描くこともできます。このような図をホフメラー図と呼びます。
ここではある緯度を選び、横軸に経度、縦軸に時間を取ったホフメラー図を描いてみましょう。
まずは、x座標（経度）・y座標（時間）の値が入った2次元配列を作ります。簡易的にデータを見る場合には937944などが入った配列timeを使っても良いですが、ここではdatetimeを使って具体的な年を見ることにしましょう。
# Lon, Time = np.meshgrid(lon, time)
Lon, Time = np.meshgrid(lon, time * timedelta(hours=1) + datetime(1900, 1, 1))
ホフメラー図を描く関数¶
緯度配列の番号を指定し、その緯度でのSSTの時間変化をホフメラー図で描いてくれる関数を定義します。
def drawhov_sst(n):
plt.figure(figsize=(6,8)) # 新しいキャンバスを立ち上げる; 見やすくするため、ちょっと縦長にする
plt.contourf(Lon, Time, sst[:,n,:], # 緯度nについて、色塗り等値線図を描く
levels=lev, # 色分けはlevに従って行う
cmap="plasma" # 色調に"plasma"を用いる
)
plt.colorbar() # カラーバーを図の右側に表示する
plt.title("latitude="+str(lat[n])) # 緯度を図タイトルにする
plt.show() # 図を表示する; 複数枚図を描くとき必要
drawhov_sst(80)
drawhov_sst(120)
北緯40度では白塗りの陸地を挟んで3つの海が見えますが、これらは黄海・日本海・太平洋です。
北緯30度は九州の南あたりで、東経130度あたりに黒潮が通っているのが見えます。
気候値からの偏差¶
上の図では、やはりそのままのSSTは夏に温かく冬に冷たいという季節変動が卓越しており、それ以外の変動は目立たなくなってしまいます。
そこで、気候値sst_climからの偏差のホフメラー図を描く関数を定義しましょう。
def drawhov_ssta(n):
lev = np.linspace(-5,5,21)
plt.figure(figsize=(6,8)) # 新しいキャンバスを立ち上げる; 見やすくするため、ちょっと縦長にする
plt.contourf(Lon, Time, ssta[:,n,:], # 緯度nについて、色塗り等値線図を描く
levels=lev, # 色分け0.5度ごとに行う
cmap="seismic" # 色調に"seismic"を用いる
)
plt.colorbar() # カラーバーを図の右側に表示する
plt.title("latitude="+str(lat[n])) # 緯度を図タイトルにする
plt.show() # 図を表示する; 複数枚図を描くとき必要
練習4-6の解答例¶
# 解答例1: 月の通し番号についてforループを回し、各月の気候値を引く
ssta = np.zeros((180, 161, 241)) # 先に入れ物を定義する
for n in range(180):
ssta[n,:,:] = sst[n,:,:] - sst_clim[n%12,:,:] # 月の通し番号を12で割った余り -> 0, 1, 2, ..., 11, 0, 1, 2,...
# 解答例2: 1~12月についてforループを回し、それぞれの月の気候値を引く
ssta = np.zeros((180, 161, 241))
for n in range(12):
ssta[n::12,:,:] = sst[n::12,:,:] - sst_clim[n,:,:] # nから始まり12とばしで終わりまで -> 0, 12, 24, ..., 1, 13, 25, ...
# 解答例3: reshapeで配列の時間軸を15x12の形状に折り畳み、ブロードキャスト（後述）を用いて引き算し、元の形状にもどす
ssta = sst.reshape((15, 12, 161, 241)) - sst_clim.reshape((1, 12, 161, 241))
ssta = ssta.reshape((180, 161, 241))
データの例 (1) 北緯32.5度¶
ここでは興味深いデータの見られた2つの緯度（北緯32.5度、北緯27.5度）の図を見てみましょう。
drawhov_ssta(110)
まず北緯32.5度は四国よりわずかに南に位置する緯度で、東経130度付近の白い陸地が九州です。
2018年の少し前から、東経135度あたりに強い負の温度偏差（冷水）が見られます。これは現在（2022年12月）も続く黒潮大蛇行の期間であり、蛇行の北側にある冷水渦を横切っているために見られているものです。
データの例 (2) 北緯27.5度¶
次はより南の北緯27.5度でのホフメラー図を見てみましょう。
drawhov_ssta(130)
真横の縞状模様に加え、左上がりの模様が見えるでしょうか。真横の線は全域が同時に暖かい・冷たいという状況を表しており、これはその月が特別暖かい・冷たいというイベントをとらえています。
一方、左上がりの線は温暖・寒冷のパターンが時間と共に西に伝わってゆくことを表しています。暖水・冷水を中心とした多数の渦が、約5cm/sの速度（ホフメラー図から確かめてみましょう、経度1度は約100km）で西向きに伝播している様子が見て取れます。
これは地球が球面であるために自転効果が緯度変化することから生じるロスビー波の性質です。
ホフメラー図はこのような「シグナルの伝播」といったパターンをとらえるのに非常に有用です。
a = np.arange(6).reshape((2,3))
b = np.array([0,10,20]).reshape((1,3))
print(a)
print(b)
aとbは形が異なるので原理的にはベクトル演算はできませんが、ブロードキャストのおかげで次のようにベクトル演算が行えます。
print(a+b)
print(a*b)
何が起こっているかわかりますか。
形状(2,3)のaと形状(1,3)のbを計算しようとしたとき、bは第1次元に関して繰り返した(2,3)の配列、つまり次のようなものとして解釈されているのです。
print(np.array([[0,10,20], [0,10,20]]))
ブロードキャストを行うには、基本的には次元数ndim（今の場合2）、要素数が1ではない次元の要素数が一致していればよいです。
そのため、(2,3)の配列に対しブロードキャストできるのは(2,3), (2,1), (1,3), (1,1)の4種類ということになります。
# (2,3)と(2,1)のブロードキャストの例
c = np.array([-5, 10]).reshape((2,1))
print(a)
print(c)
print(a+c)
2次元グラフを描く¶
# 例：2次元のグラフを描く
x = np.linspace(-5, 5, 201).reshape((1, 201))
y = np.linspace(-5, 5, 201).reshape((201, 1))
z = np.sin(x) * np.exp(-0.5*(x**2 + y**2))
X, Y = np.meshgrid(x, y) # これはどっちにしても必要
plt.contourf(X, Y, z, 41)
plt.colorbar()
z = np.sin(x) * np.exp(y)
plt.contourf(X, Y, z, 41)
plt.colorbar()
練習4-7¶
上と同様にx, yのブロードキャストを用いて、次の関数の2次元グラフを描きなさい。
- $z=xy$
- $z=\frac{x^2}{3}+\frac{y^2}{6}$
- $z=\sin\left(\frac{x^2-y^2}{2}\right)$