天気予報データを解析してみよう 〜数値予報から時系列グラフを作る〜

この記事では、気象庁が配信する天気予報の基礎資料となる数値予報データを扱って、将来の気圧の時系列グラフを作成する方法を紹介します。

はじめに

みなさんが見ている天気予報は、さまざまな機関（気象庁）や民間の気象会社がそれぞれ持つノウハウ等を注ぎ込んで将来の天気を予想した結果です。ただそのベースは共通していて、コンピュータが将来の地球表面のようすをシミュレーションによって予測しています。この手法は、機械で数値的に計算することから数値予報と呼ばれています。

日本では気象庁が数値予報の業務を行っていて、その予報結果は気象業務支援センターを通じて契約者向けに公開されています。公開されたファイルには、気温、湿度、気圧、風、日射量…といった値が一定間隔に置かれた緯度格子、経度格子ごとに収録されています。こういった、等間隔に区切った格子ごとに気象要素を収録したデータはGPVデータ（Grid Point Value Data）と呼ばれます。日本語にすると格子点値データなのでそのままですね。

GPVデータの形式は、世界気象機関（WMO）が定めた二進形式格子点資料気象通報式（第2版）形式、いわゆるgrib2形式がよく用いられます。中身はテキストデータではなくバイナリデータなので、メモ帳などで開いてそのまま値を見ることはできません。そこで、今回はPythonのライブラリである pygrib を使って、値を取り出してみます。

そうしてしまえば、後は好きに可視化するのみです。今回は簡単にExcelで時系列データをグラフにしてみましょう。

環境の準備

仮想環境の有効化

この記事の内容はAnacondaの仮想環境下で実行することを前提としています。仮想環境を作っていない場合、

ターミナル

conda create -n gpv_env

のように仮想環境を作成（仮想環境の名前にあたる gpv_env の部分はお好みで変更可能）して、

ターミナル

conda activate gpv_env

で有効化しておきます。

ライブラリのインストール

以降では、GPVデータを読み込む pygrib のインストールが必要です。

ターミナル

conda install pygrib

でインストールします。

GPVデータの用意

気象庁のGPVデータは、残念ながら気象庁HPなどで簡単に入手することは出来ません。入手には気象業務支援センターに契約する必要があり、毎月数千円〜数万円の負担が生じます。ただし二次配布が可能で、京都大学生存圏研究所や、東京大学生産技術研究所が有志でデータを公開してくれています。なお、京都大学生存圏研究所は教育研究機関向けのデータ提供となっています。リンクは以下の通りです。

• 気象庁データ: オリジナルデータ: 数値予報GPV （京都大学生存圏研究所）
• GPV Data Archive （東京大学生産技術研究所）

「GPVデータ」という単語はあくまでも格子点値データという意味合いしかないので、その詳細はさまざまです。例えば、地球全体の大気の流れなどをコンピュータでシミュレーションして粗い解像度で気温や気圧を計算した全球数値予報モデル（GSM）や、日本周辺に限定してGSMよりもっと細かい解像度で計算したメソ数値予報モデル（MSM）があります。どちらも、基準となる日時（基点日時）から1時間毎や3時間毎の予報を数時間〜数百時間後まで収録しています。

今回はデータの種類をMSMの地表面GPV、基点日時を2019年10月12日午前0時（世界標準時）とします。京都大学のサイトでは http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2019/10/12/ の Z__C_RJTD_20191012000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin が該当します¹。似ているファイル名のものがたくさんあるので注意してください。ファイルサイズは約66MBとやや大きめなので、回線に余裕のある環境でダウンロードしましょう。

GPVデータの処理

以下では、Pythonコードをファイルに書いてから実行することを前提とします。今回は gpv.py というファイルを作って、そのなかにコードを書くことにします。もしVSCodeの code コマンドが使える環境であれば、ターミナル上で

ターミナル

code gpv.py

とするのが早いです。各自でファイルを新規作成しても良いでしょう。

ファイルが作れたら、pygribを使ってGPVデータを読み込むコードを書いてみます。

gpv.py

import numpy as np
import pygrib

