どこにでもいる30代SEの学習ブログ

主にプログラミング関連の学習内容。読んだ本の感想や株式投資についても書いてます。

トップ > プログラミング > pygribで気象庁の数値予報GPVデータを読み込む

pygribで気象庁の数値予報GPVデータを読み込む

プログラミング プログラミング-Python プログラミング-Python-pygrib 気象 気象-数値予報 気象-pygrib 気象-数値予報-GPV

f:id:predora005:20200921202207j:plain ※ 2020/08/16にQrunchで書いた記事を移行しました。

気象庁の数値予報データをPythonで扱いました。 pygribというライブラリでデータを読み込み、数値予報データの中身を確認する方法を紹介します。

[1] 数値予報GPVデータ

[1-1] 数値予報とは

気象庁が行っている気象予測の一つが数値予報です。地球全体を一定の大きさの格子に区切り、各格子点の気温や風速などを予測しています。GPVはGrid Point Valueの略で、各格子点の値を指します。 詳しくは気象庁のページを見ていただいた方がわかるでしょう。

[1-2] 数値予報の種類

全球数値予報モデル(GSM), メソ数値予報モデル(MSM), 局地数値予報モデル(LFM)の3種類があります。 今回扱ったのは全球数値予報モデル(GSM)です。

名称 説明 格子点の間隔(水平方向)
全球数値予報モデル(GSM) 地球全体の大気を対象とした気象庁の数値予報モデル。Global Spectral Model。 日本域は20km四方(緯度0.2度x経度0.25度)。全球域は緯度0.5度x経度0.5度
メソ数値予報モデル(MSM) 日本及びその近海の大気を対象とした気象庁の数値予報モデル。Meso-Scale Model。 5km四方
局地数値予報モデル (LFM） 日本領域の大気を対象とした気象庁の数値予報モデル。Local Forecast Model。 2km四方

[1-3] データの形式

GRIB2という気象データ向けのバイナリフォーマットです。

[1-4] データの取得元

数値予報データは基本的には有償で、月額¥15,000ほどかかります。 しかしながら、ありがたいことに京都大学生存研究所が教育機関向けにデータを提供してくださっているので、こちらを使わせていただきました。

[2] pygrib

pygribは、GRID形式のデータを読み書きするためのPythonモジュールです。 私はAmazon Linux2を使用していますが、Amazon Linux2でのインストール方法は「pygribのインストール方法：Amazon Linux2」にまとめてあります。

[3] 全球数値予報モデル(GSM)データの読み込み

[3-1] GRIB2形式のファイルをオープン

GSMの日本域データは、地上データと気圧面データ(高層のデータ)とに分かれています。 まずは、地上データをダウンロードし、開きます。

# coding: utf-8

import os
import subprocess
import pygrib

# 2020/08/01 15:00(UTC 06:00)の地上データをダウンロードする
cwd = os.getcwd()
url_surf = 'http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2020/08/01/Z__C_RJTD_20200801060000_GSM_GPV_Rjp_Lsurf_FD0000-0312_grib2.bin'
subprocess.run(['curl', '-O', url_surf], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL, cwd=cwd)

## GRIB2ファイルを読み込む
file_surf = os.path.join(cwd, os.path.basename(url_surf))
grbs = pygrib.open(file_surf)

[3-2] 中身を表示する

# GRIB2ファイルの中身を表示する
for grb in grbs:
    print(grb)

for文でデータを取り出しprint文で表示すると、下記の結果が得られます。

1:Pressure reduced to MSL:Pa (instant):regular_ll:meanSea:level 0:fcst time 0 hrs:from 202008010600
2:Surface pressure:Pa (instant):regular_ll:surface:level 0:fcst time 0 hrs:from 202008010600
3:10 metre U wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10 m:fcst time 0 hrs:from 202008010600
4:10 metre V wind component:m s**-1 (instant):regular_ll:heightAboveGround:level 10 m:fcst time 0 hrs:from 202008010600

...(中略)...

1015:High cloud cover:% (instant):regular_ll:surface:level 0:fcst time 84 hrs:from 202008010600
1016:Total cloud cover:% (instant):regular_ll:surface:level 0:fcst time 84 hrs:from 202008010600
1017:Total precipitation:kg m-2 (accum):regular_ll:surface:level 0:fcst time 0-84 hrs (accum):from 202008010600
1018:Downward short-wave radiation flux:W m**-2 (avg):regular_ll:surface:level 0:fcst time 83-84 hrs (avg):from 202008010600

「level」は高さ(主に気圧面データで使う)を示し「fcst time」は初期時刻から何時間後のデータであるかを示しています。

GSM日本域の地上データは初期時刻(この場合は15時)から、1時間ごとに84時間後まで予測データを含んでいます。気圧面データは3時間ごとです。

[3-3] 指定した観測データを取り出す(地上データ)

selectメソッドで取り出します。parameterNameに観測データの名称を指定します。

# 指定した観測データを取り出す
prmsl = grbs.select(parameterName='Pressure reduced to MSL')
sp    = grbs.select(parameterName='Pressure')
uwind = grbs.select(parameterName='u-component of wind')
vwind = grbs.select(parameterName='v-component of wind')
temp  = grbs.select(parameterName='Temperature')
rh    = grbs.select(parameterName='Relative humidity')
lcc   = grbs.select(parameterName='Low cloud cover')
mcc   = grbs.select(parameterName='Medium cloud cover')
hcc   = grbs.select(parameterName='High cloud cover')
tcc   = grbs.select(parameterName='Total cloud cover')
tp    = grbs.select(parameterName='Total precipitation')
dswrf = grbs.select(parameterName='Downward short-wave radiation flux')

戻り値はリスト形式で、初期時刻から84時間後のデータまでが含まれます。 nameやshortNameでも指定可能ですが、'Total cloud cover'と'Total precipitation'のname, shortNameがunknownのため、parameterNameで指定しています。

[3-4] 指定した予報時刻のデータを取り出す(地上データ)

selectメソッドで取り出します。初期時刻から何時間後のデータを取り出すのかをforecastTimeに整数値で指定します。

# 初期時刻から12時間後のデータを取り出す
fc12 = grbs.select(forecastTime=12)

[3-5] 観測値を取り出す

指定した観測データの数値を取り出したい場合は、dataメソッドを使います。 戻り値として、観測値, 緯度, 経度が返ってきます。

# 観測値を取得する
prmsl_fc0 = grbs.select(parameterName='Pressure reduced to MSL', forecastTime=0)[0]
values, lats, lons = prmsl_fc0.data()
print(values)
# [[ 99642.71240234  99664.58740234  99680.21240234 ... 101829.43115234
#   101820.05615234 101809.89990234]
#  [ 99645.05615234  99666.93115234  99685.68115234 ... 101840.36865234
#   101830.99365234 101820.05615234]
#  [ 99660.68115234  99677.08740234  99695.83740234 ... 101852.08740234
#   101842.71240234 101830.99365234]
#  ...(以下省略)...

取り出し元のデータは、特定時刻の特定の観測データである必要があります。

[3-6] 指定した範囲(緯度, 経度)のデータを取り出す

dataメソッドに緯度, 経度の範囲を指定します。 緯度, 経度を指定しない場合は、全領域のデータが返されます。

# 北緯34度から36度、東経135度から140度のデータを取り出す
prmsl_fc0 = grbs.select(parameterName='Pressure reduced to MSL', forecastTime=0)[0]
values, lats, lons = prmsl_fc0.data(lat1=34, lat2=36, lon1=135, lon2=140)
print(values)
#[[101209.11865234 101213.80615234 101219.27490234 101216.93115234
#  101197.39990234 101138.02490234 101043.49365234 100913.80615234
#  100827.86865234 100809.11865234 100765.36865234 100755.21240234
#  100780.21240234 100797.39990234 100867.71240234 100953.64990234
#  101014.58740234 101058.33740234 101093.49365234 101123.96240234
#  101152.08740234]
#  ...(以下省略)...
print(lats)
# [[36.  36.  36.  36.  36.  36.  36.  36.  36.  36.  36.  36.  36.  36.
#   36.  36.  36.  36.  36.  36.  36. ]
#  [35.8 35.8 35.8 35.8 35.8 35.8 35.8 35.8 35.8 35.8 35.8 35.8 35.8 35.8
#   35.8 35.8 35.8 35.8 35.8 35.8 35.8]
#  ...(以下省略)...
print(lons)
# [[135.   135.25 135.5  135.75 136.   136.25 136.5  136.75 137.   137.25
#   137.5  137.75 138.   138.25 138.5  138.75 139.   139.25 139.5  139.75
#   140.  ]
#  ...(以下省略)...

全データ(grbs)から取り出せればよかったのですが、それは出来ないようです。

[3-7] 気圧面データを開き、中身を表示する

今度は、気圧面データを開き、中身を表示してみます。

# 2020/08/01 15:00(UTC 06:00)の気圧面データをダウンロードする
cwd = os.getcwd()
url_pall = 'http://database.rish.kyoto-u.ac.jp/arch/jmadata/data/gpv/original/2020/08/01/Z__C_RJTD_20200801060000_GSM_GPV_Rjp_L-pall_FD0000-0312_grib2.bin'
subprocess.run(['curl', '-O', url_pall], stdout=subprocess.DEVNULL, stderr=subprocess.DEVNULL, cwd=cwd)

## GRIB2ファイルを読み込む
file_pall = os.path.join(cwd, os.path.basename(url_pall))
grbs = pygrib.open(file_pall)

# GRIB2ファイルの中身を表示する
for grb in grbs:
    print(grb)

for文でデータを取り出しprint文で表示すると、下記の結果が得られます。

1:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 0 hrs:from 202008010600
2:U component of wind:m s**-1 (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 0 hrs:from 202008010600
3:V component of wind:m s**-1 (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 0 hrs:from 202008010600
4:Temperature:K (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 0 hrs:from 202008010600
5:Vertical velocity:Pa s**-1 (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 0 hrs:from 202008010600
6:Relative humidity:% (instant):regular_ll:isobaricInhPa:level 100000.0 Pa:fcst time 0 hrs:from 202008010600
...(以下省略)...

気圧面で扱えるのは、上記の6観測データです。 parameterNameを用いた取り出し方法は次の通りです。

# 指定した観測データを取り出す
gh    = grbs.select(parameterName='Geopotential height')
uwind = grbs.select(parameterName='u-component of wind')
vwind = grbs.select(parameterName='v-component of wind')
temp  = grbs.select(parameterName='Temperature')
vv    = grbs.select(parameterName='Vertical velocity (pressure)')
rh    = grbs.select(parameterName='Relative humidity')

[3-8] 指定気圧面のデータを取り出す(気圧面データ)

気圧面データは、指定気圧面と呼ばれる同じ気圧の平面上に格子点を持っています。1000hPa面, 925hPa面など。 selectメソッドでlevelに気圧面を指定します。

# 指定した気圧面のデータを取り出す
data_850hPa = grbs.select(level=850)
data_700hPa = grbs.select(level=700)
data_500hPa = grbs.select(level=500)

[3-9] 取り出したデータから緯度, 経度を取得する

dataメソッドでも取得可能ですが緯度, 経度のみを取得したい場合はlatlonsメソッドを使います。

# 緯度,経度を取得する。
rh_500hPa = grbs.select(parameterName='Relative humidity', level=500, forecastTime=0)[0]
lats, lons = rh_500hPa.latlons()
print(lats)
# [50.  50.  50.  ... 50.  50.  50. ]
#  [49.8 49.8 49.8 ... 49.8 49.8 49.8]
#  [49.6 49.6 49.6 ... 49.6 49.6 49.6]
#  ...
#  [20.4 20.4 20.4 ... 20.4 20.4 20.4]
#  [20.2 20.2 20.2 ... 20.2 20.2 20.2]
#  [20.  20.  20.  ... 20.  20.  20. ]]
print(lons)
# [[120.   120.25 120.5  ... 149.5  149.75 150.  ]
#  [120.   120.25 120.5  ... 149.5  149.75 150.  ]
#  [120.   120.25 120.5  ... 149.5  149.75 150.  ]
#  ...
#  [120.   120.25 120.5  ... 149.5  149.75 150.  ]
#  [120.   120.25 120.5  ... 149.5  149.75 150.  ]
#  [120.   120.25 120.5  ... 149.5  149.75 150.  ]]

[3-9] GRIB2ファイルを閉じる

closeメソッドで閉じます。

# GRIB2ファイルを閉じる
grbs.close()

終わりに

最初の方は取っ付きづらく感じますが、慣れてしまえば難しいことはそれほどありませんでした。

注意した方が良いのは、ファイルサイズが大きいことです。 GSM日本域データは1ファイルのサイズが地上データで25MB, 気圧面データで70MBと非常に大きいです。6時間ごとに予測が更新されるため、1日に4ファイルあります。数年分のデータなどを扱う場合には、一定以上のストレージ容量が必要です。

また、データの読み出しに結構時間がかかる点も注意が必要です。

参考文献