続・Python による気象解析入門

観測点ごとの気象データを図示する


今回は、観測点ごとの気象データを、地図上の観測点の位置に表示する。

材料としては、気象庁による観測データ (気象庁がウェブサイトで公開しているもの) をつかう。 この観測データは csv ファイルで、縦に日時、横に地点と気象要素がならんだ形で提供されている。 このような表形式のデータを Python であつかうには、Pandas というパッケージをつかうのがよい。 また、日時の情報をあつかうには、datetime というパッケージをつかう。


バグ対策速報 [2023-07-04 晩, 更新 2023-07-05 朝]


実行方法

今回のPythonプログラムの例文とそれがつかうデータは、luce の /home/masuda/mapstationobs に置いてある。

各自、つぎのコマンドによって、プログラム例をコピーしてほしい (最後の空白のあとのピリオドが重要)。

cp -r /home/masuda/mapstationobs .

そして、mapstationobs に cd して作業してほしい。


地点表 (気象庁の気象官署等)

ここでは、「気象官署等」(気象台、測候所、特別地域気象観測所、ただし空港の観測所をのぞく) をあつかう。 その他のアメダスの観測点はあつかわない。

基本的なことは、[気象観測点を地図上に示す] のページでのべたとおりである。

ただし、気象庁ウェブサイトからダウンロードした観測データのCSVファイルでは、 地点番号がふくまれておらず、地点名が日本語の漢字・かな (文字コード Shift_JIS) でふくまれているので、 地点表のほうにも、漢字表記の地点名を含む列を追加した。 また、前回の地点表には富士山がふくまれていなかったが、今回は富士山もふくめた。 (/home/masuda/mapstation/jma_stations.csv も 富士山をふくむものにさしかえた。)

その各列の説明はつぎのとおり。

なお、地点名に別名のある地点について、stname_roma と stname_kanji との対応づけをつぎのようにしている。


気象庁の地上気象観測データ

どのように取得したか

気象庁ウェブサイトの「各種データ・資料」のページの「気象」のうちの「気象観測データ」のところに、 「過去の気象データ検索」、「過去の地点気象データ・ダウンロード」という二つのリンクがある。 どちらでもだいたい同じ情報がえられるが、形式がちがう。 csv 形式でほしいときは「ダウンロード」のほうをつかう。

ダウンロードを 2023年 7月 2日に実行した。

「ダウンロード」のページでは、 「地点を選ぶ」「項目を選ぶ」「期間を選ぶ」「表示オプションを選ぶ」をしてから 「CSVファイルをダウンロード」を実行すると、data.csv というファイルをダウンロードすることができる。 ただし、一回にダウンロードできるデータ量に制限があり、選択されたデータ量がそれをこすとダウンロードに進めない。

「地点を選ぶ」は、まず都道府県をえらび、都道府県の地図の中から地点をえらぶようになっている。 ただし、他の都道府県にうつっても指定はのこるので、複数の都道府県にまたがったデータをダウンロードすることができる。 CSVファイルに記録される地点の順序はここで指定したとおりになる。

ここでは、日本全国の「気象官署等」 (地点を選択する地図上で 赤の二重まる で表示される地点) を 選択しようとしたのだが、つぎにのべる項目と期間の選択と組みあわせると、 一度にダウンロードできるデータ量の制限をこえてしまう。 そこで、都道府県をいくつかずつまとめた地方ごとにダウンロードすることにした。 (地方のわけかたは下のファイルの説明のところでのべる。)

「項目を選ぶ」は、まずデータの時間間隔を指定する。ここでは「日別値」をえらんだ。 そして「気温」の部から「日平均気温」「日最高気温」「日最低気温」、 「降水」の部から「降水量の合計」、 「日照/日射」の部から「日照時間」をえらんだ。

【ちかごろ、アメダス観測点では、湿度の観測がはじまったのにともなって日照計がなくなり、 日照時間としては衛星観測からの推定値がふくまれていることがある。 今回選択したのは「気象官署等」なので日照時間の観測は継続されていると思うが、まだ確認していない。】

「期間を選ぶ」では、2023年 6月 1日から 2022年 6月 30日 までを指定した。 6月2日ごろの全国的な大雨と、記憶の新しい時期をふくめたかったからである。

「表示オプションを選ぶ」のところは、いずれも、あらかじめ与えられているほう (デフォールト値) にしておく。 とくに「ダウンロード CSV ファイルのデータ仕様」のうち「日付の形式」は「日付リテラルで格納」とした。 【増田の過去の教材では「年月日などに分けて格納」をつかったこともある。 そのときは datetime 型をつかわず、年、月、日をそれぞれ整数として処理した。】

取得されたデータファイル

ダウンロードされた data.csv のファイル名をつけなおし、つぎのところに置いた。 (いずれも、日本語文字コード Shift_JIS。)

