PythonからGMT6.xを使って天気図を描く
はじめに
前回、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を書けるよう頑張ります!!