気象庁GPVからQGISで簡単にアニメーション付きコンター図を作るプラグイン作ったよ

はじめに

こんなプラグインを作ったよって話

プラグインファイルはgithubから
アニメーションはQGIS version 3.14からしか使えないので、プラグインを使う場合は3.14以上を使ってください

GPVとは

Grid Point Value。気象観測の結果を使って数値予報モデルという予報プログラムを回して得られた、数時間先の格子点ごとの大気の物理量の値のデータが格納されています。このGPVを元に天気予報が行われる。

GPVの種類

予報する天気予報の種類によってGPVの予報モデルは数種類あり、主に以下のものがある。

モデルの種類	予報領域	格子間隔	主な予報
全休モデル(GSM)	地球全体	約２０km(0.2度 x ０.25度)	府県天気予報、 週間天気予報など
メソモデル(MSM)	日本周辺	約５km(0.1度 x ０.125度)	降水短時間予報、 防災気象情報など
局地モデル（LFM）	日本周辺	約２km（0.02度×0.025度）	降水短時間予報、 防災気象情報など

これらのGPVは気象庁から提供されており、業務等で利用する場合は気象業務支援センターを通して購入する必要がありますが、今回はサンプルデータ（MSMモデルの地上面）を使ってプラグインの検証していきました。
気象庁　GPVサンプルデータの一覧

データの命名規則

MSMの場合、「Z__C_RJTD_yyyyMMddhhmmss_MSM_GPV_Rjp _Lsurf_FHhh-hh_grib2.bin」となっていて、MSMがメソモデル、Lsurfが地表面、FHが予報期間を示しています。LsurfではなくL-pallとなってるものは気圧面のデータであることを示してます。
Lsurfの場合1時間ごと、L-pallは3時間ごとの予測値（15時間先まで）も格納されています。

GPVをQGISで開く

GPVデータはGrib２フォーマットで提供されており、wgrib2などを使ってコマンドを叩きながらデータを取り出すのですが、複雑なことせずにQGISで開けられるって知ってますか？
そう、何を隠そうドラッグアンドドロップするだけです。

このGPVデータ。QGISで追加するとマルチバンドのラスタデータとして扱われ、追加した最初の状態ではレンダリングがマルチバンド設定になってます。各バンドには高度ごとの気圧や風、気温などの値が格納されています。
上記のGPVサンプルデータの一覧ページに、どんな要素の値が格納されているかわかります。今回はMSMの地上面のGPVを使っているので、海面更正気圧や風ベクトル、気温などの要素がつまってるわけですが、バンドを見ると何が何やらわかりません。

そこで、プロパティの「情報」を見ると、各バンドの情報を見ることができる。
例えばバンド１ならPressure reduced to MSL（海面更正気圧）であることがわかる。
また、同じく情報にFORECAST_SECONDSというのがあるのがわかる通り、予測値も別バンドに格納されている。
海面更正気圧の場合、バンド1（初期値）、11（1時間先）、23（2時間先）、35（3時間先）、47(4時間先)・・・、と14時間先までの予報値も入っているという作りになっています。

そこで、QGISで特定の要素について初期値と予測値それぞれの等値線を作って、ついでにQGIS ver.3.14から搭載されたTemporal Controllerを使ってアニメーションを簡単に作れるプラグインを作ってみました。

プラグインについて

使い方

• プラグインをgithubからダウンロードして、QGISの「プラグインの管理とインストール」からzipを指定してインストールします
  

• インストールするとアイコンが追加されます。いかにもコンターっぽいアイコンです

• アイコンをクリックするとプラグインが起動します。
  

• 対象とするGPVデータとコンターを作成する要素を選択します。
  現在のところ、地上面の場合は海面更正気圧と地上気圧、気圧面の場合は各気圧面高度を指定できます。
  

• 実行するとコンターが作られます。
  予測時間ごとのコンターが全て出力されるので、何やらごちゃごちゃしてみにくいかもしれません。
  

