この教材は、2024年度とほぼ同じですが、プログラム例の番号をつけかえています。
X Window (Windows パソコンでは Xming) を起動したうえで、 端末 (Windows パソコン では TeraTerm) で niigata にログインして、 各自のホームディレクトリで (この授業用のディレクトリをつくっている人はそこに cd したうえで) つぎのことをやってください。
「.csv」でおわる名まえのファイルが 1つと、 「.py」でおわる名まえのファイルが 11個 (変更するかもしれない) あるはず。
このうち、daily_kumagaya_202101_202505.csv は、 [気象庁の観測データを取得する] に書いた手順で、 気象庁ウェブサイトの「過去の気象データ ダウンロード」の機能をつかって、 2025年6月2日にダウンロードした、熊谷の 2021年1月から2025年5月までの日ごとの値である。 (この下からもリンクしておく。)
Pythonの変数名は ASCII コードであらわせるアルファベット・数字・ 下線文字 (「_」、アンダースコア) で構成する。はじめの字はアルファベットにする。 大文字と小文字は区別されるが、人がまちがえやすいので、 わたしの教材では、大文字にきまっているもの以外は、小文字だけをつかうことにする。 日本語文字なども Unicode であらわされるかぎりつかえるらしいが、ここではつかわないことにする。
Python 言語のスクリプトを、テキストエディタ (emacs) で作成する。 ファイル名は「.py」でおわるようにする。たとえば、「ex.py」としよう。
$ emacs ex.py
ex.py ができあがったら、端末画面にもどって、つぎのようにして実行する。
$ python ex.py
今回の Python プログラム例は、 準備にあたる 3つの部品 (単体として動作はしない) を、 「5a」「5b」「5c」ではじまり「.py」でおわるファイル名でつくってある。 実行プログラムは「51」から「58」ではじまり「.py」でおわるファイル名でつくってあるが、 それぞれのまえに「5a」「5b」「5c」の部品を置くことによって動作する形になる。
したがって、実行はつぎのようにするとよい。 「cat」コマンドを、本来の「catenate = 結合する」という目的につかう。 例として「51」をあげるが、「52」から「58」も同様である。
$ emacs 51*.py $ cat 5a*.py 5b*.py 5c*.py 51*.py > work.py $ python work.py
なお、このような構成のばあいの汎用的な実行手順は、 emacs で実行プログラムを編集するなかで c-x i によって部品をとりこんで、保存し、実行する、というものになる。
これからつくる Python 言語のプログラムの中では、あらかじめ、つぎのように、 パッケージをとりこむことを予約し、略号を指定しておく。
# 5a_import.py import pandas as pd import numpy as np import matplotlib.pyplot as plt import matplotlib.dates as mdates #####
(略号を mdates としたパッケージが必要になるのは、日付を軸にとったグラフを かくときだけなので、それ以外のばあいは省略してよい。)
Python のパッケージ Pandas には「データフレーム」というデータ構造がある。 (R のデータフレームや、リレーショナルデータベースの表に、よくにている。) おおざっぱには、縦横の表と思えばよい。 ただし、行はおたがいに同質のデータだが、列はそれぞれにちがった意味をもつことがある。
データフレームには名まえが必要で、名まえは自由につけられるが、 ここでは、だいたい1つのデータフレーム および それから派生したものをつかうので、 データフレーム名を(データフレームであることをわすれないように)「df」としてみる。
Pandas には、CSVファイルをデータフレームに読みこむ「read_csv」という関数がある。 (関数名が下線文字(「_」, underscore)をふくんでいることに注意。) これをつかって、CSVファイル 「daily_kumagaya_202101_202505.csv」の内容をデータフレーム「df」に読みこめる。
ダウンロードしたファイルには、頭から7行に説明があり、8行めから実質のデータがはじまる。 各列の項目の説明もあるが、そのままデータフレームの列名にするのには適していないので、 7行を全部読みとばし (skiprows=7)、ファイルの頭の部分から変数名を読みとらない (header=None)、という選択をした。 また、ファイルの文字コードが Shift_JIS であることを指定した。
列名が番号だとおぼえきれないので、実際の気象要素に対応する名まえにつけなおすことにする。 日本語ローマ字でもかまわないのだが、ここでは英語の略語にした。
各列の内容は、 たとえば df['tmp'] (データフレームのあとに、列名が、single quotation marks でかこまれ、さらに角かっこでかこまれている) のような形で取り出すことができる。 また、たとえば「df.tmp」(データフレーム名、「.」(ピリオド)、列名が、空白をはさまずにつなげられている) という形でも (Python文法上のほかの構文とまぎれなければ) 参照できる。
今回の例のデータでは、 列番号 (Python のつごうで、はじまりは 0)、わたしのPythonプログラムでの列名、気象データとしての項目、 表計算ソフトウェアで読んだときのワークシートの列記号は、 つぎのように対応している。
【[注意] 昨年度の例と列の順序がちがっている。 気象庁のサイトでダウンロードするとき「データの選択」で指定した順序によって、 csvファイルにはいる列の順序もちがってくる。 全部の気象要素を新たに指定すれば、平均気温、最高気温、最低気温とならぶはずなのだが、 今回は別の作業の指定が残っていてその他を追加したので この順序になったらしい。 このちがいに対応して、プログラムの read_csv をよびだすところの names のうちわけの順序を昨年度の例文から変えた。 】
| 列番号 | 列名 | 項目 | ワークシートの列 |
|---|---|---|---|
| 0 | date | 年月日 | A |
| 1 | tmp | 平均気温 | B |
| 2 | q_tmp | 平均気温の品質情報 | C |
| 3 | h_tmp | 平均気温の均質番号 | D |
| 4 | precip | 降水量 | E |
| 5 | x_precip | 降水量の現象なし情報 | F |
| 6 | q_precip | 降水量の品質情報 | G |
| 7 | h_precip | 降水量の均質番号 | H |
| 8 | tmpmax | 最高気温 | I |
| 9 | q_tmpmax | 最高気温の品質情報 | J |
| 10 | h_tmpmax | 最高気温の均質番号 | K |
| 11 | tmpmin | 最低気温 | L |
| 12 | q_tmpmin | 最低気温の品質情報 | M |
| 13 | h_tmpmin | 最低気温の均質番号 | N |
| 14 | sundur | 日照時間 | O |
| 15 | x_sundur | 日照時間の現象なし情報 | P |
| 16 | q_sundur | 日照時間の品質情報 | Q |
| 17 | h_sundur | 日照時間の均質番号 | R |
date の列は、文字型で読みこまれているので、 日付型 (ここでは pandas パッケージの Timestamp 型) に変換しておく。 【日付型としては、datetime パッケージの datetime 型がつかわれることもある。】
ここまでの作業をするスクリプトはつぎのとおり。
ただし、1つの文を複数の行にわたって書くために、 行のおわりだが文はさらに継続されるところに「\」を置いている。 (「\」は、本来は「バックスラッシュ」 (右下がりの斜線) だが、 日本語環境では 円 (yen) マーク に見えることがある。)
# 5b_read_csv.py
# データを csv ファイルから読みこみ、各列に名前をつける
df = pd.read_csv('daily_kumagaya_202101_202505.csv', \
header=None, skiprows=7, encoding='shift_jis', \
names=('date',\
'tmp', 'q_tmp', 'h_tmp',\
'precip', 'x_precip', 'q_precip', 'h_precip',\
'tmpmax', 'q_tmpmax', 'h_tmpmax',\
'tmpmin', 'q_tmpmin', 'h_tmpmin',\
'sundur', 'x_sundur', 'q_sundur', 'h_sundur') )
# 日付の列を Timestamp 型にする
df['date'] = pd.to_datetime(df['date'])
#####
5b_read_csv.py は、気象台・特別地域観測所などの「気象官署等」用である。 そうでない一般のアメダス観測点用には、 「names=」で指定した一連のものから、「'x_precip',」と「'x_sundur',」をはずす。 そこさえ修正すれば、あとは変えなくてよい。
品質情報の値を見ると、今回のデータ例ではすべて「8」であるようだが、 ときどき他の数字がまざることがある。
「0」と「1」は、観測ができておらず、有効な値がないところ (「欠測」) である。
「4」や「5」は、おそらく故障や機器点検があって、 24時間の観測値がそろっていないばあいで、 気象庁では、有効な件数をかぞえて、 なんとか利用可能と判断したばあい「5」、不適当と判断したばあいは 4以下の値にしている。
今回は、品質情報が「8」でない観測値は、統計からはずしたいので、欠損あつかいにする。 そのためのプログラムを書いた。 (ファイル名は missing values [欠損値] のつもり。)
# 5c_miss_val.py
# 品質情報「8」以外の観測値を欠損値とする
for kt in range(len(df)):
if (df.q_tmp[kt] != 8):
df.loc[kt, 'tmp' ] = np.nan
if (df.q_tmpmax[kt] != 8):
df.loc[kt, 'tmpmax'] = np.nan
if (df.q_tmpmin[kt] != 8):
df.loc[kt, 'tmpmin'] = np.nan
if (df.q_precip[kt] != 8):
df.loc[kt, 'precip'] = np.nan
if (df.q_sundur[kt] != 8):
df.loc[kt, 'sundur'] = np.nan
#####
Pandas のデータフレームの機能をつかうと、ある列の内容が特定の値であればそれを nan に変えることや、 nan であれば特定の値に変えることは、一斉にできるのだが、 「ある列の内容が特定の値であるとき、別のある列のその行の内容を nan に変える」手続きを、 わたしはまだ発見できていない。 そこで、手順だて型のプログラム言語に共通なループの機能をつかって、 行ごとに品質情報の列の値を見て観測値の列の値を変更するように、プログラムを書いた。
「!=」は「ひとしくない」という、2つの数値から論理値を得る演算子、 「=」は代入演算 (左の変数に右の計算結果を代入する)である。 (「ひとしい」は「==」とかく。)len(df) で、dfの内容の行数(各列のデータ個数)が得られる。 今回の例では、対象期間の日数である。
for ループは、変数 kt の値を、0 から 行数ひく1 まで順に変えて実行される。
【値の書きかえに、昨年度のプログラム例では、 「df.tmp[kt] = np.nan」のような構文をつかったのだが、 これだと「データフレーム自体を書きかえずそのコピーを書きかえた」というような 警告が出るので (それでも結果にはさしつかえないようだが)、 今年度は「df.loc[kt,'tmp'] = np.nan」のような構文に変更した。】
条件を、もし「8でないものは欠損」を「4以下のものは欠損」に変えるならば、 「!= 8」を「<= 4」とすればよい。
Python の matplotlib の 作図機能をつかって、データをグラフにすることができる。 matplotlib の plot は、2つの変数列を指定して散布図をつくるのが基本的な機能である。
たとえば、つぎの 51_timeplot.py のように、日付 (date) を横軸、日照時間 (sundur) を縦軸にとれば、日照時間の時系列グラフをつくることができる。 【ただし、今回の例文は、時系列グラフの手本としてはうまくないので、 graph ではなく timeplot とよぶことにした。】
(実行する際にプログラムには、このまえに 5a, 5b, 5c の内容がはいる。このあとの例ではこの注意は省略する。)
# 51_timeplot.py
# 2025-06-02 masudako
# 横軸に時間をとって観測値をプロット
#
# 図形出力をファイルに保存するか、画面表示か
savefig = False
#
# 横軸に時間をとって観測値をプロット (例、日照時間)
plt.plot(df.date, df.sundur, 'o', markeredgecolor='black', markerfacecolor='none', markersize=4)
plt.grid()
# - 横軸のラベルは、年2けた '-' 月２けたとする。
formatter = mdates.DateFormatter("%y-%m")
plt.gca().xaxis.set_major_formatter(formatter)
plt.title('Daily sunshine duration [h] at Kumagaya')
#
if savefig:
imagefilename = 'timeplot_sundur.png'
plt.savefig(imagefilename)
plt.clf()
plt.close()
else:
plt.show()
##### End of program #####
今回のプログラム例 [51]～[57] はいずれも、現状では図を画面に表示するようにしてあるが、 savefig に代入する値を True にすれば、図は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 で、記号の周の線の色を指定している。 ここでは「'black'」を指定している。 それから、markerfacecolor で、記号の内部の色を指定している。 ここでは「None」で、記号の内部には色をぬらないことを指定している。 つぎに、markersize で、記号の大きさを指定している。 この数値を大きくすれば記号は大きく、小さくすれば記号は小さくなる。
formatter という文字列をふくむ2行は、 横軸の目盛りにつける文字列として日時をどのように書くかを調節するものである。 あらかじめ「import matplotlib.dates as mdates」を実行してあることを前提としている。
気象要素を変えて作図してみよう。
[Optional] marker の種類、色、大きさなどを変えてみよう。
散布図 (英語では scatterplot) の機能をためしてみよう。 つぎのようにして、日照時間 (sundur) と 最高気温 (tmpmax) の散布図をつくることができる。
# 52_scatter.py
# 2025-06-02 masudako
# 散布図
#
# 図形出力を画像ファイルに保存するか、画面表示するか
savefig = False
#
# 散布図 (例、日照時間と最高気温)
plt.plot(df.sundur, df.tmpmax, 'o', markeredgecolor='black', markerfacecolor='none', markersize=4)
plt.title('Kumagaya, Jan. 2021 - May 2025')
plt.xlabel('Sunshine duration [h]')
plt.ylabel('Daily max. temperature [degC]')
plt.grid()
#
if savefig:
imagefilename = 'scat_sundur_tmpmax.png'
plt.savefig(imagefilename)
plt.clf()
plt.close()
else:
plt.show()
##### End of program #####
3つ以上のデータ列について、総あたりの散布図をまとめてつくる機能がある。 縦横とも列の数だけ正方形の枠がならび、 対角線のところ以外の枠に散布図がかかれた形になる。 このようなグラフ表現を、英語では scatterplot matrix という。 日本語では 【「散布図行列」という人が多いようだが、 「行列」が数学での意味でなく日本語でふつうの意味にとられると不適切なので、】 わたしは「散布図方陣」がよいと思う。
3列以上のデータフレーム df があれば、単に plot(df) とすれば、散布図方陣がつくられる。 しかし、ここでの df のうち、年月日、品質情報、均質番号、現象なし情報 の列は、総あたりからはずしたほうがよい。 そこで、subset機能をつかって df から観測値のはいった列だけをとりだして、 Pandas にふくまれた散布図方陣をつくる関数 plotting.scatter_matrix にわたすことにする。
# 53_scat_mat.py
# 2025-06-02 masudako
# 散布図方陣 (scatterplot matrix、総あたり散布図)
#
# 図形出力 画像ファイルに保存か、画面表示か
savefig = False
#
# 散布図方陣 (日平均気温、最高気温、最低気温、日照時間、日降水量 の 総あたり)
# - データフレーム df の利用する列だけをぬきだしたデータフレームをつくる
dfsub = df[['tmp', 'tmpmax', 'tmpmin', 'sundur', 'precip']]
# - 散布図方陣 (対角線のところはその列の変数のヒストグラム)
pd.plotting.scatter_matrix(dfsub)
#
if savefig:
imagefilename = 'scatterplotmatrix.png'
plt.savefig(imagefilename)
plt.clf()
plt.close()
else:
plt.show()
##### End of program #####
結果の図は、5×5の枠のそれぞれに散布図がかかれたものになる。 ただし、対角線上の枠にはその変数のヒストグラムがかかれている。
(見出しや軸の目盛りなどの調整はまだしていない。)
データの頻度分布 (度数分布) をみるのには、ヒストグラム (histogram) が基本的な方法である。 観測値のヒストグラムは、観測値を階級わけして集計し、 階級ごとの頻度を長方形の面積であらわす柱をならべた形のグラフにしたものである。
Pandas には ヒストグラムをつくる関数「hist」がある。ひとまずこれをすなおに実行して 日照時間のヒストグラムをつくってみる。
【「頻度」は英語で「frequency」というのだが、 この語は「振動数」・「周波数」という意味もあり、まぎらわしいので、 ここでは「出現回数」にあたる「Number of occurrences」とした。】
# 54_histogram.py
# 2025-06-02 masudako
# ヒストグラム (1変数の度数分布のグラフ)
#
savefig = False
#
# ヒストグラム (例、日照時間)
df['sundur'].plot.hist()
plt.xlabel('Sunshine duration [h]')
plt.ylabel('Number of occurrences')
plt.show()
#
if savefig:
imagefilename = 'hist_sundur.png'
plt.savefig(imagefilename)
plt.clf()
plt.close()
else:
plt.show()
##### End of program #####
hist にはさまざまなオプションがあって、データ値の階級のわけかたや、 グラフの縦軸・横軸の尺度を変えることができる。 1時間ごとの降水量などについてためしてみた例が [1時間ごとの観測値をみる (1) ヒストグラム] にある (このリンクさきページはちかぢか改訂予定)。
【matplotlib にも hist があって、Pandas の hist は それをよびだしている。 いまのばあい、hist をよびだしている行を「plt.hist(df.sundur)」に変えて matplotlib の hist を直接よびだしても、結果はかわらなかった。 ただし、オプションのあたえかたはちがう。】
ヒストグラムの情報を要約した表現として 「箱ひげ図」(box and whisker plot) がある。 たとえば、つぎのようにして、日平均気温、最高気温、最低気温の箱ひげ図をならべてつくることができる。
# 55_boxplot.py
# 2025-06-02 masudako
# 箱ひげ図 (1変数の度数分布の表示。上下の四分位値の間に箱)
#
savefig = False
#
# 箱ひげ図 (例、日平均気温、最高気温、最低気温 の3つの箱ひげ図をならべる)
df.boxplot(column=['tmp', 'tmpmax', 'tmpmin'])
plt.ylabel('degC')
if savefig:
imagefilename = 'boxplot_tmp.png'
plt.savefig(imagefilename)
plt.clf()
plt.close()
else:
plt.show()
##### End of program #####
箱の上下の端はデータの四分位値である。データ値の半分は箱の中にふくまれることになる。 中央値のところに横線がひかれている。 「ひげ」の上下の端はかならずしもデータの最大値・最小値ではなく、 「はずれ値」を「ひげ」からはずして個別に表示することがある。 ひとことで書けないので、Pandas の ドキュメントの boxplot のところを参照してほしい。
データフレームの日付 (date) の列から月の情報をとりだし、month の列をつくっておく。 (この部分には Pandas の データフレームについての技巧をつかっている。) そして、boxplot 関数のオプション by に month を指定すれば、 それぞれの月ごとに箱ひげ図をつくり、ならべて表示することは、簡単にできる。 (ここでは、年はとわずに月でまとめることにした。)
# 56_bpxplot_month.py
# 2025-06-02 masudako
# 箱ひげ図
# 日付のうちの月によって分類し (年は問わない)、それぞれの度数分布を見る。
#
savefig = False
#
# 日付の列 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]')
#
if savefig:
imagefilename = 'boxplot_tmpmax_bymonth.png'
plt.savefig(imagefilename)
plt.clf()
plt.close()
else:
plt.show()
##### End of program #####
【横軸は月の番号なので、「.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列にならべるプログラムをつくった。
横軸の数値の範囲は plt.xlim、縦軸の数値の範囲は plt.ylim で指定しているが、 これは、このデータセットで、しかも気象要素が日照時間のばあいについて 適切と判断した値であり、データセットや気象要素がちがうときは変更する必要があるだろう。
# 57_hist_season.py
# 2025-06-02 masudako
# 3か月ごとの四季にわけ (年は問わない)、それぞれヒストグラムをつくる。
#
savefig = False
#
# 日付の列 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)
# - 横軸の範囲 (日照時間用。他の気象要素では変更が必要。)
plt.xlim(-1, 15)
# - 縦軸の範囲 (このデータの日照時間用。他のデータではおそらく変更が必要。)
plt.ylim( 0, 200)
#
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('Number of occurrences')
#
if savefig:
imagefilename = 'hist_sundur_byseason.png'
plt.savefig(imagefilename)
plt.clf()
plt.close()
else:
plt.show()
##### End of program #####
[57] を変形し、 気温 (日平均、最高、最低のいずれか) と 降水量について、 同様に、3か月の季節ごとにまとめたヒストグラムをつくって見てみよう。
[57] のように3か月ごとにわけたそれぞれについて、 [52] のように、2つの気象要素の散布図をつくってみよう。
「plt.xlim(-1,15)」のかっこの中の数値を、 対象とする変数のありうる値をふくむように変更する。
たとえば「dfdjf['sundur'].plot.hist()」となっているところの 「( )」の中を、たとえば「(bins=20, range=(-20,20)」とする。 bins では区間の数、range では集計対象の数値の下限と上限を指定する。 上記の例では、-20 から 20 までを20 区間にわけるので、区間のきざみは 2 となる。
たとえば、降水量 1 mm 以上の日だけについて統計をとるとする。
「dfx = df[df.precip >= 1] 」を実行して、 あらたなデータフレーム変数 dfx をつくり、それ以後は df のかわりに dfx をつかえばよい。
[57] (四季) または [54] (全期間) のプログラムを変更した ヒストグラムの画像を、/home/reports/20250602 に提出。
savefig = True にすれば、図形は PNG ファイルに書き出される。
図の表題または x軸のラベルを、表示された気象要素にあわせること (英語でも日本語ローマ字でもよい。)
提出されるファイル名は学籍番号がはいっていればよいが、 気象要素名をふくめるならば、それも内容にあわせておくこと。
頻度分布の特徴を、グラフではなく、数値による要約で見たいことがある。 Pandas にはそのための describe という関数がある。 ここでは、[57] と同様な3か月ごとの「季節」について、 また、全期間について、describe をすなおに呼び出してみる。 結果を標準出力に print するようにしたので、 とくに指定しなければ端末画面に文字が出る。 (上記のように work.py をつくって実行したばあい、たとえば) 「python work.py > describe.txt」のように実行すれば、 describe.txt というテキストファイルに結果が保存される。
このプログラムは図形出力をしないので、 matplotlib は不要だが、便宜上 [5a] をとりこむことは変えていない。
# 58_describe_season.py
# 2025-06-02 masudako
# 3か月ごとの四季にわけ (年は問わない)、それぞれの記述統計の要約を書きだす。
#
savefig = False
#
# 日付の列 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)]
# 日付の列をのぞいたデータフレームをつくる
dfsub = df [['tmp', 'tmpmax', 'tmpmin', 'sundur', 'precip']]
dfdjfsub = dfdjf[['tmp', 'tmpmax', 'tmpmin', 'sundur', 'precip']]
dfmamsub = dfmam[['tmp', 'tmpmax', 'tmpmin', 'sundur', 'precip']]
dfjjasub = dfjja[['tmp', 'tmpmax', 'tmpmin', 'sundur', 'precip']]
dfsonsub = dfson[['tmp', 'tmpmax', 'tmpmin', 'sundur', 'precip']]
# 記述統計の要約を書きだす
print('DJF ----------')
print(dfdjfsub.describe())
print('MAM ----------')
print(dfmamsub.describe())
print('JJA ----------')
print(dfjjasub.describe())
print('SON ----------')
print(dfsonsub.describe())
print('All ----------')
print(dfsub.describe())
print('End ----------')
##### End of program #####
結果の「All」 (全期間) の部分はつぎのようになっている。
All ----------
tmp tmpmax tmpmin sundur precip
count 1610.000000 1603.000000 1604.000000 1607.000000 1607.000000
mean 16.121242 21.657205 11.490399 6.395706 3.157747
std 8.612120 8.720864 9.110324 4.139782 9.602766
min -1.300000 2.300000 -6.400000 0.000000 0.000000
25% 8.300000 14.100000 3.200000 2.300000 0.000000
50% 16.100000 21.500000 11.100000 7.500000 0.000000
75% 23.000000 28.400000 19.300000 9.700000 0.500000
max 33.600000 40.000000 29.600000 13.900000 129.500000
End ----------
countは件数、meanは平均、stdは標準偏差、minは最小値、maxは最大値。
25%、50%、75%は分位値 (パーセンタイル)で、 25%が下側四分位値、50%が中央値、75%が上側四分位値にあたる。
(この対象期間のとりかたは、1～5月が5年ぶん、6～12月が4年ぶんで、 季節の重みが不公平である。 この地点の年間の気象要素の頻度分布の特徴をみるためには、 月によってちがう重みづけをするべきだが、今回はその考察を見送った。)
なお、describe にオプションを指定すれば、 たとえば 90%、99% などの分位値を計算させることもできる。
気象要素の種類によって、ヒストグラムの形がちがう。 気温は左右対称でわりあい正規分布に近いが、 降水量は値の大きいほうに長いすそをひく。 日照時間は、つぎにのべる雲量ほど極端ではないが、2つの山をもつ分布をしている。
(今回のデータ例にはふくめていないが) 雲量の観測値のヒストグラムは (日本の方式でいえば) 0 と 10 の両端にピークをもつことが多い。正規分布とはまったくちがう。 (ただし、たとえば1か月の平均をとると、かなり正規分布に近くなる。)
衛星観測や気象モデル (数値天気予報モデル、気候モデルの大気部分) から得られる (広域の) 雲量は、地表面積のうちで雲におおわれた部分の 面積の比率をさすことが多い。 この雲量と、地上観測による各地点の雲量とは、おおまかには対応するが、くわしく対応するとはかぎらない。
離れた地点について、同じ気象要素の値を作図して比較してみよう。
世界の気候の分布にとってはこまかい特徴だが、 日本の冬には、冬の季節風による雲・降水が多い地方と、少ない地方がある。 熊谷は冬の雲・降水が少ない関東地方の特徴をしめす地点である。 ここからわりあい小さな距離でありながら冬の雲・降水が多いところとして、 新潟県の地点をどこかえらんで、冬の降水量や日照時間を比較してみよう。
計算結果を図にするだけでなく、数値としてファイルに書き出すことも必要になる。 データフレームを CSV ファイルに書き出すならば、 df.to_csv(csvファイル名) とすればよい。 ただし ここで df としたところにはデータフレーム型の変数名がくる。
2地点のデータを、期間をそろえていっしょにあつかい、 同時の観測値について、散布図をつくったり、 相関係数を計算したりするのもおもしろいだろう。
![[日別・月別気象要素の度数分布]](../zu/fujibe_histogram.png){kind=link}