風の速度は、方向と大きさをもつベクトル量である。 その空間分布を表示するには、矢印や、天気図の矢羽根などがつかわれる。 ここでは、それを Python で作図してみよう。 なお、風のデータとしては、観測点での観測値も、 データ同化による格子点データもつかわれる。 両者を比較できる形でいっしょに表示することをこころみる。
Python で Cartopy パッケージをつかうと、 quiver という関数で風ベクトルを矢印で地図上に表示できる。 残念ながら、まだわたしには quiver で矢印の長さや太さなどを思うように決められない。 わたしの quiver についての理解が不足しているのかもしれないし、 Cartopy が開発途上でまだ機能が不足しているのかもしれない。 ひとまず、現時点でよさそうだと思われる方法を紹介するが、 この方法が模範だと思わないで、よりよい方法をめざしてほしい。
前回の[格子点型の気象データ (JRA-55 海面気圧) の地図上表示]と同様に、 Pythonプログラムの例文は、luce の /home/masuda/shared/program/ に置いてある。
各自、自分のホームディレクトリの下の作業用ディレクトリ (python_map_practice をまた使ってもよい) に cd してから、 つぎのコマンドによって、プログラム例をそこにコピーしてほしい (最後の空白のあとのピリオドが重要)。
$ cp /home/masuda/shared/program/2*wind*.py .
JRA-55のデータは、前回と同様、/home/masuda/shared/jra55/surf125/ に 置いてあるものを読むので、自分のところにコピーする必要はない。 (もし、別のところに置かれたデータを読むばあいは、 プログラムのなかの datapathu, datapathv という変数にあたえる文字列を変更する。)
気象庁からダウンロードした観測データと、それを読む補助としてつくった「key表」は、 /home/masuda/shared/jma_obs/ に置いてあるものを読むので、自分のところにコピーする必要はない。 (もし、別のところに置かれたデータを読むばあいは、 プログラムのなかの datapathobs という変数にあたえる文字列を変更する。 とくに、データファイルを作業ディレクトリに持った場合は、 「datapathobs = ''」 [右辺は single quotation mark を2つ続けて 空文字列を表現している] とすればよい。)
プログラム例 23～25では、都道府県境の線をひいている。 そのためのデータは、「Python による気象解析入門」の「19. 地図を描いてみよう」の スクリプト 15_japan_map_with_boundaries.py でつかったのと同じ GADM の shape ファイルである。
$ ln -s /home/share_data/python_map_practice/gadm36_JPN_shp .
前回と同じだが、 作図を実行するには、あらかじめ Xming を起動してから Teraterm で luce に login し、 作業用ディレクトリに cd して、 プログラムを必要に応じて編集したあと、つぎのような形で実行する。
$ python3 なになに.py
図を画面に表示したいときは「plt.show()」を、 画像ファイルに保存したいときは「plt.savefig(...)」をふくむ4行を、それぞれ生かし、 つかわないほうは行の頭に「#」をつけてコメントにしておく。 例文のその部分が自分の意図とあっているか確認し、ちがっていたら編集で変更してから実行してほしい。 画像ファイルは、作業用ディレクトリの下の「image」というディレクトリの中につくられる。 画像ファイルに保存したいときは、あらかじめ「mkdir image」で そのディレクトリをつくっておいてから実行する。
風の東西成分と南北成分を計算機であつかうばあいの変数名や変数記号は、かならずしも統一されていない。 数式で書く場合の変数記号に u, v をつかうことが多いので、 ここでは便宜上「u」「v」としておく。 格子 (grid) データであることを明示して「ugrd」「vgrd」とすることもある。
JRA-55 の anl_surf125 には、地上 10 m (海抜ではないことに注意) の風速の 東西成分 (西風が正) と 南北成分 (南風が正) がふくまれている。 anl_surf125.ctl に記述された変数名は、ugrd10m と vgrd10m、 単位は m/s である。
ここでしめすプログラム例は、 [JRA-55のGRIBデータをGrADSで読み、単純バイナリデータに変換する]によって 変換ずみの単純バイナリデータを読む。
ここでしめすプログラム例では、u と v のデータが置いてあるディレクトリ名を、 それぞれ datapathu, datapathv という文字列変数で指定している。 今回の実習では、「/home/masuda/shared/jra55/surf125/ 」と設定してあり、 そのままつかえばよい。
データの例として、2007年 8月 16日 世界時 6時 (日本標準時 15時)の u と v の数値データをリンクしておく。
Python の matplotlib には quiver という関数がある。 これを Cartopy とくみあわせてつかうことによって、 風速を矢印で地図上に表現することができる。 凡例の矢印をしめすための quiverkey という関数も用意されている。
(quiver ということばのもともとの意味は、「矢づつ」らしい。 「振動」のような意味もあるが、ここではそうではない。)
quiver には、水平2次元ベクトルの東西成分と南北成分の数値をあたえる。 ここではそれぞれ、u, v としておく。
quiver は、複数の地点の 経度, 緯度, u, v の値を配列でまとめてあたえられて、 そのすべての点に矢印をかくことができる。 格子点データを作図するばあいは、 u, v は 2次元配列で用意し、quiver の一度の呼び出しで作図するのが のぞましいつかいかたとされている。 わたし(増田) は、そのやりかたで作図することはいちおうできたのだが、 スケーリング (ここでは、風速の大きさと地図上の矢印の長さの関係) をきめるための scale と scale_units をどう指定すれば図にどのようにあらわれるのか、 まだよくわからない。 さらに、quiverkey でかいた凡例の矢印が、期待した長さにならなかった。 (状況は[別記事]の「第1回試作」に書いた。)
今回のプログラム例では、 格子点データも観測点データも 1地点ごとに quiver を呼び出して同じスケーリングで矢印をかく。 凡例の矢印も quiverkey ではなく quiver でつくり、text で注記をつけることにする。
前回の「海面気圧の等値線図 (日本付近)」(プログラム例 19) と同じく、北緯 20～50度、東経 120～150度の領域について作図した。 速度ベクトルの大きさ (u, v からプログラム内で計算した) の等値線 (contour による) をかさねたプログラム例をしめす。 (等値線をかかないようにするには、 コメント「# 等値線をかく」から コメント「# タイトルをかく」の前までをとりのぞけばよい。)
スケーリングのための quiver への指定は、「scale=200」で、scale_unitsは指定しない。 (図が見やすいように数値を試行錯誤で調節したもので、数値の意味はよくわかっていない。)
プログラムは、前回の 18b, 18c, 19, 20 と同様に、年は固定して、月、日、時の 3重ループで 図をつくることをくりかえすように書いてある。 現状では、2007年 8月 11日 から 20日 の毎日の 世界時 (UTC) 0時 の 図をつくるように なっている。
↓ 結果の例 wind2007081606_jra_japan.png (2007年 8月 16日 世界時 6時の地上10m風)
前回の「海面気圧の等値線と格子点での数値の表示 (関東地方付近)」 (プログラム例 20) と同じく、 北緯 34.5～38度、東経 138～142度の領域について作図した。 速度ベクトルの大きさの等値線をかさねたプログラム例をしめす。
スケーリングのための quiver への指定は、「scale=50」で、scale_unitsは指定しない。
日時についてのループは、上の「日本付近」のばあいと同じにしてある。 日数をしぼって1日4回にかえてみるとよいだろう。
↓ 結果の例 wind2007081606_jra_kanto.png (2007年 8月 16日 世界時 6時の地上10m風)
観測値の例として、 (熊谷で最高気温の極値が出た 2007年 8月 16日 日本標準時 (JST) 14時をふくむ) 2007年 8月 11日から20日までの期間の、 関東7都県と福島、新潟、長野、静岡、山梨の各県にある 気象官署や特別地域気象観測所 (ただし北緯34.5度よりも南の島と富士山をのぞく) の 風向・風速の時別値を、 気象庁ウェブサイトの「過去の気象データ ダウンロード」から 2021年7月23日にCSVファイルとしてダウンロードした。
左からかぞえて最初の列に、日時が、(気象庁のいう) 「日付リテラル」の形式ではいっている。 それから、1観測点ごとに5列をつかって、観測値がはいっている。 (5列は、風速、風速の品質情報、風向、風向の品質情報、風速と風向の均質番号である。) このファイルは、日本語文字を Shift_JIS 文字コードでふくんでいること、 そのまま Pandas の列名につかえる見出しをふくむ行がないこと、など、不便ではあるが、 手作業で編集するのもてまがかかるので、 そのまま読むようにプログラムを組むことにする。
観測地点がどのような順序で収録されているかのリストを手作業でつくった。 (ダウンロードの作業で指定した順番だが、おぼえていない。 CSVファイルを表計算ソフトウェアで開いて、順に見て、記録した。 風向・風速のばあい、どの観測点も列数は同じなので、出現順序さえわかれば列番号はわかる。)
気象庁がダウンロードページとは別のページで公開しているアメダス地点表 (気象官署もアメダス地点でもある) から、各地点の緯度・経度を知り、 地点名、緯度・経度、データファイル中の列番号の対応表 (key表とよぶことにする) をつくった。 こちらは、最初の行に、Pandasの列名としてつかう文字列を入れた。 (そこまでの手順に興味のある人は [別ページ]の 「第2回試作」のところを見てほしい。)
ここからしめす、観測点の風のデータを作図するプログラムのそれぞれで、 観測データを読んだら、この処理をおこなう。
気象庁のサイトからダウンロードしたファイルには、 風向が、たとえば「北北東」のような日本語文字列で収録されている。 これを、静穏を 0、北北東を 1、北東を 2、... 北を 16 という 16方位のコード番号 (整数値、変数名 idir とする) におきかえる。
風速は m/s 単位で収録されているので、単位換算の必要はない。 しかし、quiver がほしいのはベクトルの東西成分と南北成分だから、 風向・風速から u, v への変換が必要だ。
まず、風向を、東向きを x 方向、北向きを y 方向としたベクトルの偏角 (変数名 arg とする) に換算する。 風向コードと偏角とは逆まわりであり、しかも、北風は南向きの流れなので、まちがえやすい。
どうやら、次の計算式でよいようだ。
arg = (-4 - idir) * math.pi/8
(もし風向が「度」単位の数値であたえられているならば、つぎのようになる。)
arg = (-90 - idir) * math.pi/180
そして、風速に cos(arg) をかけると風の東西成分、sin(arg)をかけると南北成分がえられる。
プログラム例 22 と同じ関東地方付近の地図に、観測点を plot によって 〇印でかき、 風速ベクトルを quiver によって矢印で表示する。
風速の値と矢印の長さとの関係は、 quiver に、「scale=50」とあたえている。scale_units は指定しない。 静穏 (風速 0.2 m/s 以下) のばあいは地点を ○ でかこむ (quiver と別に plot をよびだしている)。 凡例の矢印を、quiverkey をつかわず、quiver に 同じ scale をあたえてかいている。
時間についてのループは一重にし、 対象としたデータの範囲の時刻 (1時間きざみ) の とおし番号 (1はじまり) を kstepsel として、その値で時刻を指定する形にしてある。 日時の情報は、観測データのCSVファイルの左端の列の日付リテラルを Python の datetime によって読みだしている。
現状のプログラムの kstepsel の範囲は 121 から 144 までで、これは 8月16日の日本時間1時から24時 (翌日0時) までの24枚の図をつくるものである。
↓ 結果の例 wind2007081614_obs_kanto.png (2007年 8月 16日 日本標準時 14時の地上10m風)
【試作の過程で、(u, v)の値、風向、風速を、 それぞれ文字で表示することもこころみたが、現状のプログラムではコメントにしてある。 そのうち風向の表示は、漢字で、MS Windowsの日本語フォントを前提としているので、 luce ではそのままでは動かないだろう。】
風速ベクトルを天気図の矢羽根で表示する。 quiver のかわりに barbs をつかう。 風速の値と矢羽根との関係は、 短い羽根が 1 m/s、長い羽根が 2 m/s、旗が 10 m/s とした。 (国際式天気図の方式をまねているが、数値が同じでないことに注意。 陸上の地上風は上空や海上の風にくらべて弱いので、同じ尺度だとよく見えてこない。)
↓ 結果の例 wind2007081614_obs_barb_kanto.png (2007年 8月 16日 日本標準時 14時の地上10m風)
観測値と再解析の値とを同じ地図にしめしたい。
そのためには、データの時刻をあわせる必要がある。 気象庁のウェブサイトからダウンロードした風速の観測値のデータは、1時間ごとで、日本標準時 (JST) の1時からはじまっている。 JRA-55は、6時間ごとで、世界時 (UTC) の 0時からはじまっている。 Python の datetime パッケージを利用して、JST と UTC の変換をする。
観測値のデータを用意してある 2007年 8月 11-20日のうちから日をえらび (現状では16日をえらんでいるが、他の日あるいは連続した複数日に変更可能)、 その日の日本標準時 3, 9, 15, 21 時 (UTCでは前日の18時、当日の 0, 6, 12時) について、 観測点の位置に〇、格子点の位置に + 印をかき、 観測値の風とJRA-55の風とを色をかえた矢印で作図し、 その表題に JST と UTC の両方の日時を表示するプログラムをつくった。 画像ファイル名には JST で日時を入れたので「jst」という文字もいれておいた。
↓ 結果の例 wind2007081615jst_obs_jra.png (2007年 8月 16日 日本標準時 15時 (世界時 6時) の地上10m風)。 緑が観測値、明るい青がJRA-55再解析 (1.25度格子)。
図にしめした時刻には、観測値によれば、埼玉県付近で風向がさまざまで風速 5 m/s ぐらいの風がふいていた。 それが、JRA-55再解析では、風が弱いように表現されていた。