csv ファイルのなかみをみると、6行の見出し情報のあと、7行め (Python のかぞえかたで第6行) からの各行がそれぞれ 1つの時点 (ここでは日別値なので 1日) のデータとなっている。 行数は、6 + 時点数なので、いまのばあい 36行、列数は、1 + 地点数×17 である。 表計算ソフトウェアに入れて目で見るには異常に横長のファイルだが、 Python で Pandas のデータフレームに読みこんでつかうにはさしつかえない。

各列の内容はつぎのようになっている。

「品質情報」などの意味は、気象庁の「ダウンロード」のページの右上の「CSVファイルの形式」というリンクのさきのページに書かれている。

品質情報の値は 8 が正常な観測値、5 以上が問題はあるが気象庁が使えると判断した値、 4 ~ 2 は数値はあるがあやしいもの、1 は欠損である。

【日平均値や日合計値は、1時間ごとの報告を24回あわせて計算している。 これが全部そろっていないが 80 % 以上あるばあいは品質情報 5、もっと欠けているばあいは品質情報 4 としている。】

「現象なし情報」は、降水量ならば「無降水」と「雨量計の計量能力未満の降水」を区別するためのものである。 今回選択した地点は いずれも「気象官署等」なので降水量も日照時間も「現象なし情報」の列がある。 【無人の観測所で現象なしの判定が実際にうまくできるか疑問はあるが。】 したがって、今回のデータは地点あたりの列の数が一定である。つぎにしめすプログラムではそのことを利用している。 【一般のアメダス観測点のばあいは「現象なし情報」の列がない。したがって、「気象官署等」と一般のアメダス観測点とがまざった CSV ファイルをつくると、 地点と気象要素から列番号を知るのが、ちょっとやっかいだ。それに対応する方法は別の機会に考えたい。】


[28c,d,e] 指定された緯度経度範囲にある観測点の ひとつの時点 (日) の ひとつの気象要素 (最高気温) を地図上にしめす

地図上の観測点の位置に「+」印をかき、気象要素の観測値を数字でしめすことにした。

プログラム例は、つぎの「28c」ではじまるものである。 「28a」, 「28b」はその作成途上のプログラムである。 【地方別に分割してダウンロードした観測値のファイルをあわせてつかえるように、また 観測値のファイルに出てくる地点の順序が地点表に出てくる順序と同じでなくてもよいように、 データフレームをいくつもつくってデータをつめかえている。その手順の概略はあとでのべる。】

【[2023-07-04 授業後補足] [28d] は [28c] の改訂版である。】 【[2023-07-07 補足] [28e] はさらにその改訂版である。】

「28c」のプログラムをそのまま実行すると、 北緯 34.5 ~ 38.5 度、東経 136.5 ~ 141度の緯度経度範囲にふくまれる観測点の位置に「+」印をかき、 その上側に、2023年 6月 28日の日最高気温を青色の数字で表示するようにしてある。

その事例には出てこないが、品質情報を見て、つぎのようなあつかいをしている。

気象要素を変える

「jcol = 4」であたえている数値を変更する。 これは観測データが読みこまれたデータフレーム dfobs の、 その気象要素が最初に出てくる列の番号で、csv ファイルの列番号 (上記) をひきついでいる。 【もっとわかりやすくするべきだが、ひとまずこうしてある。】 なお、図の表題の文字列や、図を画像ファイルに保存するときのファイル名も変更する。

気象要素jcolstitleの例ファイル名中の気象要素略号例
平均気温1Daily mean temperature (degC)tmp
最高気温4Daily maximum temperature (degC)tmpmax
最低気温7Daily minimum temperature (degC)tmpmin
降水量10Daily precipitation (mm/day)precip
日照時間14Sunshine duration (hr)sundur

【なお、品質情報の列番号は、「現象なし情報」がある場合は 観測値の列番号+2、 その他の場合は 観測値の列番号+1 である。 今回のプログラム例では、jcol の値が 7以下ならば、品質情報の列番号は観測値の列番号+1、 それよりも大きければ、観測値の列番号+2 としている。 これは、今回の教材用の観測データファイルにかぎって正しい条件であり、 気象庁からダウンロードした観測データ一般になりたつものではない。】

日付を変える

krecstart, krecend にあたえている値を変える。 krec は、読みこまれたデータにふくまれた最初の日時を 0 とした日時の通し番号である。 プログラム例では、krec の値を krecstart から krecend まで 1つずつ変えて、それぞれ作図するようにしてある。 今回のデータセットでは 6月 1日 を 0 としているので、27 は 6月 28日をさしている。 【これも、何月何日という形で指定できるようにするべきであり、 datetime または calendar パッケージをつかえばできるのだが、 いまの段階ではプログラムを単純にしたかったので、ひとまずこのようにした。】

図の表題の文字列や、図を画像ファイルに保存するときのファイル名の日付の部分は、 dfobs の 第0列から日付を読み、datetime パッケージの機能 strptime と strftime をつかって変換してあたえている。

対象地域を変える

経度と緯度の下限と上限を指定する変数 lonmin, lonmax, latmin, latmax にあたえる値を変える。

