X Window (Windows パソコンでは Xming) を起動したうえで、 端末 (Windows パソコン では TeraTerm) で niigata にログインして、 各自のホームディレクトリで (この授業用のディレクトリをつくっている人はそこに cd したうえで) つぎのことをやってください。
「.csv」でおわる名まえのファイルが 1つと、 「.py」でおわる名まえのファイルが 10個 (変更するかもしれない) あるはず。
このうち、daily_kumagaya_202201_202405.csv は、 [気象庁の観測データを取得する] で、 気象庁ウェブサイトの「過去の気象データ ダウンロード」の機能をつかって、 ダウンロードした、熊谷の 2022年1月から2024年5月までの日ごとの値である。 (この下からもリンクしておく。)
Python 言語のスクリプトを、テキストエディタ (emacs) で作成する。 ファイル名は「.py」でおわるようにする。たとえば、「ex.py」としよう。
$ emacs ex.py
ex.py ができあがったら、端末画面にもどって、つぎのようにして実行する。
$ python ex.py
Pythonの変数名は ASCII コードであらわせるアルファベット・数字・ 下線文字 (「_」、アンダースコア) で構成する。はじめの字はアルファベットにする。 大文字と小文字は区別されるが、人がまちがえやすいので、 わたしの教材では、大文字にきまっているもの以外は、小文字だけをつかうことにする。 日本語文字なども Unicode であらわされるかぎりつかえるらしいが、ここではつかわないことにする。
これからつくる Python 言語のプログラムの中では、あらかじめ、つぎのように、 パッケージをとりこむことを予約し、略号を指定しておく。
# 50_import.py import pandas as pd import numpy as np import matplotlib.pyplot as plt import matplotlib.dates as mdates #####
この内容 (独立して動く Python プログラムではなく、断片である)を、 50_import.py として ex5 にふくめている。 emacs では、「c-x i ファイル名」で、編集中のバッファーに 外にあるファイルの内容をとりこむことができる。 これからつくる Python プログラムを編集するなかで、 「c-x i 50_import.py」で読みこんでほしい。
(略号を mdates としたパッケージが必要になるのは、日付を軸にとったグラフを かくときだけなので、それ以外のばあいは省略してよい。)
Python のパッケージ Pandas には「データフレーム」というデータ構造がある。 (R のデータフレームとは、よくにている。) おおざっぱには、縦横の表と思えばよい。 ただし、行はおたがいに同質のデータだが、列はそれぞれにちがった意味をもつことがある。
データフレームには名まえが必要で、名まえは自由につけられるが、 ここでは、だいたい1つだけのデータフレーム (および、それから派生したもの) をつかうので、 データフレーム名を(データフレームであることをわすれないように)「df」としてみる。
Pandas には、CSVファイルをデータフレームに読みこむ「read_csv」という関数がある。 (関数名が下線文字(「_」, underscore)をふくんでいることに注意。) これをつかって、CSVファイル 「daily_kumagaya_202201_202405.csv」の内容をデータフレーム「df」に読みこめる。
ダウンロードしたファイルには、頭から7行に説明があり、8行めから実質のデータがはじまる。 各列の項目の説明もあるが、そのままデータフレームの列名にするのには適していないので、 7行を全部読みとばし (skiprows=7)、ファイルの頭の部分から変数名を読みとらない (header=None)、という選択をした。 また、ファイルの文字コードが Shift_JIS であることを指定した。
列名が番号だとおぼえきれないので、実際の気象要素に対応する名まえにつけなおすことにする。 日本語ローマ字でもかまわないのだが、ここでは英語の略語にした。
各列の内容は、 たとえば df['tmp'] (データフレームのあとに、列名が、single quotation marks でかこまれ、さらに角かっこでかこまれている) のような形で取り出すことができる。 また、たとえば「df.tmp」(データフレーム名、「.」(ピリオド)、列名が、空白をはさまずにつなげられている) という形でも (Python文法上のほかの構文とまぎれなければ) 参照できる。
今回の例のデータでは、 列番号 (Python のつごうで、はじまりは 0)、わたしのPythonプログラムでの列名、気象データとしての項目、 表計算ソフトウェアで読んだときのワークシートの列記号は、 つぎのように対応している。
| 列番号 | 列名 | 項目 | ワークシートの列 |
|---|---|---|---|
| 0 | date | 年月日 | A |
| 1 | tmp | 平均気温 | B |
| 2 | q_tmp | 平均気温の品質情報 | C |
| 3 | h_tmp | 平均気温の均質番号 | D |
| 4 | tmpmax | 最高気温 | E |
| 5 | q_tmpmax | 最高気温の品質情報 | F |
| 6 | h_tmpmax | 最高気温の均質番号 | G |
| 7 | tmpmin | 最低気温 | H |
| 8 | q_tmpmin | 最低気温の品質情報 | I |
| 9 | h_tmpmin | 最低気温の均質番号 | J |
| 10 | precip | 降水量 | K |
| 11 | x_precip | 降水量の現象なし情報 | L |
| 12 | q_precip | 降水量の品質情報 | M |
| 13 | h_precip | 降水量の均質番号 | N |
| 14 | sundur | 日照時間 | O |
| 15 | x_sundur | 日照時間の現象なし情報 | P |
| 16 | q_sundur | 日照時間の品質情報 | Q |
| 17 | h_sundur | 日照時間の均質番号 | R |
date の列は、文字型で読みこまれているので、 日付型 (ここでは pandas パッケージの Timestamp 型) に変換しておく。 【日付型としては、datetime パッケージの datetime 型がつかわれることもある。】
ここまでの作業をするスクリプトはつぎのとおり。 ただし、1つの文を複数の行にわたって書くために、各行のおわりに「\」を置いている。 (「\」は、本来は「バックスラッシュ」 (右下がりの斜線) だが、 日本語環境では 円 (yen) マーク に見えることがある。)
# 51_read_csv.py
# データを csv ファイルから読みこみ、各列に名前をつける
df = pd.read_csv('daily_kumagaya_202201_202405.csv', \
header=None, skiprows=7, encoding='shift_jis', \
names=('date',\
'tmp', 'q_tmp', 'h_tmp',\
'tmpmax', 'q_tmpmax', 'h_tmpmax',\
'tmpmin', 'q_tmpmin', 'h_tmpmin',\
'precip', 'x_precip', 'q_precip', 'h_precip',\
'sundur', 'x_sundur', 'q_sundur', 'h_sundur') )
# 日付の列を Timestamp 型にする
df['date'] = pd.to_datetime(df['date'])
#####
この内容を、51_read_csv.py として ex5 にふくめている。 作図用の Python プログラムをつくる際に、読みこんでほしい。
今回のここからさきの例文では、いつも、50_import.py のような パッケージの予約と、 51_read_csv.py のようなデータ読みこみがあって、そのさきの作業がつづく。
51_read_csv.py は「気象官署等」用である。 そうでない一般のアメダス観測点用には、 「names=」で指定した一連のものから、「'x_precip',」と「'x_sundur',」をはずす。 そこさえ修正すれば、あとは変えなくてよい。
品質情報の値を見ると、今回のデータ例ではすべて「8」であるようだが、 ときどき他の数字がまざることがある。
「4」や「5」は、おそらく故障や機器点検があって、 24時間の観測値がそろっていないばあいで、 気象庁では、有効な件数をかぞえて、 なんとか利用可能と判断したばあい「5」、不適当と判断したばあいは 4以下の値にしている。
ここでは、品質情報が「8」でない観測値は欠損にしておこう。
Pandas のデータフレームの機能をつかうと、ある列の内容が特定の値であればそれを nan に変えることや、 nan であれば特定の値に変えることは、簡単にできるのだが、 「ある列の内容が特定の値であるとき、別のある列のその行の内容を nan に変える」手続きを、 わたしはまだ発見できていない。 そこで、手順だて型のプログラム言語に共通なループの機能をつかって、 つぎのようなプログラムを書いた。 ex5 に、52_miss_val.py としてふくめておく。 (ファイル名は missing values [欠損値] のつもり。)
# 52_miss_val.py
# 品質情報「8」以外の観測値を欠損値とする
for kt in range(len(df)):
if (df.q_tmp[kt] != 8):
df.tmp[kt] = np.nan
if (df.q_tmpmax[kt] != 8):
df.tmpmax[kt] = np.nan
if (df.q_tmpmin[kt] != 8):
df.tmpmin[kt] = np.nan
if (df.q_precip[kt] != 8):
df.precip[kt] = np.nan
if (df.q_sundur[kt] != 8):
df.sundur[kt] = np.nan
#####
「!=」は「ひとしくない」という、2つの数値から論理値を得る演算子、 「=」は代入演算 (左の変数に右の計算結果を代入する)である。 (「ひとしい」は「==」とかく。) len(df) で、dfの内容の行数(各列のデータ個数)が得られる。 今回の例では、対象期間の日数である。 この for ループは、変数 kt の値を、0 から 行数ひく1 まで順に変えて実行される。
これを実行すると、つぎのような warning (警告) メッセージが 5回出る。 ここでおこなったような代入文では、df の本体ではなく、 そのコピーが書きかえられるということだ。 しかしここからあとでつかわれるのはコピーのほうなので、実質的にまずいことはないと思う。 警告が出ないようにするにはどのようにプログラミングしたらよいか、わたし (増田) はまだ知らない。
A value is trying to be set on a copy of a slice from a DataFrame
この欠損値の処理も、 毎回データ読みこみにつづけてやるようにプログラムに組みこむべきである。 ただし、警告が出るのはうまくないし、今回のデータ例に関するかぎり 欠損の処理は必要ないので、今回の実習では取りこまなくてよいことにしよう。
Python の matplotlib の 作図機能をつかって、データをグラフにすることができる。 matplotlib の plot は、2つの変数列を指定して散布図をつくるのが基本的な機能である。
たとえば、つぎの 53_graph.py のように、日付 (date) を横軸、日照時間 (sundur) を縦軸にとれば、日照時間の時系列グラフをつくることができる。
(ここにしめしたもののまえに 51_import.py と 52_read_csv.py の内容がはいる。 このあとの例ではこの注意は省略する。)
# 53_graph.py
# 時系列のグラフ (例、日照時間)
plt.plot(df.date, df.sundur, 'o', markeredgecolor='k', markerfacecolor='none', markersize=4)
plt.grid()
formatter = mdates.DateFormatter("%y-%m")
plt.gca().xaxis.set_major_formatter(formatter)
plt.title('Daily sunshine duration [h] at Kumagaya')
plt.show()
##imagefilename = 'plot_sundur.png'
##plt.savefig(imagefilename)
##plt.clf()
##plt.close()
#####
このプログラムは、図を画面に表示するようにしてある。 「plt.show()」の行をコメントにし、 その下の4行を、「##」をはずして生かせば、図はpngファイルに出る。 pngファイルの名まえは「imagefilename =」の右辺で指定する。 図はつぎのようなものになる。
気象要素として日照時間を選択することは、 plt.plot の 2番めの位置引き数に plt.sundur をあたえることで指定している。 この sundur は df の列名である。 これを tmp, tmpmax, tmpmin, precip にかえることによって他の気象要素を作図できる。
気象要素を変更したばあいは、plt.title で指定される図の表題の文字列も変更しよう。 また、図を画像ファイルに書きだすばあいは、imagefilename で指定する 画像ファイル名も列名にあわせて変更しよう。
| 気象要素 | 列名 | 英語表題の例 |
|---|---|---|
| 日平均気温 | tmp | Daily mean temperature [degC] |
| 日最高気温 | tmpmax | Daily max. temperature [degC] |
| 日最低気温 | tmpmin | Daily min. temperature [degC] |
| 日降水量 | precip | Daily precipiation [mm/day] |
| 日照時間 | sundur | Daily sunshine duration [h] |
ひとつひとつの観測値が、まる(〇)であらわされている。 これは、plt.plot のオプションで指定している。 まず、3つめの位置引き数を「'o'」とすることで、 作図の形式は記号列であること、記号は まる であることを指定している。 つづいて、markeredgecolor で、記号の周の線の色を指定している。 「'k'」は「'black'」の省略形である。 それから、markerfacecolor で、記号の内部の色を指定している。 ここでは「None」で、記号の内部には色をぬらないことを指定している。 つぎに、markersize で、記号の大きさを指定している。 この数値を大きくすれば記号は大きく、小さくすれば記号は小さくなる。
formatter という文字列をふくむ2行は、 横軸の目盛りにつける文字列として日時をどのように書くかを調節するものである。 あらかじめ「import matplotlib.dates as mdates」を実行してあることを前提としている。
気象要素を変えて作図してみよう。
marker の種類、色、大きさなどを変えてみよう。
散布図 (英語では scatterplot) の機能をためしてみよう。 つぎのようにして、日照時間と最高気温の散布図をつくることができる。
# 54_scatter.py
# 散布図の作図 (例、日照時間と最高気温)
plt.plot(df.sundur, df.tmpmax, 'o', markeredgecolor='k', markerfacecolor='none', markersize=4)
plt.title('Kumagaya, Jan. 2022 - May 2024')
plt.xlabel('Sunshine duration [h]')
plt.ylabel('Daily max. temperature [degC]')
plt.grid()
plt.show()
##imagefilename = 'scat_sundur_tmpmax.png'
##plt.savefig(imagefilename)
##plt.clf()
##plt.close()
#####
3つ以上のデータ列について、総あたりの散布図をまとめてつくる機能がある。 縦横とも列の数だけ正方形の枠がならび、 対角線のところ以外の枠に散布図がかかれた形になる。 このようなグラフ表現を、英語では scatterplot matrix という。 日本語表現としては「散布図方陣」がよいと思う。
3列以上のデータフレーム df があれば、単に plot(df) とすれば、散布図方陣がつくられる。 しかし、ここでの df のうち、年月日、品質情報、均質番号の列は、総あたりからはずしたほうがよい。 そこで、subset機能をつかって df から観測値のはいった列だけをとりだして、 Pandas にふくまれた散布図方陣をつくる関数 plotting.scatter_matrix にわたすことにする。
# 55_scat_mat.py
# 散布図方陣 (日平均気温、最高気温、最低気温、日照時間、日降水量 の 総あたり)
dfsub = df[['tmp', 'tmpmax', 'tmpmin', 'sundur', 'precip']]
pd.plotting.scatter_matrix(dfsub)
# 図の縦横比を 1:1 にする
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
##imagefilename = 'scatterplotmatrix.png'
##plt.savefig(imagefilename)
##plt.clf()
##plt.close()
#####
結果の図はつぎのようになる。 対角線上の枠にはその変数のヒストグラムがかかれている。
(見出しや軸の目盛りなどの調整はまだしていない。)
データの頻度分布(度数分布)をみるのには、ヒストグラム(histogram)が基本的な方法である。 観測値のヒストグラムは、観測値を階級わけして集計し、 階級ごとの頻度を長方形の面積であらわす柱をならべた形のグラフにしたものである。
Pandas には ヒストグラムをつくる関数「hist」がある。ひとまずこれをすなおに実行して 日照時間のヒストグラムをつくってみる。
【縦軸の「frequency」の ここでの意味は「頻度」であり、 「振動数」や「周波数」ではない。】
# 56_histogram.py
# ヒストグラム (例、日照時間)
df['sundur'].plot.hist()
plt.xlabel('Sunshine duration [h]')
plt.ylabel('Frequency')
plt.show()
##imagefilename = 'hist_sundur.png'
##plt.savefig(imagefilename)
##plt.clf()
##plt.close()
#####
hist にはさまざまなオプションがあって、データ値の階級のわけかたや、 グラフの縦軸・横軸の尺度を変えることができる。 1時間ごとの降水量などについてためしてみた例が [1時間ごとの観測値をみる (1) ヒストグラム] にある。
【matplotlib にも hist があって、Pandas の hist は それをよびだしている。 いまのばあい、hist をよびだしている行を「plt.hist(df.sundur)」に変えて matplotlib の hist を直接よびだしても、結果はかわらなかった。 ただし、オプションのあたえかたはちがう。】
ヒストグラムの情報を要約した表現として 「箱ひげ図」(box and whisker plot)がある。 たとえば、つぎのようにして、日平均気温、最高気温、最低気温の箱ひげ図をならべてつくることができる。
# 57_boxplot.py
# 箱ひげ図 (例、日平均気温、最高気温、最低気温)
df.boxplot(column=['tmp', 'tmpmax', 'tmpmin'])
plt.ylabel('degC')
plt.show()
##imagefilename = 'boxplot_tmp.png'
##plt.savefig(imagefilename)
##plt.clf()
##plt.close()
#####
箱の上下の端はデータの四分位値である。データ値の半分は箱の中にふくまれることになる。 中央値のところに横線がひかれている。 「ひげ」の上下の端はかならずしもデータの最大値・最小値ではなく、 「はずれ値」を「ひげ」からはずして個別に表示することがある。 ひとことで書けないので、Pandas の ドキュメントの boxplot のところを参照してほしい。
データフレームの日付 (date) の列から月の情報をとりだし、month の列をつくっておく。 (この部分には Pandas の データフレームについての技巧をつかっている。) そして、boxplot 関数のオプション by に month を指定すれば、 それぞれの月ごとに箱ひげ図をつくり、ならべて表示することは、簡単にできる。 (ここでは、年はとわずに月でまとめることにした。)
# 58_boxplot_month.py
# 日付の列 date から 月の列 month をつくる
df['month'] = np.nan
for kt in range(len(df)):
df.loc[kt, 'month'] = df.date[kt].month
# 箱ひげ図 (月ごとの日最高気温)
df.boxplot(column='tmpmax', by='month')
plt.xlabel('month')
plt.ylabel('Daily max. temperature [degC]')
plt.show()
##imagefilename = 'boxplot_tmpmax_bymonth.png'
##plt.savefig(imagefilename)
##plt.clf()
##plt.close()
#####
【横軸は月の番号なので、「.0」がついているのはうまくない。 対策はあるはずだが、わたしはまだみつけていない。】
1年を季節にわけることは、それ自体、研究課題になりうる問題でもある。 しかし、多数の地点の観測値をいっしょにあつかうばあいは、 各地点ごとに季節区分を考えるわけにいかず、 同じ季節区分をいっせいにあてはめたい。
大気上端にはいってくる太陽放射は、中緯度では、夏至のころに最大、冬至のころに最小になる。 地上の日射量や日照時間も、もし雲の影響が一定ならば、夏至のころに最大、冬至のころに最小になる。 したがって、夏至・冬至などを基準とした季節のわけかたも、有意義だ。
しかし、気象観測値は、グレゴリオ暦の月ごとに整理する慣習ができている。 そこで、気象学で世界の気候・天候をあつかうときは、 1年を、12・1・2月、3・4・5月、6・7・8月、9・10・11月にわけて「季節」(英語では season)とみなすことが多い。 北半球の中高緯度の感覚で、それぞれ、冬、春、夏、秋とよぶこともある。 しかし、季節名は南半球では逆になるし、熱帯では不適切なこともあるので、 英語の月の名の頭文字をとって、DJF, MAM, JJA, SON とすることも多い。
ここでは、上の「月ごと」のばあいと同じ処理で日付の列から月の列をつくっておき、 年は問わずに、月を 3か月の「季節」ごとにまとめて、ヒストグラムをつくってみよう。
Pandas では、subset機能をつかって、データフレームの部分集合をつくることができる。 「df[条件]」で、条件にあった行だけが抽出され、あらたなデータフレームとしてつかうことができる。 (条件のうちの論理演算子が、Python ではふつう「and」「or」なのだが、 Pandas では「&」「|」をつかう。) これによって、データフレーム df を「季節」ごとにわけたデータフレーム dfdjf, djmam, djjja, djson をつくってみる。
ひとつのグラフィックス画面に複数のグラフをならべることは、 matplotlib の subplot 機能をつかうとできる。 【subplot の番号は 1 からはじまる。 range(4) の値は [0, 1, 2, 3] となるので、1をたしてつかっている。】
四季の日照時間のヒストグラムを2段×2列にならべるプログラムをつくった。 比較のためには軸の尺度をそろえたほうがよいのだが、 現状のプログラムは尺度をそろえずそれぞれの季節ごとのデータに応じてきめている。
# 59_hist_season.py
# 日付の列 date から 月の列 month をつくる
df['month'] = np.nan
for kt in range(len(df)):
df.loc[kt, 'month'] = df.date[kt].month
# 季節(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)]
stitles = ['Dec-Jan-Feb', 'Mar-Apr-May', 'Jun-Jul-Aug', 'Sep-Oct-Nov']
# 画面を 4つ (2 x 2) に 分割して作図
nrow = 2
ncol = 2
# 部分図の間隔をあける。wspace が左右の間隔、hspace が上下の間隔。
plt.subplots_adjust(wspace=0.4, hspace=0.6)
for iseason in range(4):
plt.subplot(nrow, ncol, iseason+1)
if iseason == 0:
dfdjf['sundur'].plot.hist()
elif iseason == 1:
dfmam['sundur'].plot.hist()
elif iseason == 2:
dfjja['sundur'].plot.hist()
else:
dfson['sundur'].plot.hist()
plt.title(stitles[iseason])
plt.xlabel('Sunshine duration [h]')
plt.ylabel('Frequency')
plt.show()
##imagefilename = 'hist_sundur_byseason.png'
##plt.savefig(imagefilename)
##plt.clf()
##plt.close()
#####
気温 (日平均、最高、最低のいずれか) と 降水量について、 同様に、3か月の季節ごとにまとめたヒストグラムをつくって見てみよう。
59のように3か月ごとにわけたそれぞれについて、 54のように、2つの気象要素の散布図をつくってみよう。
気象要素の種類によって、ヒストグラムの形がちがう。 気温は左右対称でわりあい正規分布に近いが、 降水量は値の大きいほうに長いすそをひく。 日照時間は、つぎにのべる雲量ほど極端ではないが、2つの山をもつ分布をしている。
(今回のデータ例にはふくめていないが) 雲量の観測値のヒストグラムは (日本の方式でいえば) 0 と 10 の両端にピークをもつことが多い。正規分布とはまったくちがう。 (ただし、たとえば1か月の平均をとると、かなり正規分布に近くなる。)
衛星観測や気象モデル (数値天気予報モデル、気候モデルの大気部分) から得られる (広域の) 雲量は、地表面積のうちで雲におおわれた部分の 面積の比率をさすことが多い。 この雲量と、地上観測による各地点の雲量とは、おおまかには対応するが、くわしく対応するとはかぎらない。
離れた地点について、同じ気象要素の値を作図して比較してみよう。
世界の気候の分布にとってはこまかい特徴だが、 日本の冬には、冬の季節風による雲・降水が多い地方と、少ない地方がある。 熊谷は冬の雲・降水が少ない関東地方の特徴をしめす地点である。 ここからわりあい小さな距離でありながら冬の雲・降水が多いところとして、 新潟県の地点をどこかえらんで、冬の降水量や日照時間を比較してみよう。
計算結果を図にするだけでなく、数値としてファイルに書き出すことも必要になる。 データフレームを CSV ファイルに書き出すならば、 df.to_csv(csvファイル名) とすればよい。 ただし ここで df としたところにはデータフレーム型の変数名がくる。
2地点のデータを、期間をそろえていっしょにあつかい、 同時の観測値について、散布図をつくったり、 相関係数を計算したりすることもやってみたい。
![[日別・月別気象要素の度数分布]](../zu/fujibe_histogram.png){kind=link}