PythonからGMT6.xを使って天気図を描く

天文と気象（とその他いろいろ）

2020年7月23日 15:34

はじめに

前回、PythonからGMT6.xを使って、天気図を作成するための環境を作った。

今回は、GSM/MSMのGPVデータから天気図を作成してみる。

ーーーーー追記ーーーーー
現在はGMTを使わずに天気図を描いています。

過去の経緯はマガジンで紹介しています。

ーーーーー追記おわりーーーーー

GPVのデータを読み取る

下記のサイトを参考にさせてもらった。

GPVファイルをオープンする。以下、GSMを例に示す。

>>> import pygrib as grib
>>> gpv = grib.open('Z__C_RJTD_20200704120000_GSM_GPV_Rgl_FD0000_grib2.bin')

gpv.message(int num)メソッドを使うと、GSMのGPVデータには、順に地上の海面更正気圧、現地気圧・・・が含まれていることが分かる。

>>> gpv.message(1)
1:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 202007041200
>>> gpv.message(2)
2:Surface pressure:Pa (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202007041200

これはwgrib2コマンドで確認できる項目と同じだ。合計103種類のデータを含んでいる。

username@Mint-Cinnamon:~$ wgrib2 Z__C_RJTD_20200704120000_GSM_GPV_Rgl_FD0000_grib2.bin
1.1:0:d=2020070412:PRMSL:mean sea level:anl:
1.2:0:d=2020070412:PRES:surface:anl:
・・・・・・
1.102:0:d=2020070412:TMP:10 mb:anl:
1.103:0:d=2020070412:VVEL:10 mb:anl:

gpv.select(name = 'name')を使うと、nameをキーワードに検索することができる。以下では、高度場で検索してる。

>>> gpv.select(name = 'Geopotential Height')

検索された項目がリスト形式で取得できる。

また、gpv.select(shortName = 'shortname')でも同様に検索できる。

>>> gpv.select(shortName = 'gh')

説明が後になったが、nameやshortNameというのは、下記に定義されているので参照して欲しい。

複数の条件で検索することもできる。以下では、500hPaの気温を検索する。

>>> gpv.select(shortName = 't', level = 500)
[44:Temperature:K (instant):regular_ll:isobaricInhPa:level 50000.0 Pa:fcst time 0 hrs:from 202007041200]

次に、GPVデータを見る。

>>> list = gpv.select(shortName = 't', level = 500)
>>> list[0]
44:Temperature:K (instant):regular_ll:isobaricInhPa:level 50000.0 Pa:fcst time 0 hrs:from 202007041200
>>> list[0].data()
(array([[257.09692383, 257.09692383, 257.09692383, ..., 257.09692383,
       257.09692383, 257.09692383],
      [256.83129883, 256.83129883, 256.83129883, ..., 256.83129883,
       256.83129883, 256.83129883],
      [257.01879883, 257.01879883, 257.01879883, ..., 257.00317383,
       257.01879883, 257.01879883],
      ...,
      [230.06567383, 230.06567383, 230.06567383, ..., 230.08129883,
       230.06567383, 230.06567383],
      [230.69067383, 230.69067383, 230.69067383, ..., 230.69067383,
       230.69067383, 230.69067383],
      [231.37817383, 231.37817383, 231.37817383, ..., 231.37817383,
       231.37817383, 231.37817383]]), array([[ 90. ,  90. ,  90. , ...,  90. ,  90. ,  90. ],
      [ 89.5,  89.5,  89.5, ...,  89.5,  89.5,  89.5],
      [ 89. ,  89. ,  89. , ...,  89. ,  89. ,  89. ],
      ...,
      [-89. , -89. , -89. , ..., -89. , -89. , -89. ],
      [-89.5, -89.5, -89.5, ..., -89.5, -89.5, -89.5],
      [-90. , -90. , -90. , ..., -90. , -90. , -90. ]]), array([[  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5],
      [  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5],
      [  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5],
      ...,
      [  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5],
      [  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5],
      [  0. ,   0.5,   1. , ..., 358.5, 359. , 359.5]]))

見難いかもしれないが、GPVデータは各グリッドポイントの値のnumpy配列、緯度のnumpy配列、経度のnumpy配列のタプルで構成されている。

このデータはGSMは全球データなので、地球全体のGPVデータを含んでいるが、日本周辺の天気図を作成する場合、全データは必要がない。

その場合には、data()メソッドのlat1/lat2/lon1/lon2パラメーターに緯度経度の最小値/最大値を代入して、必要なデータだけ取り出すことができる。

>>> list[0].data(lat1 = 10, lat2 = 70, lon1 = 70, lon2 = 180)
(array([[256.61254883, 256.76879883, 256.98754883, ..., 258.09692383,
       258.08129883, 258.00317383],
      [256.87817383, 257.12817383, 257.42504883, ..., 258.47192383,
       258.45629883, 258.37817383],
      [257.08129883, 257.40942383, 257.76879883, ..., 258.55004883,
       258.56567383, 258.53442383],
      ...,
      [269.50317383, 269.48754883, 269.48754883, ..., 268.08129883,
       267.86254883, 267.64379883],
      [269.50317383, 269.59692383, 269.65942383, ..., 267.94067383,
       267.80004883, 267.53442383],
      [269.40942383, 269.55004883, 269.65942383, ..., 267.65942383,
       267.40942383, 267.12817383]]), array([[70. , 70. , 70. , ..., 70. , 70. , 70. ],
      [69.5, 69.5, 69.5, ..., 69.5, 69.5, 69.5],
      [69. , 69. , 69. , ..., 69. , 69. , 69. ],
      ...,
      [11. , 11. , 11. , ..., 11. , 11. , 11. ],
      [10.5, 10.5, 10.5, ..., 10.5, 10.5, 10.5],
      [10. , 10. , 10. , ..., 10. , 10. , 10. ]]), array([[ 70. ,  70.5,  71. , ..., 179. , 179.5, 180. ],
      [ 70. ,  70.5,  71. , ..., 179. , 179.5, 180. ],
      [ 70. ,  70.5,  71. , ..., 179. , 179.5, 180. ],
      ...,
      [ 70. ,  70.5,  71. , ..., 179. , 179.5, 180. ],
      [ 70. ,  70.5,  71. , ..., 179. , 179.5, 180. ],
      [ 70. ,  70.5,  71. , ..., 179. , 179.5, 180. ]]))

このように、pygribを使うことで、GPVデータをデコードし、numpyの配列としてデータを取得できた。

GMTで天気図を描画する

ここでは例として、300hPaの気温と風速の天気図を描いてみる。AUPQ35の上段の図に相当する。

まずは、下記のようにインポートする。

>>> import pygrib as grib
>>> import pygmt as gmt
>>> import numpy as np

GPVのファイルオープンは同じだが、風速のデータは、U componentとV componentという2つのデータのベクトルで表されているので、下記のように2つのデータを読み出す。

>>> u_comp = gpv.select(shortName = 'u', level = 500)
>>> v_comp = gpv.select(shortName = 'v', level = 500)

タプルをアンパックし、風速、緯度、経度に分解する。

>>> u_wind, lat, long = u_comp[0].data(lat1 = 10, lat2 = 70, lon1 = 70, lon2 = 180)
>>> v_wind, _, _ = v_comp[0].data(lat1 = 10, lat2 = 70, lon1 = 70, lon2 = 180)

GPVの風速データの単位は[m/s]、天気図にはノット[KT]なので、変換して風速を計算する。

>>> speed = np.sqrt(np.square(u_wind) + np.square(v_wind)) * 1.9438

さらに、speedはnumpyの二次元配列なのだが、GMTへ渡すには一次元化する必要がある。緯度、経度も同様だ。

>>> speed.shape
(121, 221)
>>> speed = speed.flatten()
>>> speed.shape
(26741,)
>>> lat = lat.flatten()
>>> long = long.flatten()

これでデータ準備ができた。ようやくここからGMTの出番となる。

まず、surface()メソッドでグリッドデータに変換する。

>>> grid_wind = gmt.surface(x = long, y = lat, z = speed, region = '70/180/10/70', spacing = '0.1')

パラメーターの詳細はGMTのマニュアルなどを見て欲しい。基本的には、コマンドラインでGMTを実行する際のパラメーターと同じで問題ないようだ。

あとは地図上にデータを描くだけだ。

>>> fig = gmt.Figure()
>>> fig.grdimage(grid_wind, R = '100/15/180/55r', C = cptfile, J = 'U53/23c', B = 'f5g10')
>>> fig.grdcontour(grid_wind, J = 'U53/23c', B = 'f5g10', R = '100/15/180/55r', A = 20, W = '0.5p,black,.')
>>> fig.colorbar(D = '5c/-0.5c/8c/0.25ch', B = 'g60a60+l"Wind Speed(KT)"', C = cptfile)
>>> fig.coast(R = '100/15/180/55r', J = 'U53/23c', B = 'f5g10', W = '0.25p', D = 'f')
>>> fig.text(text = '2020/07/04_1200UTC_300hPa_T=0', x = 150, y = 18)
>>> fig.savefig('out.png')

GMTを使用したことがあれば、それぞれgrdimage、grdcontourなどと同じパラメーターが使えることが理解できるだろう（だから詳細な説明は端折る）。

パレットファイルcptfileはGMT4.xで使用していたものをそのまま指定して問題なかった。

高度場も追加すると、下記のような天気図ができあがる（表紙タイトルと同じ）。

さいごに

描画のところはさらっと簡単に説明してしまった。ただ、GMTを使った経験があれば、コードとPyGMTのマニュアルと付き合わせて理解できると思う。

PythonからGMTを使う利点は大きい。Bashプログラムを使っていると、2度中間ファイルを生成していたのだが、それが不要になった。また、オブジェクト指向言語なので、クラスを使えば再利用性もよくなり、メンテナンスにも有効になる。

全てをPythonに置き換えるまでには時間がかかりそうだが、上手くクラス定義をして再構築したい。

天気図の作図にあたっては、京都大学生存圏研究所の生存圏データベースのGPVデータを使用させて頂いています。

よろしければ、サポートをお願いします。より多くの方に役立つnoteを書けるよう頑張ります！！