また、観測値のCSVファイルのうちどの地方のものを読みこむかを指定している [28c] では、「for region in ('kanto_chubu', 'tohoku'):」の文を書きかえる。 (その上にコメント行にしてある 「for region in ('hokkaido', 'tohoku', 'kanto_chubu', 'kansai', 'kyushu'):」を生かして 実際につかうものだけを残せばよい。) 【[2023-07-04 授業後補足] [28c] では、2か所の「for region in (...)」を同じように書きかえないといけない。 1か所かえればよいように修正したものを [28d] とした。 [28d] では「regionlist = [...]」という形の行を書きかえればよい。】


[27c] 地点の標高を図示する

気温にはあきらかに観測地点の標高の影響がある。 地点の標高を地図に表示するプログラムをつくった。前回の補足なので 27c という番号をつけたが、 今回用のバージョンを mapstationobs にふくめてある。


[29a ~ d] ひとつの地点の観測値の時系列をグラフにする

同じデータセットにより、ひとつの地点について、横軸に時間 (日付)、 縦軸に観測値をとったグラフをつくってみる。

降水量や日照時間は、0 がその量がないことをあらわすので、棒グラフが適切である。 気温は、0 ℃ は温度目盛りの中間点にすぎないので、折れ線グラフか、 記号 (marker) によるグラフが適切である。ここでは、折れ線と記号を併用した形にしている。

プログラム例の現在のバージョン (version 1 とよぶ) では、 気象庁からダウンロードしたファイルの観測値の列にある数値をそのままグラフにしており、 品質情報の列を見ていないので、データ欠損に対処できていない。 【地図表示の [28b, c] でやったように品質情報を見て表示内容を変えることはできるが、 (昨年につづいて今年も) プログラム例作成がまにあわなかった。それができたら、version 2 としたい。】

この節の例文では、matplotlib と Pandas と datetime をつかうが、Cartopy はつかわない。

この節の例文は、観測データのCSVファイルを 1つだけ読むようにつくってある。 [28c] のような複数のCSVファイルからのデータフレームを横につなぐことはしない。 (その需要はないと思うが、もしそうするならば、kstnsel をきめるところか、データフレームをつなぐところにも、変更が必要である。)

観測データをデータフレーム dfobs に読みこむ。 dfobs の (Python でいう) 第0列には、日付が「2023/6/1」のような文字列ではいっているが、 これを、datetime パッケージでつかわれる「datetime型」に変換しておく。 こうすると、日付の列を、matplotlib のグラフの軸の変数にとることができる。

地点の選択は、プログラムのはじめに、stnamesel という変数の値として、 ローマ字表記 (頭文字は大文字) の地名であたえることにした。 例文では Kumagaya があたえてあるが、他の地点名は、jma_stations_ja_sjis.csv を見てほしい。 (日本語文字がじゃまになるならば前回の jma_stations.csv でもよい。 ただし前回の表には Fujisan をふくめていなかった。その後 /home/masuda/mapstation/jma_station.csv も Fujisan をふくめるように修正した。)

matplotlib の subplot 機能をつかっているが、[29a, b, c] では、subplot は 1個だけである。 (これだけならば subplot をつかう必要はない。) [29d] では、画面を上下に 3分割して 3つの subplot をおさめている。 (このような変更をしやすくするために、はじめから subplot をつかう形にしたのだ。)

気温のグラフには、matplotlib の plot をつかっている。 plot では、横の座標値、縦の座標値につぐ 3つめの引数に、 記号 (marker) の種類をしめす「'o', '+', '^', 'v'」などの文字を指定すると、 記号列だけが表示される。 その指定をしないと、折れ線になる。 さらに「marker=」で記号の種類を指定すると、折れ線と記号列の両方が表示される。

降水量と日照時間の棒グラフには、matplotlib の bar をつかっている。

棒グラフで、値が図の枠の上限よりも大きいときは、棒が上限でうちきられていることに 注意してほしい。

横軸につける日付の目盛りの位置は自動的にきめさせているが、 現状では「06-29」と「07-01」の文字列がかさなって読みにくくなっている。 この対策はまだわからない。わかったらプログラム例を修正したい。

地点を変えるには、kstnsel にあたえる地点名を変えればよい。 filenameobs にあたえるCSVファイル名も変える必要があることもある。


データフレームをくみたてる手順の概略

観測値のCSVファイルが1個のばあい ([28a-b], [29a-d])

観測値のCSVファイルが複数個にわかれているばあい ([28c,d])、上のものにくわわる手順

[28e] での変更 [2023-07-07 追加]

これは、同じ地点群について複数のCSVファイルがあって地点の出現順序がちがうものを処理することに そなえた変更である。しかし実際に複数のデータをつなぐプログラムにはなっていない (手続きを「関数」にうまく分割する Python プログラミング技法をわたしがまだ修得していないので)。


2023-07-02, 2023-07-04, 2023-07-05, 2023-07-07
増田 耕一 (MASUDA Kooiti)
増田の [B, F, G 合同セミナーのページ] [立正大学用教材目次] [計算機を利用した気象分野の実習教材目次] [教材目次]