今回は、観測点ごとの気象データを、地図上の観測点の位置に表示する。
材料としては、気象庁による観測データ (気象庁がウェブサイトで公開しているもの) をつかう。 この観測データは csv ファイルで、縦に日時、横に地点と気象要素がならんだ形で提供されている。 このような表形式のデータを Python であつかうには、Pandas というパッケージをつかうのがよい。 また、日時の情報をあつかうには、datetime というパッケージをつかう。
今回のPythonプログラムの例文とそれがつかうデータは、luce の /home/masuda/mapstation に置いてある。
各自、つぎのコマンドによって、プログラム例をコピーしてほしい (最後の空白のあとのピリオドが重要)。
cp -r /home/masuda/mapstation .
そして、mapstation に cd して作業してほしい。
気象庁による気象観測については、 気象庁ウェブサイトの「知識・解説」のページの「気象」のうちの 「気象に関する観測と予報の技術の解説」のうちに、「気象の観測」というリンクがあり、 それをたどると「気象観測・気象衛星」というページがある。
そのうち「地上気象観測所について」というリンクのさきの「地上気象観測」というページに 地上気象観測所の一覧表や分布地図がある。 そこからさらに「地上気象観測 地点情報履歴ファイル・統計切断情報ファイル」というリンクがある。 それをたどって「地上気象観測 地点情報履歴ファイル」(smaster.index.zip, タイムスタンプ 2022-02-25) をダウンロードし、 それを材料として、つぎのファイルを作成した。(作成手順の説明は、ひとまず省略する。)
その各列の説明はつぎのとおり。
気象庁の地上気象観測点には、このほかに、アメダス (AMeDAS) の観測点がある。 【気象台・測候所・特別地域気象観測所などではないアメダス観測点は、 WMO 地点番号をもっていない。】
ここで仮に「日本中央部領域」と呼ぶ領域は、北緯 32 ～ 42度、東経 132 ～ 142 度である。 (前回の「日本中央部領域」とは、同じ領域名をつかってしまったが、緯度の範囲をかえた。 今回は陸上をなるべくふくめたかったからである。)
地点表の csv ファイルをそれぞれ Pandas の「データフレーム」 dfstn1, dfstn2 に読みこんでいる。
作図は、まず fig のなかにひとつの subplot をつくり、そこに正距円筒図法の地図をつくり、 経線・緯線、海岸線をひく。ここまではすでにやったとおり。
さらに、都道府県境をひいている。 【[2022-08-26 改訂] 授業のときつかった例文プログラムでは、 Cartopy の機能にたよって Natural Earth というサイトが提供するデータを自動的に読みこむようにしたが、 これではユーザーごとの状況によって動作が不確実だとわかったので、 2022-08-03 から 2022-08-26 までに、 [日本の都道府県境をかく / 中国の「省」の境をかく]の ページに書いた方法をつかうように例文を変更した。 区別のため、改訂版プログラムのファイル名に「_v20220803」をふくめている。 授業のときつかった例文はそのままのファイル名で残してあるが、そちらをつかうことはおすすめしない。 luce 以外の計算機で改訂版プログラムをつかう人は、 [日本の都道府県境をかく / 中国の「省」の境をかく]の ページを見て県境のデータを取得してからプログラムを動かしてほしい。】
dftsn1, dfstn2 をそれぞれ 1行ずつ読んで、 緯度経度が図の領域内であれば、その位置に記号「+」をかき、その上に地点名 (ローマ字表記) をかくようにした。 dfstn1 の地点は赤で、dfstn2 の地点は青で表示することにした。
(text のオプションの ha は horizontal alignment (左右の位置合わせ)、ここでは center (中央) とする。 va は vertical alignment (上下の位置合わせ)、ここでは bottom (下側) とする。 「指定した点が文字列の下側にくるようにする」ので、文字列が指定した点の上にくることになる。)
ここで仮に「関東領域」と呼ぶ領域は、北緯 34.5 ～ 38.5 度、東経 137 ～ 141 度である。 【授業中に見せたページの記述はまちがっていました。授業後に修正しました。】 実際は「関東甲信越」にちかい。ただし新潟県の北のほうがすこしはずれ、静岡県はふくまれる。
つぎの作業のために、この図の範囲にふくまれた観測点の地点名 (ローマ字表記) のリストを作成した。
気象庁ウェブサイトの「各種データ・資料」のページの「気象」のうちの「気象観測データ」のところに、 「過去の気象データ検索」、「過去の地点気象データ・ダウンロード」という二つのリンクがある。 どちらでもだいたい同じ情報がえられるが、形式がちがう。 csv 形式でほしいときは「ダウンロード」のほうをつかう。
ダウンロードを 2022年 7月 16日に実行した。
「ダウンロード」のページでは、 「地点を選ぶ」「項目を選ぶ」「期間を選ぶ」「表示オプションを選ぶ」をしてから 「CSVファイルをダウンロード」を実行すると、data.csv というファイルをダウンロードすることができる。 ただし、一回にダウンロードできるデータ量に制限があり、選択されたデータ量がそれをこすとダウンロードに進めない。
「地点を選ぶ」は、まず都道府県をえらび、都道府県の地図の中から地点をえらぶようになっている。 ただし、他の都道府県にうつっても指定はのこるので、複数の都道府県にまたがったデータをダウンロードすることができる。 CSVファイルに記録される地点の順序はここで指定したとおりになるので、 多地点のデータを、複数回ダウンロードしてあわせてつかおうとすると、かなり注意が必要だ。 ここでは、上記の station_list_kanto.txt に書いた順に地点を指定した。 【実は、1地点ぬかしたままダウンロードまで進んでしまい、やりなおした。】
「項目を選ぶ」は、まずデータの時間間隔を指定する。ここでは「日別値」をえらんだ。 そして「気温」の部から「日平均気温」「日最高気温」「日最低気温」、 「降水」の部から「降水量の合計」、 「日照/日射」の部から「日照時間」をえらんだ。
【ちかごろ、アメダス観測点では、湿度の観測がはじまったのにともなって日照計がなくなり、 日照時間としては衛星観測からの推定値がふくまれていることがある。 今回選択したのは「気象官署等」なので日照時間の観測は継続されていると思うが、まだ確認していない。】
「期間を選ぶ」では、2022年 6月 1日から 2022年 7月 15日 までを指定した。 6月末から7月はじめの猛暑と、そのあとの大雨をふくめたかったからである。
「表示オプションを選ぶ」のところは、いずれも、あらかじめ与えられているほう (デフォールト値) にしておく。 とくに「ダウンロード CSV ファイルのデータ仕様」のうち「日付の形式」は (増田はふつう「年月日などに分けて格納」をつかっているが) ここでは「日付リテラルで格納」とした。
【7月2日から4日に、気象庁が契約していた au (KDDI) の通信網の故障で、データの欠損が生じている。 気象庁ウェブサイトのダウンロードのページには つぎのようなお知らせが出ている。
7月2日から4日のアメダスデータは、通信障害によりデータが欠落している場合があります。 復旧まで2週間程度かかる見込みですので、ご利用の際はご注意下さい。（2022.07.04）
今回ダウンロードしたデータを見ると、この事情による欠損があらわれている。 気象庁はリアルタイムでとどかなかった観測値を回収しようとしているのだが、 7月16日現在、まだ回収できていなかったのだ。】
【[2022-08-26] その後に回収がすすんだようなので、あらたにダウンロードしたものを別のファイル名で置く。】
ダウンロードされた data.csv のファイル名をつけなおし、つぎのところに置いた。
csv ファイルのなかみをみると、6行の見出し情報のあと、7行めからの各行がそれぞれ 1つの時点 (このばあいは 1日) のデータとなっている。
各列の内容はつぎのようになっている。
地点数が 40 なので、列の数は 1 + 17 × 40 で、681 である。 表計算ソフトウェアに読みこむと、列記号が A, ... Z, AA, AB, ... と進んで「ZE」にいたる。 このような横長の表を表計算であつかうのはむずかしいが、 Python の Pandas であつかうのはむずかしくない。
「品質情報」などの意味は、「ダウンロード」のページの右上の「CSVファイルの形式」というリンクのさきのページに書かれている。
品質情報の値は 8 が正常な観測値、5 以上が問題はあるが気象庁が使えると判断した値、 4 ～ 2 は数値はあるがあやしいもの、1 は欠損である。
【日平均値や日合計値は、1時間ごとの報告を24回あわせて計算している。 これが全部そろっていないが 80 % 以上あるばあいは品質情報 5、もっと欠けているばあいは品質情報 4 としている。】
「現象なし情報」は、降水量ならば「無降水」と「雨量計の計量能力未満の降水」を区別するためのものである。 今回選択した地点は いずれも「気象官署等」なので降水量も日照時間も「現象なし情報」の列がある。 【無人の観測所で現象なしの判定が実際にうまくできるか疑問はあるが。】 したがって、今回のデータは地点あたりの列の数が一定である。つぎにしめすプログラムではそのことを利用している。 【一般のアメダス観測点のばあいは「現象なし情報」の列がない。したがって、「気象官署等」と一般のアメダス観測点とがまざった CSV ファイルをつくると、 地点と気象要素から列番号を知るのが、ちょっとやっかいだ。】
[27b] と同じひろがりの地図上の観測点の位置に、気象要素の観測値を数字でしめすプログラムを用意した。
気象要素としては、まず最高気温をあつかう。 ほかの気象要素にかえるには「jcol = 4」であたえている数値を変更する。 これは観測データが読みこまれたデータフレーム dfobs の、 その気象要素が最初に出てくる列の番号で、csv ファイルの列番号 (上記) をひきついでいる。 【もっとわかりやすくするべきだが、ひとまずこうしてある。】 なお、図の表題の文字列や、図を画像ファイルに保存するときのファイル名も変更する。
【[2022-07-19 授業後補足] 授業中に luce に置いてあった [28a] と つぎの [28b] のプログラムは、 降水量と日照時間の図が正しくできませんでした。 品質情報の列番号がまちがっていました。 授業後の19時すぎに、luce と、このウェブページのプログラムを 改訂版にさしかえました。 改訂では、jcol の値が 7以下ならば、品質情報の列番号は観測値の列番号+1、 それよりも大きければ、観測値の列番号+2 としました。 これは、今回の教材用の観測データファイルにかぎって正しい条件であり、 気象庁からダウンロードした観測データに共通ではありません。】
どの日をとりあげるかは、データフレーム dfobs の行番号を「krec = 30」のように指定するようにしてある。 この行番号は、観測データのある最初の日 (今回のデータセットでは 6月 1日) を 0 としているので、30 は 7月 1日をさしている。 【これも、何月何日という形で指定できるようにするべきであり、 datetime または calendar パッケージをつかえばできるのだが、 いまの段階ではプログラムを単純にしたかったので、ひとまずこのようにした。】 図の表題の文字列や、図を画像ファイルに保存するときのファイル名の日付の部分は、 dfobs の 第0列から日付を読み、datetime パッケージの機能 strptime と strftime をつかって変換してあたえている。
プログラムの中では、まず jma_observatories.csv と jma_stations2.csv を読み、データフレームを連結して dfstn とする。 さらに、station_list_kanto.txt も読み、「関東領域」内の地点だけのデータフレーム df をつくる。
また、観測値の csv ファイルを データフレーム dfobs に読みこむ。
データフレーム df に、観測値の列 val と、品質番号の列 qval を追加する。 そこに入れる内容が、対象となる時点 (日) によってちがう。 val の列に、dfobs の、krec で指定された行、jcol から 17 ずつふえていく番号で指定された列の値を入れていく。 qval の列に、そのつぎの列の内容を入れていく。
【[2022-07-24 補足] プログラム例では、val と qval を、それぞれリストとして作成し、 その内容を df の列に追加してからつかったのだが、 計算の流れをみなおしたら、df の列にせずにリストのままつかえばよく、 そのほうがプログラムが簡単になることがわかった。 そのように改訂した [28a] と [28b] の改訂版を、このページに追加しておく。 (luce の mapstation の下にも置く。)】
ここから作図。だいたい [27b] と同じ。 ただし、地点の上側に表示する文字列が、地名ではなく、つぎのものである。
[28a] と同じ図を複数の日についてつくる。 データを読みこむのは一度でよいが、 作図作業は、krec を krecstart から krecend まで変えてループする。
同じデータセットにより、ひとつの地点について、横軸に時間 (日付)、 縦軸に観測値をとったグラフをつくってみる。
降水量や日照時間は、0 がその量がないことをあらわすので、棒グラフが適切である。 気温は、0 ℃ は温度目盛りの中間点にすぎないので、折れ線グラフか、 記号 (marker) によるグラフが適切である。ここでは、折れ線と記号を併用した形にしている。
プログラム例の現在のバージョン (version 1 とよぶ) では、 気象庁からダウンロードしたファイルの観測値の列にある数値をそのままグラフにしており、 品質情報の列を見ていないので、データ欠損に対処できていない。 【地図表示の [28a, b] でやったように品質情報を見て表示内容を変えることはできるが、 プログラム例作成が今回の授業日程に まにあわなかった。 それができたら、version 2 としたい。】
【[2022-07-19 授業後補足] 授業中に置いたプログラムには「Year 2002」という文字列がありましたが、 これは「Year 2022」の書きまちがいでした。 授業後の 19時ごろ、luce に置いたものも、このウェブページのリンクも、訂正版にさしかえました。】】
この節の例文では、matplotlib と Pandas と datetime をつかうが、Cartopy はつかわない。
まず、data_kanto_daily_20220601_20220715.csv をデータフレーム dfobs に読みこむ。 dfobs の (Python でいう) 第0列には、日付が「2022/6/1」のような文字列ではいっているが、 これを、datetime パッケージでつかわれる「datetime型」に変換しておく。 こうすると、日付の列を、matplotlib のグラフの軸の変数にとることができる。
地点名のローマ字表記を得るために、 station_list_kanto.txt を、もうひとつのデータフレーム dfstnsel に読みこむ。
地点の選択は、dfstnsel の行番号によることにした。 熊谷ならば 11 である。(日常用語でいえば 12行めだが、 Python では最初の行を 0 とするので 11 になる。 「kstn = 12-1」のように、1をひくことを計算機にやらせてもよい。)
matplotlib の subplot 機能をつかっているが、[29a, b, c] では、subplot は 1個だけである。 (これだけならば subplot をつかう必要はない。) [29d] では、画面を上下に 3分割して 3つの subplot をおさめている。 (このような変更をしやすくするために、はじめから subplot をつかう形にしたのだ。)
気温のグラフには、matplotlib の plot をつかっている。 plot では、横の座標値、縦の座標値につぐ 3つめの引数に、 記号 (marker) の種類をしめす「'o', '+', '^', 'v'」などの文字を指定すると、 記号列だけが表示される。 その指定をしないと、折れ線になる。 さらに「marker=」で記号の種類を指定すると、折れ線と記号列の両方が表示される。
降水量と日照時間の棒グラフには、matplotlib の bar をつかっている。
棒グラフで、値が図の枠の上限よりも大きいときは、棒が上限でうちきられていることに 注意してほしい。