[計算機を利用した気象分野の実習の教材]

格子点型の気象データ (JRA-55 地上風) の地図上表示 (Pythonによる試作)

[格子点型の気象データ (JRA-55 海面気圧) の地図上表示 (Pythonによる試作)]のつづき。

格子点型の気象データによるベクトル場の表示の例として、 風速ベクトルを矢印で地図上に表示する。


JRA-55の地上風

JRA-55 の anl_surf125 には、 地上 10 m (海抜ではないことに注意) の風速の 東西成分 (西風が正) と 南北成分 (南風が正) がふくまれている。 anl_surf125.ctl に記述された変数名は、ugrd10m と vgrd10m、 単位は m/s である。

風の東西成分と南北成分を計算機であつかうばあいの 変数名や変数記号は、かならずしも統一されていない。 数式で書く場合の変数記号に u, v をつかうことが多いので、 ここでは便宜上「u」「v」としておく。 【地上風と等圧面の風が両方でてくるばあいや、 再解析の風と観測値の風が両方でてくるばあいには、区別した変数名をつける。】

【なお、気象の研究のためのデータによくつかわれている netCDF というデータ形式でデータを記録するときの 変数名などをきめた CF規約 (CF Convention) の標準変数名は、 x_wind または eastward_wind と、y_wind または northward_wind になっている。 (Version 77 (2021-01-19) https://cfconventions.org/Data/cf-standard-names/77/build/cf-standard-name-table.html )。 CF規約は、気象学以外の分野にもつかわれており、とくに海洋学との用語の共通化が 考えられたので、このような表現になったのだと思う。】

ここでしめすプログラム例は、 [JRA-55のGRIBデータをGrADSで読み、単純バイナリデータに変換する]によって 変換ずみの単純バイナリデータを読む。

データの例として、第2回試作の図の例に対応する、2007年 8月 16日 世界時 6時 (日本標準時 15時)の u と v の数値データをリンクしておく。

これからしめすプログラム例では、u と v のデータが置いてあるディレクトリ名を、 それぞれ datapathu, datapathv という文字列変数で指定している。 利用者は自分がつかう計算機上のデータ置き場に応じてその変数の値をかきなおしてつかうとよい。


第1回試作 (2021-07-09)

Python の matplotlib にある quiver という関数によって、 2次元平面上の多数の位置にあるベクトルをまとめて矢印で表示することができる。 それを Cartopy の地図投影とくみあわせて、地図上に表示することができる。

ここでは、[格子点型の気象データ (JRA-55 海面気圧) の地図上表示] で 海面気圧のデータを読むかわりに u, v のデータを読み、 「日本付近の海面気圧の分布図 (等値線)」と同様に、 北緯 20~50度、東経 120~150度の領域について、 quiver によって風速ベクトルの矢印をかくプログラムをつくった。

quiver は、ベクトル量を2次元配列としてあたえるばあい、 位置座標も同じ寸法の2次元配列としてあたえるのが原則である。 ここではそうしないで、経度、緯度をそれぞれ1次元配列であたえている。

u, vの配列全体をあたえると、つぎの警告が出る。

cartopy\mpl\geoaxes.py:1841 UserWarning: Some vectors at source domain corners may not have been transformed correctly
u, v = self.projection.transform_vectors(t, x, y, u, v)

これは、北極・南極の格子点の風の地図投影変換がさだまらないということだとわかった。 quiver にあたえる配列から、北極・南極にあたる行をはずしておけばよい。 それには、lat, u, v のところをそれぞれ「lat[1:-1]」, 「u[1:-1,:]」, 「v[1:-1,:]」とすればよい。 (つぎのウェブ記事を参考にした。)

https://qiita.com/convection/items/6d5fa546d81fd5c9b452
@convection が 2020年02月16日に作成
cartopyとmatplotlibを使って緯度経度座標のデータ(気象データとかを想定)の可視化

図を読む人が風速を読みとれるためには、凡例となる矢印をしめす必要がある。 そのために quiverkey という関数があるので、それをつかって 10 m/s の矢印をしめした。 ところが、quiverkey にあたえる矢印の長さの数値を 10 にしたのでは、結果の矢印が短すぎる。 わたしは、quiverkey で矢印の長さをどう決めているか、まだ理解できていない。 quiver で矢印の長さをどう決めているかも、あまりよく理解できていない。 quiverkey には、ひとまず目分量で、つぎにしめす風速の等値線図を参考に、14という数値をあたえてある。

プログラムは、 同じ年のうちで、月、日、時の3重ループで、つづけて図をつくるように構成した。 ただし、現状では、1つの時刻 (2007年 1月 7日 世界時 0時) の図をつくるようになっている。 3つの for 文を それぞれ コメントにしてあるものと入れかえると、 複数の月をつづけてあつかうこと、月の全部の日について作図すること、毎日4回の時刻について作図することなどが できるようににしてある。

プログラムの終わり近くの「plt.show()」を生かすと、図は画面に出る。 「plt.savefig(...)」をふくむ4行を生かすと、図は画像ファイルに保存される。

風の東西成分・南北成分から風速の大きさを計算して、その等値線をかさねてみた。

↓ 結果の例 wind2007010700_contour_japan_v1.png (2007年 1月 7日 世界時0時の地上10m風)

wind2007010700_contour_japan_v1.png

結果をみると、陸上の風速が、海上にくらべてけたちがいに弱い。 実際に陸面のほうが摩擦が大きくて風速が弱いのだろうが、JRA-55 の「地上 10 m の風速」が陸上で弱すぎる可能性もあると思う。 観測値と比較して評価してみることにしたい。


第2回試作 (2021-07-24)

[気象観測データによる風の空間分布の作図]の第2回試作ができたので、 JRA-55の地上風を、それとかさねて比較できるように作図することにした。

観測データによる作図では、1地点ごとに quiver を呼び出して矢印をかいた。 凡例の矢印も、同様に ひとつの (u, v) の組をあたえて quiver を呼び出してかいた。 今回は、JRA-55 の格子点データによる作図でも、1地点ごとに quiver を呼び出して矢印をかく方法をとることによって、 風速と矢印の長さとの関係を統一することにした。 【quiver 関数の本すじのつかいかたは (u, v) を 2次元配列であたえるやりかただと思うが、 わたしにはそのばあいの矢印の長さや太さなどのきまりかたがよくわからないので、ひとまずそれをさけた。】

日本付近

第1回と同じく、北緯 20~50度、東経 120~150度の領域について作図した。 速度ベクトルの大きさの等値線をかさねたプログラム例をしめす。 (等値線をかさねないばあいは、 コメント「# 等値線をかく」から コメント「# タイトルをかく」の前までをとりのぞけばよい。)

↓ 結果の例 wind2007081606_contour_japan_v2.png (2007年 8月 16日 世界時 6時の地上10m風)

wind2007081606_contour_japan_v2.png

関東地方付近

北緯 34.5~38度、東経 138~142度の領域について作図した。 速度ベクトルの大きさの等値線をかさねたプログラム例をしめす。

↓ 結果の例 wind2007081606_contour_kanto_v2.png (2007年 8月 16日 世界時 6時の地上10m風)

wind2007081606_contour_kanto_v2.png

観測値と再解析の値とを同じ地図にしめす

観測値と再解析の値とを同じ地図にしめしたい。

そのためには、データの時刻をあわせる必要がある。 気象庁のウェブサイトからダウンロードした風速の観測値のデータは、1時間ごとで、日本標準時 (JST) の1時からはじまっている。 JRA-55は、6時間ごとで、世界時 (UTC) の 0時からはじまっている。 Python の datetime パッケージを利用して、JST と UTC の変換をする。

まず、[気象観測データによる風の空間分布の作図]の プログラム例 map_hourly_wind_arrow_kanto.py を変形して、 観測値の作図の対象期間である 2007年 8月 11-20日の、 日本時間 3, 9, 15, 21 時だけをえらんで、観測値の風を矢印で作図し、 その表題に JST と UTC の両方の日時を表示するプログラムをつくった。

それと、上の wind_vector_contour_map_kanto_6hourly_v2.py とをくみあわせて、 観測値と再解析値の両方を表示するプログラムをつくった。

↓ 結果の例 map_wind_obs_jra_2007081615jst.png (2007年 8月 16日 日本標準時 15時 (世界時 6時) の地上10m風)。 緑が観測値、明るい青がJRA-55再解析 (1.25度格子)。

map_wind_obs_jra_2007081615jst.png

図にしめした時刻には、観測値によれば、埼玉県付近で風向がさまざまで風速 5 m/s ぐらいの風がふいていた。 それが、JRA-55再解析では、風が弱いように表現されていた。


2021-07-09, 2021-07-24
増田 耕一 (MASUDA Kooiti)
[計算機を利用した気象分野の実習の教材] [立正大学用教材目次] [教材目次]