X Window (Windows パソコンでは Xming) を起動したうえで、 端末 (Windows パソコン では TeraTerm) で niigata にログインして、 各自のホームディレクトリで (この授業用のディレクトリをつくっている人はそこに cd したうえで) つぎのことをやってください。
「.csv」でおわる名まえのファイルが 1つと、 「.py」でおわる名まえのファイルが 4つあるでしょう。
これは、観測値にもとづいてつくられたデータセットだが、直接の観測値ではない。 陸上は地上気温、海上は海面水温の観測値にもとづき、 長期平均からの偏差を推計したものである。 長期平均としては、1991-2020年の平均が採用されている。(対象全期間でないことに注意。)
ここでは、気象庁が作成した csv ファイルを Linux (niigata) 上で読んで、気象庁ウェブサイトにあるグラフを再現してみることにする。
グラフ作成に matplotlib (のうちの pyplot) をつかう。 (「import matplotlib.pyplot as plt」で予約している。) 今回も、前回にひきつづき、fig と ax をわける方法でなく、単純に plt の機能を呼びだす方法でつかう。
プログラム例 41 では、 ひとまず、最近11年ぶんの値を人が読んで、プログラムに書きこんでおいた。 リスト変数「year」が年、「tmp」が全球平均・年平均の地表温度である。
matplotlib.pyplot の plot 関数では、第1, 第2の位置引き数であたえたものを、それぞれ、横軸の座標値、 縦軸の座標値とみる。
つぎの3つのうち 1つを生かしてほかはコメントにしてある。 コメントの印の「#」(ハッシュマーク。音楽の「シャープ」とは厳密には同じではないが同一視する人もいる) を消したり入れたりしてためしてみてほしい。
x軸のラベルとして 'year' を入れ、y軸のラベルは省略し、 図の上の表題として、データが何であるかを (日本語ローマ字でもよいのだが、例として) 英語で説明している。 (日本語の漢字・かなのあつかいは (OS によってかわるところがあるので) この授業ではあつかわない。)
このまま実行すると図は X Window の画面表示になるが、 plt.show をコメントにして plt.savefig をふくむ4行を生かすと 画像ファイル出力になる。 この点はこのあとの例でも同様。 今回の実習では、44 まで行ってから画像ファイル出力にしよう。
csvファイルのすべてではないが、縦の列と横の行からなる表の形にととのった内容のものは、 Python の pandas というパッケージ の read_csv という関数で データフレーム (DataFrame) という型の変数に読みこんであつかうのが便利である。
プログラム例 42 では、an_wld_d20240601.csv を read_csv で dftmp というデータフレームに読みこんで、41 と同じ手順でグラフにする。
気象庁からダウンロードした csv ファイルには、 最初の行に、各列の表題としてつかえる文字列がふくまれているが、 日本語 (Shift-JIS 文字コード) である。そのままではあつかいにくいので、 ここでは read_csv のオプション「skiprows=1, header=None」で最初の行を無視し、 names で列の表題を英語の略語であたえた。 gl [ジー エル] は global の この場かぎりの略語 [注], nh は Northern Hemisphere, sh は Southern Hemisphere の一般的略語である。 year, gl, nh, sh の 4つの列からなるデータフレームができる。
このプログラム例では、year を 横軸、gl を縦軸にとったグラフをかく。 gl を nh や sh にかえれば、北半球平均や南半球平均の地表温度のグラフができる。 そのばあいは、plt.title で指定する文字列の「global」を「NH」や「SH」に、 また、画像ファイル出力をするばあいはファイル名中の「global」も「nh」や「sh」に変えておこう。 (このような頭文字略語は目で見るときは大文字がよい。 しかしファイル名は小文字で統一したほうがあつかいやすいとわたしは思う。)
速い (周期の短い) 変動をふくむ時系列から 遅い (周期の長い) 変動をとりだすために、 移動平均 (英語では running mean あるいは moving average) という方法がよくつかわれる。 (同様な目的に、ディジタルフィルター (digital filter) のうちのローパスフィルター (low pass filter) がつかわれることもある。 このばあいの low pass は、振動数 (周波数) が低い成分をとおす、という意味である。 移動平均はローパスフィルターの一種だともいえる。)
移動平均は、平均をとる区間の幅を一定にしながら 時系列のうちで区間をうごかしていって、 得られた平均値をならべた時系列を作る操作である。
プログラム例 43 では、気象庁の図と同じように、5年移動平均をつくる。 そのばあいの入力データの項数 5 を、変数 stride の値であたえるようにした。 この値を、たとえば 11 にかえれば、11年移動平均を計算できる。 (移動平均をとる区間の中央をデータがあたえられた時刻のひとつと一致させたいので、 stride の値は奇数であることを想定している。)
データ期間の両端からそれぞれ (stride-1)÷2 個 (5年移動平均ならば2年ぶん) の値は さだまらない。プログラムでは、Numpyパッケージで「数値ではない」ことを意味する np.nan がはいるようにしてある。 (pandas で「有効な値ではない」ことを意味する pd.NA が適切だと思ったのだが、 matplotlib の plot でエラーになってしまった。)
pandas パッケージの rolling という機能をつかうと、 つぎの1行で、データフレーム dftmp の列 gl の移動平均を計算して 列 ma に入れることができる。
dftmp['ma'] = dftmp['gl'].rolling(stride, center=True).mean()
rolling にたよらず、自分で移動平均を計算して、結果を列 ma2 に入れるプログラムもつくった。
ma と ma2 の結果は、(小数点下3けたを書かせてくらべたかぎりでは) 一致していた。 したがって現状のプログラムでは ma2 はグラフ表示していない。 実用的には ma を計算する部分と ma2 を計算する部分はどちらか一方を残して他方を消すのがよい。
なお、神山 (2024) の本では 第2巻の第3章で移動平均をあつかっていて、 Numpy パッケージの convolve という関数をつかっている。 (この本では pandas をつかっていない。)
グラフでは移動平均を青の折れ線でしめした。 折れ線の意味を説明する凡例 (legend) もつけた。 凡例にかく文字列は、plot のなかで「label」で指定している。 (label が空文字列 (長さ0の文字列) だと、 その項目の凡例が省略されるようになっているようだ。)
これも 42 と同様に gl を nh や sh に変更することが可能にしてある。
時系列全体に、回帰 (regression) という方法で、直線をあてはめることもよくおこなわれる。 直線回帰は、x から y をもとめる推定式「y = a x + b」をたて、 y の推定値と観測値の差の2乗を最小にするように係数 a, b の値をもとめることにあたる。 ここでは x を時間、y を温度とする。
直線回帰は、多項式回帰の特殊なばあいとして、Numpy パッケージの polyfit 関数で 係数 a, b を計算することができる。 (この方法は、神山 (2024) の 第1巻 第5章「線型トレンドとその除去」にならった。
回帰の定義にさかのぼれば、a, b は、x と y の平均・分散・共分散から 計算することができる。(ここでは dftmp の "year" の列を x、"gl" の列を y とする。) 林 (1985) の第4章「相関について」の第1節「相関係数について」の 記述を参考にして、係数を aa, bb として得るプログラムを書いた。
a, b と aa, bb はじゅうぶんな精度で一致しているので、実用にはどちらか一方でよい。
43 の移動平均を表示するプログラムに追加して、 気象庁の図と同様に赤で、回帰直線を表示するようにした。 凡例のところに、気象庁の図と同様に、a の値を表示した。b の値は省略している。