気象庁ウェブサイトの 「過去の気象データ・ダウンロード」から取得できる 熊谷の 2022年7月の全期間の 「時別値」の全気象要素を (2022年 8月 3日に) ダウンロードした。
この期間の熊谷では、雲量と日射量は、観測されておらず、品質情報が全部「1」になっている。
Pandas の read_csv で、df というデータフレームに読む。
Shift_JIS なので「encoding=shift_jis」がふつうだが、 Microsoft拡張がつかわれているので「encoding=cp932」にする。 このデータファイルには、 そのまま列名につかえる見出し行はなく、最初の6行は観測値ではないので、 「header=None, skiprows=6」とする。 また、そのまま index列 につかえる列もないので「index_col=None」とする。
各列に列名をつけるには、read_csv の中で names を指定してもよいのだが、 長くなるので、別に「df.columns」でやることにした。ここでの df はデータフレームにつけた名まえ。
各列が、文字型 (str)、整数型 (int)、浮動小数点型 (float、実数の近似なので「実数型」といわれることもある) であることを指定する。 「天気」(コード番号なので、ここでは weathercode とした) を整数型にしようとしたら、 この列にはデータ欠損 (csv では空、つまり文字数ゼロになっている) があり、 pandas はデータ欠損を NaN としようとするが整数型には NaN がないので、エラーになる。 weathercode は浮動小数点型にしておく。
「print(df)」をやると、行・列ともに途中が省略されてしまうが、 見えるかぎりではうまく読めているようだ。
日付は datex という列に読みこんだが、「日付リテラル」という文字列である。 datetime パッケージをつかって、datetime の日付型に変換する。 (datex の列の内容を書きかえてしまう。)
風向は winddir という列に読みこんだが、日本語 (Shift_JIS) の文字列である。 16方位なので、0～17の整数値に変換しておく (winddir16 という列にする。)
風向以外の気象要素について、品質情報が 8 でないばあいは データ欠損とみなし、気象要素の列の値を np.nan に変えておく。
気温は (日本の地点では) あきらかに日変化があるので、 日のなかの時刻が同じばあいにかぎって気温時別値の頻度分布をみることにする。 例として、正午 (12時) をあつかう。
datetime の機能をつかって、datex の列から「時」をとりだし、「hour」という列にしておく。
df のうちで hour が 12 である行をとりだし、別のデータフレーム「dfnoon」とする。
matplotlib の plt.hist の既定値で、dfnoon のヒストグラムをかいてみる。
matplotlib の fig と ax をわける方法をつかい、 ひとまず「subplotが1つだけある」形にして、 ax.hist の既定値で、dfnoon のヒストグラムをかいてみる。 結果は histtmp1 と同じ。
昼の気温は天気によってちがうことが予想されるので、 天気コード別に分類して (weathercode の値によって別々のデータフレームをつくり) ヒストグラムを subplot として上下にならべてみる。
天気コードの意味は、気象庁の「過去の気象データ・ダウンロード」のページの 右上の「CSVファイルの形式」のリンクさきのページの 「天気記号の対応表」のところに書いてある。
データを見ると、大部分の時刻の天気コードは 2, 4, 10, 15 のいずれかなので、 その 4つのばあいのヒストグラムをつくった。
それぞれ ax.hist の既定値の処理によったので、 図の横軸の範囲も、集計区間も、4つがばらばらである。
Pandas のデータフレームの hist メソッドをつかってみる。 これは中で matplotlib の hist をよびだしているにちがいなく、 結果は histtmp3 と実質的にちがわない。 ただし、格子線がひかれた。 (どちらの hist でも格子線をひくかどうかえらべるのだが、既定値がちがっているのだ。)
matplotlib の機能の ax.set_xlimit によって、 4つの subplot の横軸の範囲をそろえた。
ヒストグラムはそれぞれ既定値 (bins=10) で計算されているので、 集計区間はそろっていない。
hist に「bins=20, range=(20, 40)」と指定することによって、 20℃ から 40℃ まで、 1℃ きざみに統一した。 なお、横軸の目盛りは、既定値だと 2.5℃ 間隔になるのだが、 「22.5」などが表示されるとわかりにくいと判断して、5℃間隔を指定した。
降水量は、日本のばあいには日変化が大きくないので (熱帯では日変化が大きいことが多いのだが)、 件数をふやすことを優先して、全部の時刻をつかうことにする。
pandas の hist を既定値でつかってみた。 最小値と最大値のあいだを 10個の bin にわけているらしい。 降水量 約 4 mm 以下の区間の頻度が圧倒的に多いことしかわからない。
hist に 「bins=50, range=(0,50)」と指定し、 0 から 50 mm/hour まで 1 mm/hour きざみにしてみた。
これでも、0 をふくむ区間の頻度が圧倒的に多いので、ほかはよくわからない。
0.5 mm/hour 以上の降水があったときだけを対象とすることにする。 precip の値が 0.5 以上のばあいだけをとりだしたデータフレーム dfx をつくり、 dfx の precip のヒストグラムをつくる。 0 から 50 mm/hour まで、1 mm/hour きざみ。 この図は有用な情報をもっている。
降水量ゼロをはずしたので、降水量の対数をとることができる。 ここでは、常用対数 (10を底とする対数) をとって (NumPy の log10 によった)、 dfx に列 logprecip として追加し、 そのヒストグラムをすなおにつくってみた。
降水量自体よりはいくらか正規分布に近いとはいえるが、 ヒストグラムは極大値の左右に対称というわけではなく、正規分布とはだいぶちがう。
「bins=25, range=(-0.5, 2.0)」と指定して、区間きざみをこまかくしてみた。 (このきざみが妥当かどうかの検討はできていない。)
pandas の hist に オプションを指定しないで、風速のヒストグラムをかいてみた。
「bins=20, range=(0,10)」と指定して、0.5 m/s きざみにした。
(負になりえない量のヒストグラムではありがちなのだが) 左右非対称で、大きな値のほうにすそをひく分布になっている。
風速の平方根を計算し、df に 列 sqrtwindsp として追加したうえで、 そのヒストグラムをかいてみた。 左右対称に近く、正規分布に似た分布になった。
ここで平方根をとったのは、分布を左右対称にちかづける変数変換としてえらんだだけで、 物理的あるいは気象学的にとくに意味があるわけではない。
(つぎの項目は別のページに移動するかもしれない。)
pandas には describe() という、基本統計量を計算するメソッドがある。 これを、ひとまず、時間降水量について、まずそのまま、 それから 0.5 mm/時 以上の時にかぎって、適用してみる。 quantile としては標準の中央値と四分位値のほかに、0.9, 0.99, 0.999 を計算してみる。