grib = pygrib.open("Z__C_RJTD_20191012000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin")

次に気象データを取得します。たとえば、

gpv.py（つづき）

prmsl_grib = grib.select(shortName="prmsl")

で取得できます。ここで、 grib には複数の気象要素、複数の予報時間のデータが全部詰まっているので、 select() で内容を絞っています。 shortName="prmsl" は気象要素に海面更正気圧（Pressure Reduced to Mean Sea Level）を指定したことになります。

ここで、一旦 prmsl_grib の中身を確認してみましょう。prmsl_grib はリストなのでfor文で以下のように1つ1つ見るコードを書いて、

gpv.py（つづき）

for x in prmsl_grib:
    print(x)

スクリプトを実行すると、

ターミナル

python gpv.py

出力結果（ターミナル上）

1:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 201910120000
11:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 1 hrs:from 201910120000
23:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 2 hrs:from 201910120000
35:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 3 hrs:from 201910120000
47:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 4 hrs:from 201910120000
59:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 5 hrs:from 201910120000
71:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 6 hrs:from 201910120000
83:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 7 hrs:from 201910120000
95:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 8 hrs:from 201910120000
107:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 9 hrs:from 201910120000
119:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 10 hrs:from 201910120000
131:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 11 hrs:from 201910120000
143:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 12 hrs:from 201910120000
155:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 13 hrs:from 201910120000
167:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 14 hrs:from 201910120000
179:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 15 hrs:from 201910120000

というのが返ってきました。リストには何やら16個のデータがあって、それぞれ fcst time n hrs の n の部分が違うようです。fcst time というのは forcast time の略で、予報時間を意味しています。更に後ろには from 201910120000 とあることから、これは2019年10月12日00時00分（協定世界時）を基準として n 時間後の気象状況を表している、という意味合いになります。従って、このリストには2019年10月12日00時ちょうど〜15時ちょうどまでの気圧データが1時間毎にあるようです。

さて、今回は特定地点の気圧を取得してみます。MSMのGPVデータは緯度が0.05度ごと、経度が0.0625度ごとに収録されているので、お試しで北緯35.70度、東経139.75度を対象にします。地図で確認するとどうやら飯田橋くらいらしいです。

先ほどはリストの1個1個のデータを確認しましたが、それに対して data() メソッドを使うことで実際の数値を取得出来ます。 data() には lat1 lat2 lon1 lon2 を指定することで、 lat1 から lat2 まで、 lon1 から lon2 までの範囲内の数値のみが返ってきます²。

gpv.py（つづき）

iidabashi_lat = 35.7
iidabashi_lon = 139.75

for x in prmsl_grib:
    # 予報時刻
    dt = x.validDate

    # 予報値
    dat = x.data(lat1=iidabashi_lat, lat2=iidabashi_lat+0.01, lon1=iidabashi_lon, lon2=iidabashi_lon+0.01)
    
    # 気圧（単位がパスカルなので0.01倍してヘクトパスカルに変換）
    prmsl = dat[0][0][0] * 0.01

    # 北緯
    lat = dat[1][0][0]

    # 東経
    lon = dat[2][0][0]
    
    print(f"{dt} {lon:.2f} {lat:.2f} {prmsl:.1f}")

気圧などを取得するときに dat[0][0][0] のように配列の配列を扱っているように見えるのは、経緯度を制限する前の x.data() が、緯度ごとの配列、経度ごとの配列という配列の入れ子で構成されていた名残です。これを実行した結果は

出力結果

2019-10-12 00:00:00 139.75 35.70 1001.9
2019-10-12 01:00:00 139.75 35.70 999.9
2019-10-12 02:00:00 139.75 35.70 999.9
2019-10-12 03:00:00 139.75 35.70 997.6
2019-10-12 04:00:00 139.75 35.70 994.5
2019-10-12 05:00:00 139.75 35.70 991.6
2019-10-12 06:00:00 139.75 35.70 988.8
2019-10-12 07:00:00 139.75 35.70 986.2
2019-10-12 08:00:00 139.75 35.70 983.6
2019-10-12 09:00:00 139.75 35.70 981.7
2019-10-12 10:00:00 139.75 35.70 979.5
2019-10-12 11:00:00 139.75 35.70 977.3
2019-10-12 12:00:00 139.75 35.70 975.0
2019-10-12 13:00:00 139.75 35.70 972.5
2019-10-12 14:00:00 139.75 35.70 970.0
2019-10-12 15:00:00 139.75 35.70 976.9

