[気象庁の観測データを使う (序論)] (板倉用教材) のページで、 気象庁ウェブサイトの「過去の気象データ ダウンロード」の機能をつかって、 例として、館林の2年ぶんの日ごとの値をダウンロードした。 (2020年12月13日にダウンロードしたものを、つぎのところにも置いておく。)
[気象観測データのグラフ作成と基本的な統計処理 (Rによる) (1)] (板倉用教材) のページで、 そのデータを R で 読んで、図化や基本的な統計処理をした。
ここでは、同じデータを Python で読んで、図化や基本的な統計処理をこころみる。 【説明は Windows 上 の Anaconda にくみこまれた Python 3.8.5 での 実行結果によるが、他の OS でも同様にできるはずである。このページの例文の実行には Anaconda は必要ではない。】 【Anaconda は便利だが、ディスクを大量に消費するので (標準的インストールをしただけで 3 ギガバイト)、 最小版の Miniconda を入れて、必要に応じてパッケージを追加していったほうがよいかもしれない。】
作業用ディレクトリをつくっておく。 気象庁からダウンロードしたCSVファイルを、 ファイル名を「daily_tatebayashi.csv」として、作業用ディレクトリに置く。 (Windows ならば Anaconda Powershell Propmt または Anaconda Prompt の中で) cd コマンドを実行して、 作業用ディレクトリをえらんでおく。
Python 言語のスクリプトを、テキストエディタで作成する。 ファイル名は「.py」でおわるようにする。たとえば、「ex.py」としよう。 これを(Windows ならば Anaconda Powershell Propmt または Anaconda Prompt の中で) つぎのようにして実行する。
python ex.py
【 [2021-05-05 補足] Anaconda (わたしがインストールしたのは Anaconda3-2020.11-Windows-x86_64.exe ) を インストールすると、Spyder というプログラムがふくまれている。これには、 コンソール窓、エディタ窓、plot 窓 など複数の窓があって、 Python プログラムを編集しながら対話型で実行することができる。 すすめられると思うが、わたしはまだあまり経験がないので紹介を書けない。】
これからつくる Python 言語のスクリプトの中では、あらかじめ、つぎのように、 呼び出すパッケージを予約し、略号を指定しておく。
import pandas as pd import numpy as np import matplotlib.pyplot as plt
Python のパッケージ Pandas (Anaconda には常備されている) には「データフレーム」というデータ構造がある。 R のデータフレームとよくにている。 おおざっぱには、縦横の表と思えばよい。 ただし、行はおたがいに同質のデータだが、列はそれぞれにちがった意味をもつことがある。
データフレームには名まえが必要で、名まえは自由につけられるが、 ここでは、だいたい1つだけのデータフレーム (および、それから派生したもの) をつかうので、 データフレーム名を(データフレームであることをわすれないように)「df」としてみる。
Pandas には、CSVファイルをデータフレームに読みこむ「read_csv」という関数がある。 (関数名が下線文字(「_」, underscore)をふくんでいることに注意。) これをつかって、CSVファイル「daily_tatebayashi.csv」の内容をデータフレーム「df」に読みこめる。
ダウンロードしたファイルには、頭から6行に説明があり、7行めから実質のデータがはじまる。 各列の項目の説明もあるが、そのままデータフレームの列名にするのには適していないので、 6行を全部読みとばし(skiprows=6)、ファイルの頭の部分から変数名を読みとらない(header=None)、という選択をした。 また、ファイルの文字コードが Shift_JIS であることを指定した。
【Pythonの変数名は ASCII コードであらわせるアルファベット・数字・ 下線文字 (「_」、アンダースコア) で構成する。はじめの字はアルファベットにする。 (下線文字ではじめることも可能であり、 わたしのプログラム例では、ある種類の定型処理で、参考にした例文の表現をひきついで、下線文字だけの変数名をつかうことがある。) 大文字と小文字は区別されるが、人がまちがえやすいので、 わたしの教材では、大文字にきまっているもの以外は、小文字だけをつかうことにする。 日本語文字なども Unicode であらわされるかぎりつかえるらしいが、ここではつかわないことにする。】
列名が番号だとおぼえきれないので、実際の気象要素に対応する名まえにつけなおすことにする。 日本語ローマ字でもかまわないのだが、ここでは英語の略語にした。
各列の内容は、 たとえば df['year'] (データフレームのあとに、列名が、single quotation marks でかこまれ、さらに角かっこでかこまれている) のような形で取り出すことができる。 また、たとえば「df.year」(データフレーム名、「.」(ピリオド)、列名が、空白をはさまずにつなげられている) という形でも (Python文法上のほかの構文とまぎれなければ) 参照できる。
[気象庁の観測データを使う (序論)]のページでつくった館林の日ごとの値のばあい、 列番号 (Python のつごうで、はじまりは 0)、わたしのPythonプログラムでの列名、気象データとしての項目、 表計算ソフトウェアで読んだときの列記号は、 つぎのように対応している。
| 列番号 | 列名 | 項目 | 表の列 |
|---|---|---|---|
| 0 | year | 年 | A |
| 1 | month | 月 | B |
| 2 | day | 日 | C |
| 3 | temp | 平均気温 | D |
| 4 | q_temp | 平均気温の品質情報 | E |
| 5 | h_temp | 平均気温の均質番号 | F |
| 6 | tempmax | 最高気温 | G |
| 7 | q_tempmax | 最高気温の品質情報 | H |
| 8 | h_tempmax | 最高気温の均質番号 | I |
| 9 | tempmin | 最低気温 | J |
| 10 | q_tempmin | 最低気温の品質情報 | K |
| 11 | h_tempmin | 最低気温の均質番号 | L |
| 12 | precip | 降水量 | M |
| 13 | q_precip | 降水量の品質情報 | N |
| 14 | h_precip | 降水量の均質番号 | O |
| 15 | sundur | 日照時間 | P |
| 16 | q_sundur | 日照時間の品質情報 | Q |
| 17 | h_sundur | 日照時間の均質番号 | R |
ここまでの作業をするスクリプトはつぎのとおり。 ただし、1つの文を複数の行にわたって書くために、各行のおわりに「\」(本来はバックスラッシュだが、 日本語環境では円(yen)マークのことがある) を置いている。 ここからさきの例文では、いつも、さきにのべたパッケージの予約と、 ここでのべたデータ読みこみがあって、そのさきの作業がつづく。
# データをファイルから読みこみ、各列に名前をつける
df = pd.read_csv('daily_tatebayashi.csv', \
header=None, skiprows=6, encoding='shift_jis', \
names=('year', 'month', 'day', 'temp', 'q_temp', 'h_temp', \
'tempmax', 'q_tempmax', 'h_tempmax', 'tempmin', 'q_tempmin', 'h_tempmin', \
'precip', 'q_precip', 'h_precip', 'sundur', 'q_sundur', 'h_sundur') )
日照時間の品質情報 df.q_sundur の値を見ると、ほとんど「8」なのだが、ときどき「4」がまざっている。 その日には、(おそらく機器点検があったので) 24時間の観測値がそろっていないのだ。 df.sundur にはその日の数値がはいっているが、それは統計処理からはずしたほうがよいだろう。 そこで、 df.q_sundur の値が 8 でない日の df.sundur の値を、欠損 (NumPy のもつ定数 nan) に変えておくことにする。
気温の品質情報には、ときどき「5」がまざっている。 ここでは、気温や降水量についても、品質情報が「8」でない観測値は欠損にしておこう。
Pandas のデータフレームの機能をつかうと、ある列の内容が特定の値であればそれを nan に変えることや、 nan であれば特定の値に変えることは、簡単にできるのだが、 「ある列の内容が特定の値であるとき、他の列のその行の内容を nan に変える」手続きを、 わたしはまだ発見できていない。 そこで、手順だて型のプログラム言語に共通なループの機能をつかって、 つぎのようなプログラムを書いた。
for kt in range(len(df)):
if (df.q_sundur[kt] != 8):
df.sundur[kt] = np.nan
if (df.q_temp[kt] != 8):
df.temp[kt] = np.nan
if (df.q_tempmax[kt] != 8):
df.tempmax[kt] = np.nan
if (df.q_tempmin[kt] != 8):
df.tempmin[kt] = np.nan
if (df.q_precip[kt] != 8):
df.precip[kt] = np.nan
「!=」は「ひとしくない」という、2つの数値から論理値を得る演算子、 「=」は代入演算 (左の変数に右の計算結果を代入する)である。 (「ひとしい」は「==」とかく。) len(df) で、dfの内容の行数(各列のデータ個数)が得られる。 実際の数は、うるう年を1つふくむ2年間の日数なので、365×2+1 である。 この for ループは、変数 kt の値を、0 から 行数ひく1 まで順に変えて実行される。
これを実行すると、つぎのような warning (警告) メッセージが 5回出る。 わたしはまだその原因や対策を理解していない。 いまのところ、warning が出る以外は、まずいことはないように見える。
A value is trying to be set on a copy of a slice from a DataFrame
この欠損値の処理も、毎回データ読みこみにつづけてやるようにプログラムに組みこむことにする。
Python の matplotlib の 作図機能をつかって、データをグラフにすることができる。 matplotlib の plot は、2つの変数列を指定して散布図をつくるのが基本的な機能である。 しかし、1つの変数列だけを指定すれば、行を順にかぞえた番号を横軸にとったグラフをつくることができる。
たとえば、つぎのようにすれば、日照時間の時系列グラフをつくることができる。 (気温や降水量についても変数名がちがうだけで同様である。)
# 日時を明示しない時系列のグラフ (例、日照時間) plt.plot(df.sundur, 'o', markeredgecolor='k', markerfacecolor='none', markersize=4) plt.grid() #plt.show() imagefilename = 'plot_sundur.png' plt.savefig(imagefilename) plt.clf() plt.close()
このプログラムは、図を png 画像ファイルに書き出すようにしてある。 コメントにしてある「plt.show()」を生かし、その下の4行をコメントにすれば、図は画面に出る。
図はつぎのようなものになる。
ここで、横軸はデータの行を順にかぞえた番号で、 daily_tatebayashi.csv では、データのはじめの日 (2018年12月1日) を 0 (ゼロ) とした とおしの日付にあたる。 【 R による同様なプログラムでは、データのはじめの日を 1 とした とおしの日付をつかった。 それとは図の横座標値がちがう。】 日時をわかりやすく表現しなおすこともやりたいのだが、 わたしがまだ Python での日付のあつかいに熟練していないので、 ひとまずこのままにしておく。
ひとつひとつの観測値が、まる(〇)であらわされている。 記号(ここでは「まる」)の大きさを変えるには、プログラム中の markersize を変える。
散布図 (英語では scatterplot) の機能をためしてみよう。 つぎのようにして、日照時間と最高気温の散布図をつくることができる。
# 散布図の作図 (例、日照時間と最高気温)
plt.plot(df.sundur, df.tempmax, 'o', markeredgecolor='k', markerfacecolor='none', markersize=4)
plt.title('館林, 2018年12月～2020年11月', fontname='MS Gothic')
plt.xlabel('日照時間 (時)', fontname='MS Gothic')
plt.ylabel('最高気温 (℃)', fontname='MS Gothic')
plt.grid()
#plt.show()
imagefilename = 'scat_sundur_tempmax.png'
plt.savefig(imagefilename)
plt.clf()
plt.close()
【日本語 Windows 上で実行したばあい、matplotlib の title などの関数は 日本語文字を表示することができるが、それには日本語フォントを指定してやる必要がある。 それで、「 fontname='MS Gothic'」と入れてある。 他の OS 上 での実行では、この部分は変える必要があるだろう。】
3つ以上のデータ列について、総あたりの散布図をまとめてつくる機能がある。 縦横とも列の数だけ正方形の枠がならび、 対角線のところ以外の枠に散布図がかかれた形になる。 このようなグラフ表現を、英語では scatterplot matrix という。 日本語表現としては「散布図方陣」がよいと思う。
3列以上のデータフレーム df があれば、単に plot(df) とすれば、散布図方陣がつくられる。 しかし、ここでの df のうち、年月日、品質情報、均質番号の列は、総あたりからはずしたほうがよい。 そこで、subset機能をつかって df から観測値のはいった列だけをとりだして、 Pandas にふくまれた散布図方陣をつくる関数 plotting.scatter_matrix にわたすことにする。
# 散布図方陣 (日平均気温、最高気温、最低気温、日照時間、日降水量 の 総あたり) dfsub = df[['temp', 'tempmax', 'tempmin', 'sundur', 'precip']] pd.plotting.scatter_matrix(dfsub) #plt.show() imagefilename = 'scatterplotmatrix.png' plt.savefig(imagefilename) plt.clf() plt.close()
結果の図はつぎのようになる。(見出しや軸の目盛りなどの調整はまだしていない。) 対角線上の枠にはその変数のヒストグラムがかかれている。
データの頻度分布(度数分布)をみるのには、ヒストグラム(histogram)が基本的な方法である。 観測値のヒストグラムは、観測値を階級わけして集計し、 階級ごとの頻度を長方形の面積であらわす柱をならべた形のグラフにしたものである。
Pandas には ヒストグラムをつくる関数「hist」がある。ひとまずこれをすなおに実行して 日照時間のヒストグラムをつくってみる。
# ヒストグラム (例、日照時間)
df['sundur'].plot.hist()
plt.xlabel('日照時間 (時)', fontname='MS Gothic')
plt.ylabel('頻度', fontname='MS Gothic')
#plt.show()
imagefilename = 'hist_sundur.png'
plt.savefig(imagefilename)
plt.clf()
plt.close()
hist にはさまざまなオプションがあって、データ値の階級のわけかたや、 グラフの縦軸・横軸の尺度を変えることができるが、わたしはまだためしていない。
【matplotlib にも hist があって、Pandas の hist は それをよびだしている。 いまのばあい、hist をよびだしている行を「plt.hist(df.sundur)」に変えて matplotlib の hist を直接よびだしても、結果はかわらなかった。 ただし、オプションのあたえかたはちがう。】
【なお、「頻度」にあたる英語は「frequency」だが、ここでの意味は「振動数」や「周波数」ではない。】
【hist 関数をつかわずに、集計と作図を別々の手つづきによって、ヒストグラムをつくることもできる。 しかし、その手順の説明をまだ用意していない。】
ヒストグラムの情報を要約した表現として 「箱ひげ図」(box and whisker plot)がある。 たとえば、つぎのようにして、日平均気温、最高気温、最低気温の箱ひげ図をならべてつくることができる。
# 箱ひげ図 (例、日平均気温、最高気温、最低気温)
df.boxplot(column=['temp', 'tempmax', 'tempmin'])
plt.ylabel('℃', fontname='MS Gothic')
#plt.show()
imagefilename = 'boxplot_temp.png'
plt.savefig(imagefilename)
plt.clf()
plt.close()
箱の上下の端はデータの四分位値である。データ値の半分は箱の中にふくまれることになる。 中央値のところに横線がひかれている。 「ひげ」の上下の端はかならずしもデータの最大値・最小値ではなく、 「はずれ値」を「ひげ」からはずして個別に表示することがある。 ひとことで書けないので、Pandas の ドキュメントの boxplot のところを参照してほしい。
データを月 (month の列の値) によって分類し、それぞれの月の箱ひげ図をならべて表示することは、 boxplot 関数のオプション by によって簡単にできる。
# 箱ひげ図 (月ごとの日照時間)
df.boxplot(column='sundur', by='month')
plt.xlabel('月', fontname='MS Gothic')
plt.ylabel('日照時間 (時)', fontname='MS Gothic')
#plt.show()
imagefilename = 'boxplot_sundur_bymonth.png'
plt.savefig(imagefilename)
plt.clf()
plt.close()
1年を季節にわけることは、それ自体、研究課題になりうる問題でもある。 しかし、多数の地点の観測値をいっしょにあつかうばあいは、 各地点ごとに季節区分を考えるわけにいかず、 同じ季節区分をいっせいにあてはめたい。
大気上端にはいってくる太陽放射は、中緯度では、夏至のころに最大、冬至のころに最小になる。 地上の日射量や日照時間も、もし雲の影響が一定ならば、夏至のころに最大、冬至のころに最小になる。 したがって、夏至・冬至などを基準とした季節のわけかたも、有意義だ。
しかし、気象観測値は、グレゴリオ暦の月ごとに整理する慣習ができている。 そこで、気象学で世界の気候・天候をあつかうときは、 1年を、12・1・2月、3・4・5月、6・7・8月、9・10・11月にわけて「季節」(英語では season)とみなすことが多い。 北半球の中高緯度の感覚で、それぞれ、冬、春、夏、秋とよぶこともある。 しかし、季節名は南半球では逆になるし、熱帯では不適切なこともあるので、 英語の月の名の頭文字をとって、DJF, MAM, JJA, SON とすることも多い。
ここでは、CSVファイルをつくる際に日付を年、月、日にわけておいたから、 その「月」のところを読みとって、年は問わずに、「季節」ごとにまとめてヒストグラムをつくってみよう。
Pandas では、subset機能をつかって、データフレームの部分集合をつくることができる。 「df[条件]」で、条件にあった行だけが抽出され、あらたなデータフレームとしてつかうことができる。 (条件のうちの論理演算子が、Python ではふつう「and」「or」なのだが、 Pandas では「&」「|」をつかう。「&&」「||」ではない。) これによって、日照時間のデータを「季節」ごとにわけることができる。
ひとつのグラフィックス画面に複数のグラフをならべることは、 matplotlib の subplot 機能をつかうとできる。
四季の日照時間のヒストグラムを2段×2列にならべるプログラムをつくった。 比較のためには軸の尺度をそろえたほうがよいのだが、 現状のプログラムは尺度をそろえずそれぞれの季節ごとのデータに応じてきめている。
# 季節(3か月)ごとに抽出
dfdjf = df[(df.month == 12) | (df.month == 1) | (df.month == 2)]
dfmam = df[(df.month >= 3) & (df.month <= 5)]
dfjja = df[(df.month >= 6) & (df.month <= 8)]
dfson = df[(df.month >= 9) & (df.month <= 11)]
nrow = 2
ncol = 2
# 部分図の間隔をあける。wspace が左右の間隔、hspace が上下の間隔。
plt.subplots_adjust(wspace=0.4, hspace=0.6)
plt.subplot(nrow, ncol, 1)
dfdjf['sundur'].plot.hist()
plt.title('12・1・2月', fontname='MS Gothic')
plt.xlabel('日照時間 (時)', fontname='MS Gothic')
plt.ylabel('頻度', fontname='MS Gothic')
plt.subplot(nrow, ncol, 2)
dfmam['sundur'].plot.hist()
plt.title('3・4・5月', fontname='MS Gothic')
plt.xlabel('日照時間 (時)', fontname='MS Gothic')
plt.ylabel('頻度', fontname='MS Gothic')
plt.subplot(nrow, ncol, 3)
dfjja['sundur'].plot.hist()
plt.title('6・7・8月', fontname='MS Gothic')
plt.xlabel('日照時間 (時)', fontname='MS Gothic')
plt.ylabel('頻度', fontname='MS Gothic')
plt.subplot(nrow, ncol, 4)
dfson['sundur'].plot.hist()
plt.title('9・10・11月', fontname='MS Gothic')
plt.xlabel('日照時間 (時)', fontname='MS Gothic')
plt.ylabel('頻度', fontname='MS Gothic')
#plt.show()
imagefilename = 'hist_sundur_byseason.png'
plt.savefig(imagefilename)
plt.clf()
plt.close()
気象要素の種類によって、ヒストグラムの形がちがう。 気温は左右対称でわりあい正規分布に近いが、 降水量は値の大きいほうに長いすそをひく。 地上観測の雲量 (全天のうち雲でおおわれた部分のわりあい)は、0と1の両端の頻度が高い。 (ただし、平均期間を長くしていくと、正規分布にちかづく。) 日照時間も、雲量ほど極端ではないが、2つの山をもつ分布をしている。
藤部(2014)の本の図6.7に日ごとのいくつかの気象要素のヒストグラムがある。 【このうち東京の日降水量のヒストグラムで「40 mm」のすぐ右に 背の高い棒があるのは、40mm以上の値をひとつの柱にまとめてしまったための見かけのもので、 ここに頻度分布の山があるわけではない。】 図6.8には、月ごとの値のヒストグラムがある。 同じ図は、藤部(2015)の講義資料の5ページにもある。 (この講義資料は気象につかわれる統計学の専門的内容もふくんでいるが、 わかるところだけ読んでも有用だと思う。)
離れた地点について、同じ気象要素の値を作図して比較してみよう。
世界の気候の分布にとってはこまかい特徴だが、 日本の冬には、冬の季節風による雲・降水が多い地方と、少ない地方がある。 館林は冬の雲・降水が少ない関東地方の特徴をしめす地点である。 ここからわりあい小さな距離でありながら冬の雲・降水が多いところとして、 新潟県の地点をどこかえらんで、冬の降水量や日照時間を比較してみよう。
計算結果を図にするだけでなく、数値としてファイルに書き出すことも必要になるが、 わたしはまだ経験不足でうまく説明できない。わかったら補足したい。
データフレームを CSV ファイルに書き出すならば、to_csv がある。
(2地点のデータを、期間をそろえていっしょにあつかい、 同時の観測値について、散布図をつくったり、 相関係数を計算したりすることを、追加しようと思っていますが、 このページとは別にします。)