はじめに¶
前回は描画ライブラリmatplotlibの力を借りて、時系列データの折れ線グラフを描きました。 描画のために1年間・日ごとの気温データなどをCSV形式から読み込み、配列として保存しましたが、その際に用いたのがNumPyモジュールのloadtxt関数でした。
NumPyは、上記のような多数の数値データを扱うのに長けた機能が集まったライブラリで、Pythonを用いた科学計算の中核をなすものです。中でも、ndarrayと呼ばれる配列の型（リストを機能拡張したものと思えばよいでしょう）は1次元的・多次元的に敷き詰められた数値データを高速に・簡易に扱う機能が詰まっています。 今回は、配列を中心としたNumPyの機能の基本を学びましょう。
慣習として、NumPyはnpという名前でインポートすることが多いです。後で使うので、前回と同様にmatplotlibもインポートしておきます。
import numpy as np
import matplotlib.pyplot as plt
print(np.pi) # 円周率
print(np.e) # ネイピア数
print(np.inf) # 正の無限大
print(np.nan) # Not a number = 数ではない
初等関数¶
三角関数や指数関数、それらの逆関数などの初等関数ももちろん入っています。
非常によく使いますので、np.とセットで列記します。
- 三角関数:
np.sin,np.cos,np.tan - 指数・対数関数:
np.exp(exponential=指数関数),np.log(自然対数),np.log10(底10),np.log2(底2) - 逆三角関数:
np.arcsin,np.arccos,np.arctan,np.arctan2(y/xではなくy,xを引数にとり、180度の不定性を除去する) - 双曲線関数・逆双曲線関数:
np.sinh,np.cosh,np.tanh,np.arcsinh,np.arccosh,np.arctanh - 平方根:
np.sqrt(square root=平方根)、ほかの冪関数はa**bまたはnp.power(a,b)で表記する
print(np.sin(np.pi/4)) # sin, cos, tanはラジアン表記
print(np.sqrt(1/2)) # ルート2分の1, sin(pi/4)と等しいはず
print(np.exp(2))
print(np.sqrt(-1)) # マイナスの数の平方根は、虚数ではなくnanになります
print(np.log(-1)) # マイナスの数の対数も、nan
print(np.exp(1j*np.pi)) # 複素数を引数にもできる。e^(i*pi)=?
数を扱うその他の関数¶
- 絶対値:
np.abs（実数・複素数どちらにも使える） - その数以下の最大整数:
np.floor - その数より大きい最小整数:
np.ceil - 複素数の実部・虚部・角度（ラジアン）・複素共役:
np.real,np.imag,np.angle,np.conj
print(np.abs(4))
print(np.abs(-3))
print(np.floor(1.5), np.ceil(1.5))
print(np.floor(-1.5), np.ceil(-1.5))
統計¶
リストや（次に説明する）ndarrayなどの配列を対象にとって、合計・平均・分散・中央値などの統計量を計算する関数を紹介します。
# 統計の関数の例
import numpy as np
x = [6, 70, 39, 16, 56, 36, 19, 9, 24, 40, 22, 91, 80, 78, 84, 7, 29, 90, 86,
12, 93, 77, 32, 71, 95, 69, 80, 54, 14, 27, 74, 28, 21, 30, 92, 83, 96,
58, 51, 22, 22, 75, 81, 98, 57, 41, 13, 1, 48, 90]
print(np.sum(x)) # 合計
print(np.mean(x)) # 平均
print(np.var(x)) # 分散: 自由度Nだがオプション"ddof=1"で自由度N-1にもできる
print(np.std(x)) # 標準偏差: 自由度は分散と同様
print(np.median(x)) # 中央値
print(np.max(x)) # 最大値
print(np.min(x)) # 最小値
配列にNaNが混じっている場合、統計値はNaNを返します。NaNを除いた統計値を返す関数として、np.nansum, np.nanmean, np.nanmax, np.nanmin...などがあります。
NumPyの配列¶
Numpyの最大の長所が、配列の型であるndarrayオブジェクトです。 C言語やFortranの「配列」に相当する型で、Python固有のリストにも似ています。
ndarrayとリストの特徴を比べると、以下のようにまとめられるでしょう。
- ndarray
- 予め大きさ（要素数）を決めておく
- 要素の型は統一（浮動小数点数、整数、複素数、論理型など）
- まとめて計算する際のベクトル演算が非常に高速
- 多次元データの操作性に長ける
- リスト
- 要素数を付け足したり減らしたり柔軟に変更できる
- 異なる型の要素も含めることができる
- まとめて計算する際にはforループなどを使用する必要があり、遅い
- データは原則1次元構造
このように、時系列データやシミュレーション結果など、膨大な量の数値の羅列を扱う際には極めて強力な機能が実装されています。ここでは、配列の生成方法やベクトル演算・メソッドなどの配列操作の基本を解説します。
datalist = [1,3,5,7,9]
print(type(datalist))
print(datalist)
dataarray = np.array(datalist)
print(type(dataarray))
print(dataarray)
2. 外部データを読み込む¶
前回KEOの時系列データをnp.loadtxt関数で読み込んだように、CSVファイルなどをプログラム内に読み込みます。
後に、「データ入出力」で解説します。
3. 生成関数を用いる¶
NumPyには、ある決まった規則で配列を生成する関数があります。ここでは、よく用いるものを紹介します。
細かい使い方は、関数名?と入力するか、公式リファレンスを確認しましょう。
a = np.zeros(5) # 要素数5の、0ばかりの配列をつくる
print(a)
b = np.ones(10) # 要素数10の、1ばかりの配列をつくる
print(b)
c = np.arange(50) # rangeと同じく、start=0, end=50, step=1の等差数列
print(c)
d = np.arange(1, 50, 3) # start=1, end=50, step=3
print(d)
e = np.linspace(-2, 2, 11) # -2から2までを**両端点を含めて**等分した等差数列
print(e)
# d = [1, 4, 7, ..., 49], 17要素
f = np.zeros_like(d) # dと同じ大きさの、0ばかりの配列
print(f)
g = np.ones_like(d) # dと同じ大きさの、1ばかりの配列
print(g)
練習3-1¶
以下の配列を生成しなさい。生成関数を用いること。
- [0, 1, 2, 3, 4, 5]
- [0, 0, 0, 0, 0, 0]
- $[0, 0.2\pi, 0.4\pi, ..., 1.8\pi, 2\pi]$
要素の取り出し・スライス¶
ndarrayはリストと同様に、配列名[要素番号]や配列名[start:end]、配列名[start:end:interval]のような形で要素を取り出せます。
print(d)
print(d[2]) # 0から数え始めて2番目
print(d[-2]) # 後ろから2番目
print(d[5:8]) # 5から7番目
print(d[1::2]) # 1番目から最後まで、2個間隔で
a = np.arange(100) # [0, 1, 2, ..., 99]
b = np.arange(100, 0, -1) # [100, 99, 98, ..., 1]
c = np.zeros(100) # 入れ物を先に用意する
for i in range(100):
c[i] = a[i] + b[i]
print(c)
上と同じ計算を次のように一気に計算できます。 ベクトルの成分ごとの足し算$c_i=a_i+b_i (i=0,...,99)$に対して、ベクトル丸ごとの足し算$\textbf{c}=\textbf{a}+\textbf{b}$を行うイメージで、ベクトル演算と呼ばれます。
a = np.arange(100) # [0, 1, 2, ..., 99]
b = np.arange(100, 0, -1) # [100, 99, 98, ..., 1]
c = a + b # [100, 100, ..., 100]のndarray
Pythonのforループは処理が非常に遅いので、可能な限りベクトル演算を使うことが重要です。
a = np.arange(1000000) # 非常に大きい配列
b = np.roll(a, 500000) # aと同じ大きさの、別の配列
print(a)
print(b)
%%timeit
# 時間計測
c = np.zeros(1000000)
for i in range(1000000):
c[i] = a[i] + b[i]
%%timeit
d = a + b
関数を作用させる¶
NumPy配列は配列丸ごとの四則演算が行えるほか、np.sinやnp.expといった関数を作用させることも可能です。
実行したい計算をできるだけループを使わずにベクトル演算の形式で表記できるかが、Pythonプログラムの効率化に重要です。
# ベクトル演算の例
a = np.linspace(0., 1, 11) # 0 ~ 1..., 10等分した端点11個
b = a + 1 # 成分それぞれに1を足す
print(a)
print(b)
print(a * b) # 成分ごとに掛け算
print(2 ** b) # 第i成分について、2のb[i]乗
print(np.cos(2*np.pi*a)) # 成分ごとにコサインを計算する
print(np.exp(a)) # eのa[i]乗
print(np.log(b)) # log(b[i])
print(np.sqrt(a+b))
練習3-2¶
以下の配列を作りなさい。2.以降ではベクトル演算形式を用いること。
- -2, -1.9, ..., 1.9, 2の41要素をもつ配列（
xの名前で定義すること） xの2乗- 2の
x乗 - $e^{-x^2}$
y = np.linspace(-10, 10, 101) # -10から10までの100等分
print(y.size) # 配列のサイズ
print(y.dtype) # 格納されている数値の型: 今はfloat64=8バイトの浮動小数点数
メソッド¶
ndarray型の変数には属性に加えて、その配列を対象とする操作を行うための関数も含まれています。
そういった、変数に付属した関数をメソッドと呼びます。属性と似た形で、変数名.メソッド名(引数)でアクセスできます。
ndarrayには、合計・平均・分散・標準偏差・最大最小といった統計値を計算するメソッドがあるのが非常に便利です。
勿論np.meanなどの関数を呼んでもよいのですが、関数の対象を省略できるのでコードの短縮につながります。
np.median, np.nanmeanなどの関数に対応するメソッドはないので、その場合はおとなしく関数を呼びましょう。
print(x) # 上の「統計」で定義した、ランダムな長いリスト
x = np.array(x) # リストをndarrayに変換
print(x.mean()) # 平均
print(x.max()) # 最大値
print(x - x.mean()) # xからxの平均を引いた配列
実例：KEOデータの解析¶
配列の基本を学んだところで、実例形式でより多くの使い方を学んでいきましょう。
airt = np.loadtxt(r"C:\Users\yfuji\Desktop\lesson\data\data_KEO\KEOairt_2019daily.csv")
sst = np.loadtxt(r"C:\Users\yfuji\Desktop\lesson\data\data_KEO\KEOsst_2019daily.csv")
wind = np.loadtxt(r"C:\Users\yfuji\Desktop\lesson\data\data_KEO\KEOwind_2019daily.csv")
print(airt)
バイナリファイルの入出力¶
CSVファイルはテキストファイル（人間が読める、文字列の羅列）なので、数字一桁（0~9の10種類）を記録するのに1バイト必要です。 しかし、本来1バイト=8ビットの情報量では、2の8乗=256種類の数値が記録できるはずなので、テキストファイルは記録効率が悪いと言えるでしょう。
このように、（テキストとしては読めないものの）高効率なデータ保存形式をバイナリファイルと呼びます。大気・海洋の大規模データは記録容量の効率化のためバイナリファイルが用いられることがほとんどです。バイナリの中にも様々な保存形式があり、中でもNetCDF形式が良く使われますが、これは別の機会に説明することにします。
NumPyにはいろいろなバイナリ形式を読む関数が含まれているほか、NumPy独自の形式で配列をバイナリファイルとして出力・読み込みできる関数np.save, np.loadを持っています。これについて紹介します。
# 保存するファイルのパス文字列。拡張子".npy"を用いる
filepath = r"C:\Users\yfuji\Desktop\lesson\data\data_KEO\KEOairt_2019daily.npy"
np.save(filepath, airt) # 配列airtをfilepathに保存
%ls C:\Users\yfuji\Desktop\lesson\data\data_KEO\KEOairt_2019daily.*
上ではJupyterlabの機能で、KEOairt_2019daily.○○というファイルのサイズを羅列しています。
2190, 3048というのがファイルサイズで、一見.npyのほうがファイルが大きくなっているように見えますが、CSVのほうは十の位～小数点以下1桁までの有効数字3桁で記録しているのに対し、.npyのデフォルトが8バイト実数のバイナリ形式なので有効数字約16桁でデータが記録されています。このように、データ量に対する容量の効率が飛躍的に向上しています。
実際の気温測定精度を考えると有効数字16桁で記録するのは無駄なので、4バイト実数（有効数字8桁）・2バイト実数（同4桁）あたりを使うのが普通です。
では、この.npyファイルを読み込んでみましょう。
保存前のairtと同じ内容が読み込めているか、論理演算をベクトル形式で行って確かめてみます。
airt_new = np.load(filepath) # filepathは上で定義したパス文字列
print(airt_new)
print(airt_new == airt)
また、np.saveで生成される.npyファイルは一つの配列を保存しますが、np.savez関数を用いると複数の配列を.npzという拡張子の1つのファイルに保存することができます。自分で調べて、使ってみましょう。
データの概観¶
今回のように大きなデータでは、個々の数字を見ていてもその全体像がよくわかりません。ざっとその全体像をつかむには、統計量を見たりグラフを描くのが良いでしょう。matplotlibとndarrayのおさらいです。
time = np.arange(len(airt))
plt.plot(time, airt, label="airt")
plt.plot(time, sst, label="SST")
plt.legend()
plt.plot(time, wind, label="wind speed")
plt.legend()
print("Air temp:", airt.mean(), airt.max(), airt.min()) # 平均、最大、最小
print("SST:", sst.mean(), sst.max(), sst.min()) # 平均、最大、最小
print("Wind speed:", wind.mean(), wind.max(), wind.min()) # 平均、最大、最小
マスク¶
ここまでは1年間を通したデータを見てきましたが、夏季のみ・冬季のみや、風速が一定値以上の場合のみなど特定の条件下のデータを取り出して、その条件下での平均的描像を調べることを考えます。そうした解析手法をコンポジット（合成）解析と呼びます。
そういったときに便利なのがマスクです。 ndarrayは自身と同じ大きさで、論理型 (True/False) が入った配列を要素番号の代わりに用いることで、Trueの要素だけを抜き出すことができるのです。
# 例
a = np.array([1,2,3,4,5])
mask = np.array([True, False, False, True, True])
print(a)
print(mask)
print(a[mask]) # Trueの要素だけが抜き出されている！
上の例ではマスクを手動で作りましたが、マスクはベクトル形式の論理演算で容易に作成できます。
ただし、「かつ」「または」「否定」はand or notではなく、以下のような書き方をする必要があります。
xかつy：np.logical_and(x, y)、x & y、x * yxまたはy：np.logical_or(x, y)、x | y、x + yxではない：np.logical_not(x)、~x
また、要件を満たしてTrueになった個数を数えたいときにはmask.sum()のようにすれば、Trueを1, Falseを0として単純に足し算をしてくれます。
a = np.arange(30) # 0~29
mask2 = (a % 2 == 0) # 2で割った余りがゼロ = 偶数
mask3 = (a % 3 == 0) # 3の倍数
print(mask2) # 偶数マスクの例
print(a[mask2])
print(a[mask3])
print(a[mask2 & mask3]) # 偶数かつ3の倍数=6の倍数
print(a[mask2 | mask3]) # 偶数または3の倍数
print((mask2|mask3).sum()) # 偶数または3の倍数の個数
さて、KEOの時系列に戻りましょう。 大気・海洋でよく使う定義として、夏季を6,7,8月、冬季を12,1,2月とします。英語の月名の頭文字から、前者をJJA、後者をDJFと表記します。論文でも説明なしに使われるので、覚えておきましょう。
配列timeは1月1日を0として以下1日ずつカウントしていますが、6月1日は151日目、9月1日は243日目ですので、JJAのマスクは(time>=151) & (time<243)として作れます。
同様に、3月1日は59日目、12月1日は334日目ですので、DJFのマスクは(time<59) | (time>=334)となります。
ここではJJAの平均風速、DJFの平均風速を調べてみましょう。
mask_jja = (time >= 151) & (time < 243)
mask_djf = (time < 59) | (time >= 334)
print(wind[mask_jja].mean())
print(wind[mask_djf].mean())
2次量の算出¶
ここまで気温・水温・風速だけを見てきましたが、これらから算出される量もベクトル形式で簡単に算出できます。
例として、顕熱フラックスを計算してみましょう。フラックスとは、単位時間あたりに、ある物理量が単位面積の面を通過する量を表し、ここでは特に1平方メートルの海面を1秒間に通過する熱（単位は$\left[\frac{{\rm J}}{{\rm m}^2\ {\rm s}}\right]=[{\rm W}/{\rm m}^2]$）のことを表します。
海面の顕熱フラックスは海面付近の大気・海洋境界層乱流によって決まるため厳密な計測は困難ですが、以下のような形で大まかに近似されます（海洋分野ではバルク法と呼びます）。空気の密度を$\rho_a\approx1.3\ [{\rm kg}/{\rm m}^3]$、空気の定圧比熱を$C_p\approx1004\ [{\rm J}/({\rm K}\ {\rm kg})]$、顕熱フラックスを$H\ [{\rm W}/{\rm m}^2]$と表し、大気から海洋への熱の移動を正とします。
$H=\rho_aC_pC_hU_{10}\Delta T,\ \Delta T=T_{\rm air}-{\rm SST}$
ここで$C_h$はバルク係数と呼ばれます。バルク係数は観測やモデルに基づき様々な値が提案されており、風速に依存して変化することが示唆されていますが、ここではざっくりと$C_h=1.2\times 10^{-3}$を使うことにしましょう。
$U_{10}$は高度10mでの水平風速です。KEOの風速計は高度4mに設置されているので、風速を10m高に補正してやらないといけません。ここでは非常におおざっぱに、風速が高度の1/7乗に比例する、という経験則をもちいて補正することにしましょう。
$\left(\frac{U_{10}}{U_{4}}\right) = \left(\frac{10}{4}\right)^{1/7},\ U_{10} = (10/4)^{1/7} \times U_4$
rhoa = 1.3 # 空気の密度
Cp = 1004 # 空気の定圧比熱
Ch = 1.2e-3 # バルク係数
u10 = (10/4)**(1/7) * wind
H = rhoa * Cp * Ch * u10 * (airt - sst) # ベクトル形式でフラックスを算出
plt.plot(time, H, label="sensible heat flux")
plt.legend()
print("all", H.mean())
print("summer", H[mask_jja].mean())
print("winter", H[mask_djf].mean())
これらの解析から、以下のようなことがわかります。
- 平均的に、熱フラックスは負 = 海から大気に熱が放出
- 夏季平均ではほとんどゼロ、わずかに大気から海に熱が渡されている
- 冬季は海から多量の熱が大気に放出されている
積分と条件付き積分¶
KEOでは平均的に海から大気へと熱が放出されていることがわかりました。 第$k$日目に$1{\rm m}^2$あたり1日に放出される熱量は$H_k\Delta t\ [{\rm J}/{\rm m}^2]$（$\Delta t=1{\rm day}=86400{\rm s}$）なので、年間では$\sum_{k=0}^{364}H_k\Delta t\ [{\rm J}/{\rm m}^2]$の熱量が放出されることになります。 これは、$\int_0^{365} H dt$ （時間の単位は日）の近似ですね。
dt = 24 * 60 * 60 # 1日の秒数
print((H * dt).sum())
このように、2019年にKEOでは$6.23\times10^8\ [{\rm J}/{\rm m^2}]$の熱が海から大気に放出されたことになります。
さて、熱の放出は風が強いほど多くなります。前回見たように、風は低気圧や台風の通過に伴い強弱を繰り返します。 そこで「年間の熱放出のうちの強風イベントの寄与」を解析してみましょう。具体的には、「風速$U_0$[m/s]以上の日による熱放出が、全熱放出の何%を占めるか？」ということを調べます。 数式で書くと、次の通りです。
$\left.{\int_{U\geq U_0}H dt}\right/{\int_{\rm all}H dt}$
いろんな風速について調べたいので、$U_0$を引数とする関数で書きましょう。マスクを使えば、条件付き積分も簡単にできます。
def heatflux_fraction(u0):
dt = 24 * 60 * 60
num = (H[wind >= u0] * dt).sum() # 分子 (numerator)
denom = (H * dt).sum() # 分母 (denominator)
return num / denom
print(heatflux_fraction(8))
print(heatflux_fraction(10))
print(heatflux_fraction(12))
これによって、風速8m/s以上の日による熱放出が全体の54%, 10m/s以上による寄与が26%, 12m/s以上による寄与が11%であることがわかりました。ですが、それぞれの風速がどれほどの頻度で起こるかの情報がないと、まれな強風イベントがどれほど全体に重要なのかきちんとわかったことにはなりません。そこで、上と同様に「風速$U_0$[m/s]以上の日がどれほどの頻度で起こるのか」を計算する関数を作りましょう。
def wind_fraction(u0):
return (wind >= u0).sum() / len(wind) # 365分の〇日
print(wind_fraction(8))
print(wind_fraction(10))
print(wind_fraction(12))
12m/s以上を例にとると、上位4%の強風の日による熱放出が全体の11%の熱放出を占めた、ということになります。 実社会では世界の1%の富裕層が全資産の37%を保有する（ソース）そうですが、それに比べると海からの熱放出は極端イベントの相対的な重要性が低いですね。
横軸に風速の下位からの割合、縦軸にそれらの日に放出された熱の全体に対する割合を描くと、次のようになります。
u_array = np.linspace(0., 20., 100) # 調べる風速の列
windfrac_array = np.zeros_like(u_array)
heatfrac_array = np.zeros_like(u_array)
for n in range(len(u_array)):
# u_arrayの風速に対応する割合を計算していく
windfrac_array[n] = 1 - wind_fraction(u_array[n]) # 与えた風速より*弱い*日の割合
heatfrac_array[n] = 1 - heatflux_fraction(u_array[n]) # 与えた風速より*弱い*日の熱放出の寄与
plt.figure()
plt.plot(windfrac_array, heatfrac_array)
plt.plot([0,1], [0,1], 'k--') # x=yの線
plt.grid() # 格子を描く
plt.show()
このような図をQ-Qプロットと呼びます。 例えば線が(0.8, 0.6)の点を通っていますが、これは上位20%の強風が全体の熱放出の40%を担っている、と見ることができます。
その他の便利な関数¶
その他NumPyの便利な関数を、ざっと紹介します。人によっては必要ない機能も多いと思いますが、必要になった時に「こんな機能あったな」と覚えていてもらえたら十分です。
- ベクトル内積
np.dot - 線形連立方程式の解
np.linalg.solve - 最小二乗法
np.linalg.lstsq - 行列の固有値分解
np.linalg.eig,np.linalg.eigh - 乱数生成
np.random.rand,np.random.randn,np.random.randint - 離散フーリエ変換
np.fft.fft,np.fft.ifft,np.fft.rfft,np.fft.irfft - ソート・検索
np.sort,np.argwhere,np.argmax,np.argmin - 線形内挿
np.interp