となって、確かに北緯35.7度、東経139.75度の気圧値（単位はヘクトパスカル）を取得できたようです。

先ほどのコードを

gpv.py

iidabashi_lat = 35.7
iidabashi_lon = 139.75

with open("prmsl.txt", mode="w") as f:
    # 見出し
    f.write(",東経,北緯,海面更正気圧[hPa]\n")

    for x in prmsl_grib:
        # 予報時刻
        dt = x.validDate

        # 予報値
        dat = x.data(lat1=iidabashi_lat, lat2=iidabashi_lat+0.01, lon1=iidabashi_lon, lon2=iidabashi_lon+0.01)

        # 気圧（単位がパスカルなので0.01倍してヘクトパスカルに変換）
        prmsl = dat[0][0][0] * 0.01

        # 北緯
        lat = dat[1][0][0]

        # 東経
        lon = dat[2][0][0]

        f.write(f"{dt},{lon:.2f},{lat:.2f},{prmsl:.1f}\n")

と置き換えると、 prmsl.txt に結果がカンマ区切りで出力されます。良かったですね。

出力されたファイルは、以下のようになっているはずです。

prmsl.txt

,東経,北緯,海面更正気圧[hPa]
2019-10-12 00:00:00,139.75,35.70,1001.9
2019-10-12 01:00:00,139.75,35.70,999.9
2019-10-12 02:00:00,139.75,35.70,999.9
2019-10-12 03:00:00,139.75,35.70,997.6
2019-10-12 04:00:00,139.75,35.70,994.5
2019-10-12 05:00:00,139.75,35.70,991.6
2019-10-12 06:00:00,139.75,35.70,988.8
2019-10-12 07:00:00,139.75,35.70,986.2
2019-10-12 08:00:00,139.75,35.70,983.6
2019-10-12 09:00:00,139.75,35.70,981.7
2019-10-12 10:00:00,139.75,35.70,979.5
2019-10-12 11:00:00,139.75,35.70,977.3
2019-10-12 12:00:00,139.75,35.70,975.0
2019-10-12 13:00:00,139.75,35.70,972.5
2019-10-12 14:00:00,139.75,35.70,970.0
2019-10-12 15:00:00,139.75,35.70,976.9

可視化

最後に、Excelで可視化してみましょう。先ほどの prmsl.txt をコピペして貼り付けると、

となります。「データ」タブにある「区切り位置」のメニューをクリックして、データがカンマ区切りであることを教えると、

ちゃんと区切られました。

線付きの散布図で良い感じにグラフ化すると、以下のようなものが描けました。この日のお昼過ぎにかけて970ヘクトパスカル程度まで気圧が下がったあと、上昇し始めるくらいの予想までがグラフに載っています。

ちなみにこの2019年10月12日というのは2019年台風19号（令和元年東日本台風）が日本に上陸した日で、その影響で気圧がぐっと下がっているわけですね。

おわりに

この記事では、Pythonライブラリ pygrib を用いて数値予報GPVデータを解析したあと、Excelで可視化を行うところまでを解説しました。
今回使用したGPVファイルは1つなので、解析したのはある時刻時点における予報データです。ただ、後から当時を振り返るような用途では、予報ではなく実況（あるいは推定）データが必要になります。このような場合は複数のGPVファイルを使用する必要がありますが、その方法についてはまた記事を書くかもしれません。

1. Lsurfは地表面を、FH00-15は予報時間（Forcast Hour）が0時間後〜15時間後までであることを意味します。 ↩

2. おそらくlat1"以上"、lat2"以下"のように以上以下の組合せだとは思うのですが、以上なのか超なのか、以下なのか未満なのかの正確な確認はしていません。 ↩