技術面は [2022年12月の大雪に関連して地上気象観測値の空間分布を見る] の つづきなのだが、 事例を、2023年2月10日の関東地方の大雪にかえたので、別のページにする。
気象庁のアメダスの地点リストは (最新のものではないが) [気象観測データ (多地点) の空間分布と時系列の作図 (Python による試作)] の 教材ページで説明した、つぎのファイルをつかうことにする。 これは、気象庁ウェブサイトから 2020年12月19日にダウンロードしたファイルを増田が編集したものである。
アメダスのなるべく多くの観測点のデータをつかうためには、 気象庁ウェブサイトの 「各種データ・資料 > 過去の気象データ・ダウンロード」の 「地点を選ぶ」で都道府県 (北海道のばあいはそれを分割した地区) の地図を表示させたうえで その左上にある「..県全地点」を選択するのがよさそうだ。
今回は、関東地方全都県 (茨城県、栃木県、群馬県、埼玉県、千葉県、東京都、神奈川県)、 山梨県、長野県、静岡県、新潟県の全地点の、 日平均気温、日最高気温、日最低気温、 降水量の日合計 (ここでは「降水量」とよぶ)、日照時間、 日最深積雪 (「積雪深」)、降雪量の日合計 (「降雪深」)の、 2023年2月9日から2023年2月11日までの日別値をダウンロードした。
ほぼ、気象庁のダウンロードページのデフォルト (とくに指定しなかったばあいに選択されるオプション) にしたがったが、 「表示オプションを選ぶ」の「その他」のうちの「都道府県名を格納」にチェックを入れた。 (アメダス観測点名は全国で一意的についておらず、同じものがある。 これまでに、「府中」 (東京都と広島県)、「美山」 (福井県と京都府) があることがわかっている。 それを都道府県名で区別しようと考えたからである。 ただし、今回試作したプログラムでは、都道府県名を見ていない。)
ところが、作業をはじめてから、プログラムがエラーでとまってしまう原因をしらべるうちにわかったのだが、 ダウンロードされた csv ファイルには、同じ都道府県のうちに同じ地点名の列群が複数あって、 わたしがいま作成しているプログラムの理屈のままでは、ただしく処理できない。 ただし、その一方の列群には、実質的な観測値がはいっていない。 気象庁のダウンロードページは、長期にわたるデータをあつかうので、 最近観測をしている地点の値がほしいばあいにも むかし観測をしていた地点の値がほしいばあいにも対応しようとして、 このような形になっているのだ。 ダウンロードされた csv ファイルから、観測値がふくまれないほうの列群を削除すればよいのだが、 Excel で列を削除してふたたび csv で保存したものを同じプログラムで読んだら、期待した動作にならなかった。 おそらく Excel での保存の際の指定をかえれば対策できると思うが、それがわからない。
ここでは、気象庁のダウンロードのページの「地点を選ぶ」にもどって、 各都道府県の地図あるいは「選択された地点」を見て、 同じ地点名のうちで、観測がおこなわれていないほう (地点を選択すると見える説明に「観測終了」という記述がふくまれるほう) の指定をはずして、ダウンロードすることにした。
該当する地点名は、つぎのものである。 (気象庁が観測点の移転としてあつかったたばあいは、これにはふくまれていない。) 今回の対象とした都道府県についてはこれで全部だと思うが、 ほかの都道府県にもあるだろう。
そのようにしてダウンロードした data.csv の名まえをつけかえたものが、つぎの csv ファイルである。
(ファイル名の daily は、日別値であることを明示するために入れた。 gkanto は、関東を主としてもうすこし広い地域なので、"greater Kanto" のつもりで入れた。 with_pref は、都道府県名をふくめたことのおぼえのために入れた。 pref は「県」の英語訳の prefecture のつもりである。 cleaner は、同じ名まえの古い観測点をはずしたことのおぼえのために、 うまい表現を思いつかなかったのでくるしまぎれで入れた。 「-er」は比較級のつもりである。 「with_pref」や「cleaner」は、いつもそのようにしているばあいは不要だろう。)
このファイルのなかみには、観測点名はあるが、観測点の緯度経度も、地点番号もふくまれていない。 緯度経度を知るには、地点リストとつきあわせる必要がある。 観測値ファイルに観測点が出現する順序は、地点リストと同じではない。 なお、アメダスの地点名は、日本全国で一意的ではなく、同じ名まえのことがある。 都道府県内では一意的にきまるようだ。 【ただし、あとでみるように、ちがう時代の同じ地点名がちがう位置に対応していることはある。】 そこで、ひとまず、amedas_stn.csv から、対象とする都道府県の地点だけをぬきだしたファイルをつくり (いちばん上の見出し行は残した)、 それと照合することにした。
また「気象官署等」とされる地点はアメダス地点でもあるが、その地点のデータには「現象なし情報」の列がある。 それ以外のアメダス地点のデータには「現象なし情報」の列がない。 そこで、1地点あたりの列の数をかぞえておくことにする。 観測値ファイルと地点リストとの地点の照合をこころみたプログラムがつぎのものである。 このプログラムは、地点ごとに、観測値ファイルでのその地点の最初の列、列の数、 地点名、緯度、経度を、標準出力に書きだす。
このプログラムは、「都道府県名を格納」を指定してダウンロードしたデータにあわせたので、 pd.read_csv() のオプション で skiprows=3 とした。 都道府県名を格納しないばあいは、skiprows=2 とすればよい。
ダウンロードのサイトでつくられるcsvファイルに書かれる地点名と、 地点リストの地点名の文字列が同じでないところについて、あらかじめ対策が必要である。 これまでに気づいた次の2地点については、プログラムで対応した。
| amedas_stn.csv の 地点名 | 気象庁からダウンロードされた csvファイル の 地点名 | jma_stations.csv の 地点名 |
|---|---|---|
| つくば | つくば（館野） | Tateno |
| 奥日光 | 奥日光（日光） | Nikko |
標準出力をファイルに redirect することができるが、 Anaconda Powershell Prompt で実行すると、結果のファイルの文字コードが UTF-16 になる。 その文字コードをあつかえるエディタなどで読みこんで、UTF-8 または shift_jis に変えて 保存するとよい。UTF-8 にコード変換した結果を、つぎにつける。 Spyder で実行したばあいどうなるか、わたしはまだ確認していない。
今回取得したアメダス地点の気象要素のいずれかの日別値を、 まえの例と同じ仮称「関東」の範囲の地図に表示するプログラム例をつくった。 (今回の地点の選択に対しては西の端は東経 137.5 度のほうがよかったが、 まえと変えないことにした。)
データファイルに記録のある最初の日時を 0 として、 krecstart から krecend までの日時についてそれぞれ作図するようにしてある。
jvar にあたえる値によって、気象要素をえらぶようにしてある。
降水量のデータ欠損は一時的なものとおもわれるので「NA」と表示することにしたが、 その他の気象要素のデータ欠損は、ほとんどが観測機器が設置されていないところなので 地点の表示もしないことにした。
都道府県境をかいているが、その方法は、 [日本の都道府県境をかく / 中国の「省」の境をかく] のページからリンクした ne_10m_admin_1_japan.zip をダウンロードして展開してあることを前提としている。 map_amedas_daily_gkanto_202302case.py の変数 shapefilepath の値として与える文字列を、 自分がそのデータを展開したフォルダ名 (スラッシュでおわる) に変更してほしい。 【都道府県境が不要ならば、 プログラムの「都道府県境のデータを準備する」のところと「都道府県境を描く」のところを消すなり コメントにするなりすればよい。そのデータをダウンロードする必要もない。】
プログラムのなかでは、観測値がはいった 1つのファイルについて read_csv を 2回実行している。 1回めは地点名を読んで上の chkstn と同様な処理をするため、 2回めは観測値を読むためである。
このプログラムは、「都道府県名を格納」を指定してダウンロードしたデータにあわせた。 都道府県名を格納しないばあいは、2回の pd.read_csv() のオプションの skiprows の値をそれぞれ変更すればよい。
jvar 不要。 「# タイトルをかく」 の jvar で分岐しているところを 1つの選択肢だけのこす。
joffset 不要。「# 観測値のデータフレームを構築する」の 「if dfstn.ncols[kstn] ... 」 から 「joffset = 3*4 + ...」 まで不要。 「kcol = ...」の文は「kcol = dfstn.kcolb[kstn]」となる。
地点あたりの列の数は 3 または 4 なので、 「# 品質情報 (色わけにつかう)」のところをつぎのようにする。
if (dfstn.ncols[kstn] == 4):
qval[kstn] = int(dfobs.iloc[krec, kcol+2])
else: # (... == 3)
qval[kstn] = int(dfobs.iloc[krec, kcol+1])