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