観測点ごとの気象データを地図上の観測点の位置に表示する作業のつづきとして、 風のデータをあつかってみる。 材料は前回の [観測点ごとの気象データを図示する] とおなじく 気象庁のウェブサイトからダウンロードしたデータだが、 風は日平均風速だけではおもしろくないので、時別値 (1時間ごとの観測値) をつかう。 風を矢印や矢羽根で表示するには、格子点データのばあいと同様に、 Cartopy の quiver や barbs をつかう。 (この機能は格子点データむきにつくられているので、 観測点ごとのデータに最適ではないようだが。)
今回のPythonプログラムの例文とそれがつかうデータは、luce の /home/masuda/mapstationwind に置いてある。
各自、つぎのコマンドによって、プログラム例をコピーしてほしい (最後の空白のあとのピリオドが重要)。
cp -r /home/masuda/mapstationwind .
そして、mapstationwind に cd して作業してほしい。
なお、前回の授業の材料のデータを読むことがある。 プログラム例では、 前回の授業のディレクトリ mapstation が mapstationwind の横ならび (同じ親ディレクトリの下) にあることを想定しているところがある。 もし別の場所にあるばあいは、 プログラム中で read_csv, read_table などをよびだしているところに 書かれているファイル名の「../mapstation/」という部分を修正してほしい。
気象庁ウェブサイトの「各種データ・資料」のページの「気象」のうちの「気象観測データ」のところにある 「過去の地点気象データ・ダウンロード」から、 2022年 7月 18日に取得した。
地点の選択は、指定順序もふくめて、前回と同じにした。
「項目を選ぶ」で、データの時間間隔を「時別値」とし、「風向・風速」だけをえらんだ。
「期間を選ぶ」では、2022年 7月 1日から 2022年 7月 10日 まで (7月上旬) を指定した。 【6月上旬、6月中旬、6月下旬、7月中旬 (7月21日に取得) のデータも 同じ形式でダウンロードしたので、下にリンクしておく。 ただし、luce の mapstationwind には 7月上旬のデータだけを置いた。】
「表示オプションを選ぶ」のところは、いずれも、あらかじめ与えられているほう (デフォールト値) にしておく。 その「ダウンロード CSV ファイルのデータ仕様」のうち「日付の形式」は「日付リテラルで格納」である。
7月2日から4日のデータには、 気象庁が契約していた au (KDDI) の通信網の故障によるデータの欠損が生じていて、 7月24日現在、まだじゅうぶんに回復されていないようだ。
【その後、回収が進んだようなので、ダウンロードしなおして、別のファイル名で持つことにした。[2022-08-26]】
ダウンロードされた data.csv のファイル名をつけなおし、つぎのところに置いた。
【[2022-07-27 補足, 2022-08-26 改訂] luce では、wind_kanto_hourly_20220701_20220710.csv だけを /home/masuda/mapstatonwind の下に置き、 その他のデータファイルは /home/masuda/jmaobs の下に置いたので、 関心のある人はさしかえてつかってみてほしい。】
csv ファイルのなかみをみると、6行の見出し情報のあと、7行め (Python流にかぞえた行番号は 6) からの各行がそれぞれ 1つの時点 (このばあいは 1時間) のデータとなっている。
各列の内容はつぎのようになっている。
地点数が 40 なので、列の数は 1 + 5 × 40 で、201 である。 表計算ソフトウェアに読みこむと、列記号が A, ... Z, AA, AB, ... と進んで「GS」にいたる。
前回の [27b] [28a] [28b] と同じひろがりの地図上の観測点の位置に、風速の観測値を数字でしめす。
どの時刻をとりあげるかは、データフレーム dfobs の行番号を 「krecstart = 53」 「krecend = 65」 のように指定するようにしてある。 この行番号は、観測データのある最初の時刻 (今回のデータセットでは 7月 1日 日本標準時 01時 [00時ではない]) を 0 としている。 【datetime パッケージの関数をつかって日時と行番号との換算をすればよいのだが、 まだくみこんでいない。】
例文に書きこんだ 53 は 7月 3日 日本時間 6 時、65 は その日の 18 時である。 【セミナーの教材としてこの日時をえらぶ必然性はないのだが、 教員3人が1年生の巡検で秩父に行った日の状況を見たかったのと、 データ欠損のはいりかたを見たかったので、この日にした。】
プログラムの中では、まず jma_observatories.csv と jma_stations2.csv を読み、データフレームを連結して dfstn とする。 さらに、station_list_kanto.txt も読み、「関東領域」内の地点だけのデータフレーム df をつくる。
また、観測値の csv ファイルを データフレーム dfobs に読みこむ。
表示する時刻ごとに、観測値のリスト val と、品質番号のリスト qval を作成する。 (前回の教材では、このふたつのリストを、df の列として追加してつかったが、 プログラムの流れを見なおした結果、df の列にする必要はないので、リストのままつかうことにした。) val に、dfobs の、krec で指定された行の、2 から 5 ずつふえていく番号で指定された列の値を入れていく。 qval に、そのつぎの列の内容を入れていく。
ここから作図。だいたい [28b] と同じ。 ただし、地点の上側に表示する文字列が、風速である。 品質情報のあつかいを、前回の気温などのばあいと同様に、つぎのようにした。
気象庁からダウンロードした CSV ファイルには、 風向は 16方位で、「西北西」のように日本語文字 (文字コード Shift_JIS) で記録されている。 luce の Python で図にふくめられることが確実な文字は ASCII コードにかぎられるので、 風向を「WNW」のような英語の頭文字におきかえて、各地点の上側に表示する。
風速の 16 方位には、北北東を 1 として、東のほうにまわる順に番号がつけられている。 静穏を 0 とするので、北風は 16 とされている。
気象学では、x 軸 を東西方向 (東が正)、y 軸 を 南北方向 (北が正) にとるのがふつうである。 風速ベクトルは {u = dx/dt, v = dy/dt} の形になる。 u は 風の東西成分で、西風 (東にむかう流れ) が正、 v は 風の南北成分で、南風 (北にむかう流れ) が正である。
2次元ベクトルの直交直線座標と極座標とのあいだの変換はつぎのようにかける。
「風速・風向」から「(u, v)」へ変換するとき、 ベクトルの大きさは、風速をそのままつかえばよいのだが、 偏角と風向との関係は、つぎのことに注意する必要がある。
プログラムでは、つぎのような換算式をつかっている。 (4が直角、8が平角 (2直角) にあたる。)
arg = (-4 -idir) * math.pi / 8
u[kstn] = sp * math.cos(arg)
v[kstn] = sp * math.sin(arg)
静穏のばあいは別あつかいにして、u と v の両方に 0 をいれることにした。 (静穏はかならずしも風速が厳密に 0 ではないのだが、風向がわからないとき、 風速ベクトルの東西・南北成分に 0 以外の値をあたえる理由がない。)
地点の左上に u、右上に v の値を数字でかくようにしてみた。
格子点データの風の表示でもつかった quiver によって、(u, v) の値にもとづいて矢印をかく。
quiver は複数の矢印をいっぺんにかけるのだが、 ここでは 1地点ごとによびだしている。
【しかも、quiver にあたえる x, y と u, v のデータの型がちがうとうまくいかず、 試行錯誤のすえ、x, y, u, v ともに、要素が 1個だけ (0番だけ) の配列変数にいれてから わたすことにした。 ここは熟練者ならばもっとうまくやれるだろうと思う。】
風速と図上の矢印の長さとの関係は、quiver のパラメータ scale によっている。 scale の意味はまだよくわからないが、試行錯誤で 50 ときめた。 凡例として、風速の値をあたえた矢印を (図の左上に黒で) しめすことにした。
静穏のばあいは、矢印でなく まる をかくことにした。 他方、静穏ではないのだが風速が弱いばあいは、矢印をかくのだが、 小さすぎて読みとれないこともある。
quiver のかわりに barbs をつかうと、風を天気図でつかわれる矢羽根であらわすことができる。
矢羽根の線の数は、国際式天気図のきまりがあるが、それは高層天気図むきにつくられていて、 地上、とくに陸上の風ではわずかな数しかかかれない。 ここでは、つぎのように変更した。 (ただし、まだ凡例を図にふくめていない。)
[30f] の観測点の風の図に、 7月12日の授業 [格子点型の気象データによる風の矢印 (および、矢羽根) 表示] でやったのと同様な JRA-55 格子点データの地上 (10 m) 風をいっしょに表示してみる。 (観測点は「〇」印、格子点は「×」印、 観測点の風の矢印は濃い緑、格子点の風の矢印は青、と区別する。)
今回のプログラムの、格子点データの風を矢印で表示する部分は、 7月12日の教材の [25a] とほぼ同じである。 ただし、風速の等値線をかくことはしない。
luce に置いた JRA-55 のデータは、7月12日の授業のときは 2021年の 1年ぶんだったが、 その後、2022年 7月 21日のぶんまで追加したので、今回の観測点データと同じ期間がふくまれている。 ただし、JRA-55 データの時間間隔は 6時間で、世界時 [UTC] の 0, 6, 12, 18時 (日本標準時 [JST] では 9, 15, 21, 3時) である。
いずれも日本の気象庁が提供するデータだが、 日時が、日本語ウェブサイトから提供されている観測データでは日本標準時 (JST)、 JRA-55 では 世界時 (UTC) でしめされている。 (JRA-55 は全球データだから当然だが、対象地域が日本国内でも、 専門知識のある人むけのデータ提供では、日時が UTC でしめされていることが多い。)
7月12日のプログラム例 [25a] では、日時を UTC であつかい、タイトルの表示のところだけ JST への換算を追加した。 今回は、日時を JST であつかっている [30f] をもとにしたので、 JRA-55 のファイルを読むところだけ、UTC に換算した日時をつかうことにする。
また、変数名の衝突をさけるため、格子点データに関する変数名を u, v, lon, lat から ugrid, vgrid, longrid, latgrid に変えた。
プログラム例の現状では、「krecstart = 50」 「krecend = 68」として、 2022年 7月 3日の日本時間 3時から 21時まで、6時間おきに作図するようにしてある。
例文では、地上気象観測のある領域について JRA55 の値も表示したが、 JRA55 は地上気象観測のないところにも値があるので、 海上をふくめた気象状況をつかむのに有用である。
この例文では、「# 緯度経度」と注記したところの lonmin, lonmax, latmin, latmax にあたえた数値を変えれば、地図の領域を変えられる。 ただし、領域をひろげると、軸の目盛りにつけた文字がこみあいすぎることがある。 そのときは、目盛りの間隔を指定する dlon, dlat の値を大きくするとよい。 なお、lolaeps は、数値範囲の上限をちょっとこえる値を指定するために つかっている変数で、dlon, dlat がこれよりも小さくならなければ変更しなくてよいだろう。
1つの地点について、 画面を 4段にわけ、それぞれ横軸に日時をとって 風向 (16方位)、風速 (m/s)、東西成分 (m/s)、南北成分 (m/s) の 点グラフをかくプログラムを作成した。
地点の指定は、変数 kstn の値として、 station_list_kanto.txt の行番号 (ただし 0 はじまり) をあたえることによる。 現状では「kstn = 11」で、熊谷のグラフがつくられる。
(今回は、データ欠損 または 品質情報が 2以下のところは、欠損として処理した。 しかし、欠損箇所を明示する手つづきはくみこまなかったので、 欠損箇所はただ見えないだけになっている。)
前回と今回の関東領域の地図をかく例文では、県境の線をかくようにしたのだが、 授業のときにつかった例文プログラムでは、ユーザーによっては県境が表示されなかった。 (県境のデータが自動的にダウンロードされなかったようだ。)
【2022-08-03 から 2022-08-26 までに、 [日本の都道府県境をかく / 中国の「省」の境をかく]の ページに書いた方法をつかうように例文を変更した。 区別のため、改訂版プログラムのファイル名に「_v20220803」をふくめている。 [この部分、2022-08-26 改訂]】
なお、県境が不要ならば (あるいは、別の方法で県境をかくならば)、つぎの部分は不要となるので、消すかコメントにしたほうがよいだろう。