• QGIS ver.3.14から搭載されたTemporal controllerを起動して、アニメーションを設定する
  Stepを、Lsurfの場合は30minutes、L-pallの場合1.5hourにしてるといい感じのアニメーションになると思います。
  

• 再生ボタンを押すとアニメーションが再生されます。
  

出力コンターのレイヤ名

出力されるコンターのレイヤ名は以下のようになっています。

PyQGISでやっていること

コンターの一括作成

コンターを作るのはQGISのツールバーにある「ラスタ」→「抽出」→「等高線（コンター）」で可能ですが、pyQGISで書くとこのようになります。

processing.run("gdal:contour", {'INPUT': 入力ラスタレイヤ, 'BAND': コンター化するバンド, 'INTERVAL': interval, 'FIELD_NAME': 'Contour','CREATE_3D': False, 'IGNORE_NODATA': False, 'NODATA': None, 'OFFSET': 0, 'EXTRA': '', 'OUTPUT': 'TEMPORARY_OUTPUT'})

BANDで指定するバンドを、コンター化する要素の初期値と予測値のバンドでループしてプロセッシングを実行します。
例えば、メソモデルの地上面の場合、
海面更正気圧のバンド1,11,23,35,47,59,71,83,95,107,119,131,143,155,167,179
地上気圧のバンド1,11,23,35,47,59,71,83,95,107,119,131,143,155,167,179
気圧面の場合
1000hPa高度のバンド1,93,185,277,369,461を指定するようにしています。

アニメーションの設定

出力コンターにはTemporal Controllerの設定がなされています。

vlayer = レイヤを指定
mode = QgsVectorLayerTemporalProperties.ModeFixedTemporalRange
tprops = vlayer.temporalProperties()
tprops.setMode(mode)
tprops.setFixedTemporalRange(QgsDateTimeRange(
　　　　QDateTime(QDate(YYYY, MM, DD), QTime(start hh, 0, 0, 0), Qt.UTC), 
　　　　QDateTime(QDate(YYYY, MM, DD), QTime(end hh, 0, 0, 0), Qt.UTC)))
tprops.setIsActive(True)

modeはModeFixedTemporalRange（ConfigurationでいうFixed Time Range）の他、ModeFeatureDateTimeInstantFromField（Single Field with Date/Time）、ModeFeatureDateTimeStartAndEndFromFields（Separate Fields for Start and End Date/Time）など指定できます。

ラベル

等圧線には20hPaごと、気圧面の等高線には60mごとのラベルが設定されています。

# set labeling
pal_layer = QgsPalLayerSettings()
pal_layer.isExpression = True
pal_layer.fieldName = labelSetting[0]
pal_layer.enabled = True
pal_layer.placement = QgsPalLayerSettings.Line
# Create and append a new rule
root = QgsRuleBasedLabeling.Rule(QgsPalLayerSettings())
rule = QgsRuleBasedLabeling.Rule(pal_layer)
rule.setDescription('ラベルの説明')
rule.setFilterExpression(ラベルフィールドを指定)
root.appendChild(rule)
# Apply label configuration
vlayer.setLabelsEnabled(True)
rules = QgsRuleBasedLabeling(root)
vlayer.setLabeling(rules)
vlayer.triggerRepaint()

シンボル設定

コンター線を薄い灰色に設定しています。
現状は全てのコンター線が同じ太さになってますが、今後は20hPaごとのコンター線を太くしたい。

# set symbol
single_symbol_renderer = vlayer.renderer()
symbol = single_symbol_renderer.symbol()
# Set stroke colour
symbol.setColor(QColor.fromRgb(100, 100, 100))
# Refresh
vlayer.triggerRepaint()
iface.layerTreeView().refreshLayerSymbology(vlayer.id())

今後に向けて

GPVはきちんと業務で使う場合は気象業務支援センターから購入しないといけませんが、研究目的なら以下などから利用できるようです。

GPVには気圧のほか、風ベクトルや気温や雲量やいろいろな値が入ってるので、それらの値の可視化も実装していきたいと思います。