[気象庁の観測データを使う (1)](板倉用教材)、 [気象観測データのグラフ作成と基本的な統計処理 (Python による試作)]に つづいて、複数の観測点の気象データをつかって現象の空間分布をみることにする。 【[2021-05-05] Cartopy をつかって 地図上にデータを表示した。ただし地図投影方法は正距円筒図法だけである。】
この教材のもとになった R による教材をつくりはじめた2020年12月17日に、北関東の大雪のニュースがあった。 「藤原」という地名を聞いた。 (東武線の新藤原のあたり(これは栃木県)かと思ったが、そうではなかった。) 群馬県みなかみ町の藤原だ。ここには気象庁のアメダスの観測所がある。 わたしはこの事例を専門的に追いかけてはいないが、 報道や気象庁ウェブサイトで知ったことを紹介しておきたい。
気象庁ウェブサイトの「防災情報」の中の「気象情報」のページには、 気象庁やその地方気象台などが発表した、注意すべき気象状況の情報についての文章が収録されている。 全国のものと、府県別のものがある。
群馬県のところを見ると、前橋地方気象台が出した「大雪に関する群馬県気象情報」がある。 2020年12月20日夜現在最新のものは、12月20日10時59分発表の「第6号」である。 この「号」の番号はときどき1からつけなおされている。 12月16日17時00分発表の(前回の)「第6号」に、つぎのような部分がある。 (箇条書きの部分の数値は、かっこの外が12時間降雪量、かっこの中が積雪深である。 なお、わたしはこの「降雪量」を「降雪深」とよぶことにしている。 別記事[雪について観測されている量](板倉用教材) も見てほしい。)
［雪の実況］ １６日１６時現在の１２時間降雪量（積雪の深さ）は みなかみ町藤原 ３２センチ（１６３センチ） みなかみ町幸知 ２２センチ（ ７６センチ） となっています（アメダス速報値）。また、みなかみ町藤原では１２月１５日の積雪差日合計が１１７センチとなり、１９８９年の統計開始以来の極値を更新しました。（これまでの極値は１９９１年１月６日の８７センチ）
なお、この「みなかみ町幸知」はアメダスの地点「みなかみ」と思われる。 ただし、前橋地方気象台のウェブサイト (https://www.jma-net.go.jp/maebashi/) のお知らせ記事によれば、 アメダス「みなかみ」は2020年11月26日に 幸知から湯原に移転しているので、たぶん正しくは「みなかみ町湯原」での観測値だろう。 (なお、下の教材中でふれる気象庁のアメダス地点表のファイル(12月19日にダウンロードしたもの)では、みなかみ観測点の所在地は湯原となっている。)
また、2020年12月18日の毎日新聞の報道 (有料記事だがhttps://mainichi.jp/articles/20201218/ddl/k10/040/079000c?pid=14509)にはつぎのようにあった。
前橋地方気象台によると、７２時間降雪量は、みなかみ町藤原で同日午後４時までに２１９センチに達し、統計開始以来最多を更新した。 ... 同町みなかみでも７２時間降雪量が１５９センチとなり観測史上１位を記録。同町藤原の積雪量は１７日午後５時現在で２メートル超となった。
なお、今回のこの地点の積雪深は新記録ではない。 藤原の積雪深の極値は 2006年1月26日の 301 cm、 みなかみは 2006年1月28日の 275 cm である (気象庁の「過去の気象データ検索」サイトで検索した)。
そこで、時空間分布を見る対象を、この大雪の事例をふくめるように、 2020年12月の、群馬県北部を中心とする約 150 km四方の領域にふくまれる地点の、 気象庁のアメダスの観測データとした。
観測値の空間分布を図示するためには、観測地点の位置を知る必要がある。 地球上の位置は、ふつう、緯度・経度 (および 高さ) で表現される。
前回の教材でとりあげた気象庁ウェブサイトの「過去の気象データ ダウンロード」では、 地点を選択する際に、その地点の緯度・経度がポップアップの形で表示される。 しかし、ダウンロードされたCSVファイルには緯度・経度の情報ははいっていない。 ポップアップを書きとるのはうまくない方法だ。 気象庁の観測地点の緯度・経度をふくむ一覧表をつかうべきだ。
気象庁ウェブサイトの 「知識・解説」の「気象」のうち「気象に関する観測と予報の技術の解説」のうちに「気象の観測」というリンクがある。 そのページの表題は「気象観測・気象衛星」となっている。 そのページの「地上・地域気象観測」のところに、 「地上気象観測所について」と「地域気象観測システム(アメダス)について」がある。 気象台などの(いわゆる)気象官署での観測は「地上気象観測所について」のほうに書かれているのだが、 ほとんどの気象官署はアメダスの観測点にもなっているので、 ここではアメダスのほうにしぼってみることにする。
「地域気象観測システム(アメダス)について」 http://www.jma.go.jp/jma/kishou/know/amedas/kaisetsu.html のページには、アメダスではどんな観測をしているかの解説とともに、 「観測所一覧」がある。 このうち、CSV形式の「地域気象観測所一覧 [ZIP圧縮形式]」のファイル http://www.jma.go.jp/jma/kishou/know/amedas/ame_master.zip をつかう。 2020年12月19日にダウンロードしたものを、つぎのところにコピーしておく。 ファイル形式は zip で、展開すると ame_master.csv というCSVファイル (日本語文字コードはShift_JIS)がはいっている。
「観測所一覧」のページには PDFファイルもあるが、その内容は、このファイルと、もうひとつの「地域気象観測所一覧（雪） 」との内容を そのままあわせたもののようだ。
このファイルに収録された地点は、2020年12月現在のアメダスの観測点である。 過去に観測していたが現在は観測していない地点のことや、観測所の移動経歴などを知りたいときは、 同じページからリンクされている「アメダス地点情報履歴ファイル」を見る必要がある。
アメダスの地点には5けたのID番号がついている。この番号は国際地点番号ではない。
【[注意] 2021年にはいって、アメダスで湿度の観測がはじまっている (地点ごとに開始時期がことなる) が、 ここにとりあげたファイルは、湿度の観測がはじまるまえの状況を反映したものである。】
気象庁の ame_master.csv は、そのまま表計算ソフトウェアで表として見ることもできるし、 Python の Pandas の read_csv で「データフレーム」に読みこむこともできる。
Pandas の read_csv では header=0 と指定すれば入力データの最初の行の文字列を列名にすることができるが、 ame_master.csv の最初の行に書かれた日本語の列名は変数名にむかないので、 header=None, skiprows=1 でよみとばし、names(...) で列名をあたえることにした。
【 read_csv による読みこみでは、こまごまとした困難があった。 あとの作業で出あったものをふくめて、おぼえがきとして列挙しておく。
しかし、わたしがやりたい作業にとってはちょっと不便な形をしているので、あらかじめ組みかえておくことにする。
組みかえたおもな点は、緯度・経度の単位である。 気象庁の地点表の緯度・経度は、それぞれ「度」と「分」の列にわかれている。 「度」のほうには整数値がはいり、「分」のほうには小数点下1けたまでの値がはいっている。 緯度・経度の表現は、60進法の度・分・秒 (秒よりこまかいところが必要ならば秒の10進小数)か、 度とその10進小数か、がふつうだと思うが、気象庁の地点表は度・分と分の10進小数なのだ。 わたしは、「度」の単位であつかいたいので、「度」の列の値に「分」の列の値の60分の1をたして換算する。 【なお、日本国内の地点ならばすべて北緯・東経なので簡単なのだが、 世界の地点をあつかうときは南緯・西経のあつかいをきめておく必要がある。 南緯・西経を負にするならば、度の列だけでなく分の列の値も符号をかえる対象になる。】組みかえのもうひとつは、地名のカタカナ表記である。 観測所名の表記は漢字が基本となっているが、それだけでは読みかたがわからない。 気象庁の地点表にはカタカナ表記がある。 ただし、1バイトのカタカナの文字コード (漢字の半分の幅で表示されることが多かったので「半角カナ」とよばれることが多い) が つかわれている。 その体系では、濁点・半濁点はそれがつくカナ文字とは独立した文字のようにあつかわれている。 いまでは1バイトカタカナには対応していないソフトウェアが多いので、 漢字をふくむ日本語文字コードのカタカナ (いわゆる「全角カナ」) に変換しておきたい。 ウェブ検索でみつけた horio さんのブログに書かれた R プログラム をつかった。 Python では、horio さんのプログラムをもとに Python に書きかえたものをくみこんだ。
積雪深の観測地点は別の表になっているが、ame_master.csv の「備考1」の列に 積雪深の観測点としての地点IDがはいっている。 例外的に積雪深だけの観測点があるかもしれないが、積雪深を観測するアメダスの地点のほとんどはこれでわかる。 積雪深を観測していない地点は「－」(ダッシュ)がはいっているが、NA にかえておくことにする。
ame_master.csv には「所在地」の列がある。住所にあたる地名 (番地まではない) が書かれている。 それにくわえて、気象官署でもある地点のばあいは、官署名がはいっている。 (くぎり文字として「全角」の空白がはいっている。) 官署名だけ拾い出しておくことにする。
ame_master.csv には IDも地名も同じ地点が複数行ふくまれていることがある。 観測機器がちょっと離れたところにわかれて置かれているのだ。 緯度・経度・高さがちがうといってもわずかなので、さきに出てきたほうだけを残すことにする。
以上の方針で組みかえをした。
観測点の地理的分布を見るには、地図上に表現したほうがよい。 Python では、Cartopy というパッケージをつかえばできる。
しかし、ここでは、ひとまず Cartopy をつかわず、 一般的な 2次元空間に分布するデータの表示の例をかねて、 横軸に経度、縦軸に緯度をとった x-y 平面上の位置として、観測地点の位置を表示することにする。 これも一種の地図投影方法とみなすことができ、そのばあいは「正距円筒図法」とよばれる。 【地図投影方法については、2000年ごろ地理情報学・地図学の教材としてつくった [地図投影]のページを見てほしい。】実際の地球では、一定の経度間隔が、高緯度ほど短い距離間隔に対応する。 正距円筒図法の地図上の距離はそうなっていないことに注意する必要がある。 緯度と経度を同じ尺度で画面に表示すると、赤道に近いところ以外では、現実よりも横長になる。 つぎにしめす例は、経度幅 2度、緯度幅 1.5度をあつかっていて、 画面上の図の寸法は matplotlib の default 値にまかせてしまったので、 尺度のたてよこ比をよく考えていないのだが、やや横長だと思う。
アメダスの地点のうち 北緯 36～37.5度、東経 138.5～140.5度の範囲にあるものの位置と 地点名を、amedas_stn.csv のデータにもとづき、Pandas のデータフレームに読みこんで、 matplotlib の機能で表示した。
観測している気象要素を、記号の形と色の両方でしめした。
(R による例では図の下側に記号の凡例を入れたのだが、この Python による例ではそれがまだできていない。)
地理的背景として、海岸線、地形、県境などをいっしょにしめすことがのぞましいのだが、 残念ながら準備ができていない。 地図帳、または、国土地理院の「地理院地図」( https://maps.gsi.go.jp/ ) などを併用して見てほしい。
ここでつかった地点表ファイルは日本全国をふくんでいるので、
これとちがった地域についての図にかえることは、
プログラムの緯度経度範囲 (xlimit, ylimit の値として与えている) を変えればすぐできる。
(ただし、図の下の凡例をしめす位置も、あわせて変える必要がある。)
しかし、このあとしめす気象データは、上の図の緯度経度範囲だけダウンロードしたものなので、
他の地域についてやるには、ダウンロードからやりなおす必要がある。
Cartopy というパッケージをつかうと、緯度・経度であたえられた地点の情報を、 地図投影変換をして表示することができる。
ここでは、上と同じアメダスの地点分布図を、地図上に表示して、海岸線と県境をいれてみた。 いろいろな地図投影法をつかえることがのぞましいのだが、 いまのところ成功したのは正距円筒図法 (英語では「Plate Carree」) だけである。 (メルカトル図法をこころみたらエラーになった。)
海岸線と都道府県境のデータは、Cartopy が標準的につかっている Natural Earth というウェブサイトのものをつかった。 ここでは海岸線のほかに世界の多くの国のおもな行政区画の境界線データセットをつくっている。 これを読んで作図する方法は、つぎのウェブページを参考にした。 これは中国の大行政区画 (省、自治区) の境界を表示したものだが、 China を Japan に変更すれば日本の都道府県境につかえる。
https://geocat-examples.readthedocs.io/en/latest/gallery/Overlays/NCL_overlay_11b.html GeoCAT-examples NCL_overlay_11b.py (C) Copyright 2021, University Corporation for Atmospheric Research. Revision 4e24dbc1.
緯度経度の値はそれぞれの軸の両端にも入れようとしたのだが、なぜかうまくいかなかった。
空間的にも時間的にも、広く見たい、こまかく見たいという欲求もあり、データ処理能力がかぎられているという事情もある。 ここでは、空間については、上記の図と同じ 北緯 36～37.5度、東経 138.5～140.5度の範囲にふくまれるアメダス観測点(気象官署もふくむ)にかぎることにした。 時間については、2020年12月1日から19日までの日値をとりあげることにした。 (2020年12月20日現在、大雪の可能性のある状態がつづいていたので、リアルタイムにデータを更新できるとよかったのだが、 わたしの作業手順では、急にデータをさしかえるのはむずかしい。)
前回と同じく、気象庁の「過去の気象データ ダウンロード」のサイトから CSV 形式のデータを取得した。 一度にダウンロードできるデータの量の制約はあるが、時間方向の件数がすくなければ観測地点数をふやすことができる。 地点選択では都道府県別の地図の画面で地点を指定する必要があるが、複数の都道府県の地点を指定していっしょにダウンロードすることができる。 この方法で、降水量、積雪深、降雪深 [注]、日照時間、日平均気温の、 それぞれ対象とした緯度経度範囲にふくまれる地点を指定してデータをダウンロードした。 地点数は降水量がいちばんおおく、積雪深と降雪深の地点、日照時間と気温の地点は、それぞれ、降水量の地点の部分集合にあたる。
ダウンロードしたファイル(ファイル名を data.csv からかえた)を教材ウェブサイトに置く。気象庁サイトが混雑しないように、実習にはこれをつかってみてほしい。 (もしあたらしい時期のデータをくわえたければ、ダウンロードしなおすのもよいが、 そのときは、下にのべる対応表づくりもやりなおす必要が生じるだろう。)
[補足] CSVファイルの読みとばすべき行数は、日最大積雪深のファイルでは 6行、その他の4つのファイルでは 5行である。
作業を始めてみて、このダウンロードデータの形式に由来するむずかしさにぶつかった。 観測地点と、その地点のデータがふくまれるCSVファイルの列との対応を、機械的につけることがむずかしいのだ。
対策として、ダウンロードしたCSVファイルを見て、手作業で、地点名と、CSVファイルでその地点がはいっている列番号との対応表を つくっておくことにした。 ただし、「CSVファイルを見る」作業は表計算ソフトウェアでやるのが便利で、そこでは列は A～Z, AA, AB, ... とアルファベットで表示されている。 R や Python での処理では列の番号が整数でほしいのだが、表計算の列記号から番号への変換は機械にやらせることにする。 【Python の Pandas のデータフレームの列の番号は 0 からはじまるが、わたしの作業は当初 R を想定していたので、 対応表の列の番号は 1 からはじまっている。Python ではそれをつかうプログラムの側で 1 をひくことにする。】
その変換のついでに、地点表 amedas_stn.csv を読んで、地点の緯度・経度も列番号とならべて表にしておくことにする。
【対応表をつくったあと、地点名のローマ字表記がほしくなった。 カタカナ表記からの機械的変換には不満があるので、手作業で追加した。】
ここまで準備すれば、ダウンロードしたデータファイルと、key 表の2つを読んで、 日ごとの観測値の空間分布を図化することができる。
図には、観測点の位置に記号(まる)をかき、その近く (matplotlib の text 関数のdefault動作によれば右上) に数字で観測値をしめすことにした。 ただし、観測値が欠けているばあいは「NA」と表示することにした。 また、観測が不完全なばあい (「品質情報」が「8」でないばあい) は、観測値を表示する文字の色を赤にした。
図形出力先は、画面と、画像ファイル(png)とのどちらでもできるようにしたが、 ここにリンクしたプログラムは、画像ファイルに出すものにしてある。 画面に出したいときの変更点は、「plt.savefig(...)」「plt.clf()」「plt.close()」の行をコメントにする(行の頭に「#」をいれる)ことと、 コメントにしてある「plt.show()」をいかす(行の頭にいれてある「#?」などをけす)ことである。
【まえにつくった R によるプログラム例はデータの期間のはじまりの日が 2020年12月1日であることを前提としているところがあったが、 このページの Python によるプログラム例では、日付はデータファイルから読みとるようにし、 プログラムで指定するのはデータファイルに記録された最初の日を 1 とした相対番号にしたので、 データの期間がかわってもファイル名以外の変更はしなくてよいはずである。】
プログラムでは19日ぶんの図をつづけてつくるが、ここではそのうちある1日ぶんの図をのせておく。 (12月17日の降水量が「NA」となっている地点は藤原で、雪があまりに多くて雨量計が正常に動作しなかったのだ。)
対象期間の全部の日の図を一覧表のかたちでつぎのページに置いた。 ブラウザ上で縮小画像をクリックすると大きな画像が表示されるようにした。
観測地点分布図のばあいと同様に、Cartopy をつかって、正距円筒図法の地図上に 毎日の観測値の分布を表示し、海岸線・県境をかさねた。
プログラムはつぎのとおり。
ある1日の図はつぎのとおり。
対象期間の全部の日の図を一覧表のかたちでつぎのページに置いた。 ブラウザ上で縮小画像をクリックすると大きな画像が表示されるようにした。
データのあるすべての観測点についてそれぞれ、時間(ここでは日単位なので日付)を横軸、 観測値を縦軸にとったグラフをつくるプログラムを作成した。
このプログラム例は、対象の日が全部12月中であることを前提に、年月日の日を横軸にとっている。 【前に R でつくったプログラム例では、はじめの日付を 1 日にきめていたが、 今回はデータのはじめの日を読むことにした。】 複数の月にまたがるばあいは、横軸の設定のしかたを変える必要がある。
降雪深、最大積雪深、日照時間、降水量は、数値のゼロに絶対的な意味があるので、 ゼロから観測値までの棒と、観測値の位置での 〇(まる)の両方で表示することにした。 ただし、観測が不完全な(「品質情報」が 8 でない)ばあいは、まるをかかないことにした。
気温は、℃単位のゼロは約束上の原点にすぎず絶対的な意味がないので、 棒はつかわず、二重のまるで表示することにした。 ただし、観測が不完全なばあいは、大きなまるをかかないことにした。
1画面の上下に同じ横軸のグラフを2つ配置することにした。 日降雪深と日最大積雪深、日照時間と日平均気温は、同じ地点の両要素を同じ画面に配置する。 日降水量は、画面の上半分だけに図がある形になる。 【R でのプログラム例では 2地点の降水量を上下に表示することもしたが、ここではそのしかけは省略した。】
図形出力先は、画面と、画像ファイル(png)とのどちらでもできるようにしたが、 ここにリンクしたプログラムは、画像ファイルに出すものにしてある。 画面に出したいときの変更点は、「plt.savefig(...)」「plt.clf()」「plt.close()」の行をコメントにする(行の頭に「#」をいれる)ことと、 コメントにしてある「plt.show()」をいかす(行の頭にいれてある「#?」などをけす)ことである。
画像ファイル名に 日本語文字(漢字など)がはいると OS によってはただしく動作しないので、 ファイル名のなかの地点名は ASCIIコードのローマ字ではいるようにした。 【 R のプログラム例ではここまでだったが、この Python のプログラム例では】さらに、 OS によっては大文字・小文字だけがちがうファイルが同一視されたりされなかったりするので、 念のため、ファイル名のなかの地点名はローマ字の小文字に統一した。
プログラムではデータのある全地点の図をつづけてつくる。 ここではそのうち1地点の図をのせておく。藤原の降水量にはデータ欠損や不完全な観測が見られる。
全地点の図を一覧表のかたちでつぎのページに置いた。 ブラウザ上で縮小画像をクリックすると大きな画像が表示されるようにした。