GPVデータファイルからMetPyで風向風速を求める
この記事は、SDGsの目標１３「気候変動に具体的な対策を」をテーマに、気象データの利用方法の普及と気象関連ソリューション開発を促進するためのナレッジを提供することを目的として作成しています。
前回の記事では、気象庁メソモデル（MSM）から地上の風向と風速を自力で求める方法について紹介しました。
今回は、MetPyを使って気象庁メソモデル（MSM）から地上の風向と風速を求める方法について解説していきます。
目次
GPVデータファイルからMetPyで風向風速を求める
MetPyは、気象データの読み取り、視覚化、および物理計算が可能なPythonのモジュールです。
Pythonのバージョン3.9以上で利用でき、BSD 3-Clause Licenseのもとで提供されているため、商用利用を含む様々なプロジェクトで無料で利用可能です。
MetPyをpipでインストールする
Pythonの仮想環境を作成して、起動した仮想環境内にMetPyをpipでインストールします。
(MetPy) [reafnex@devsrv01 MetPy]$ pip3.12 install metpyインストール完了後、インストールされたPythonモジュールをリストして確認してみます。
(MetPy) [reafnex@devsrv01 MetPy]$ pip3.12 list
Package Version
------------------ --------
certifi 2024.2.2
charset-normalizer 3.3.2
contourpy 1.2.0
cycler 0.12.1
fonttools 4.49.0
idna 3.6
kiwisolver 1.4.5
matplotlib 3.8.3
MetPy 1.6.1
numpy 1.26.4
packaging 23.2
pandas 2.2.0
pillow 10.2.0
Pint 0.23
pip 23.2.1
platformdirs 4.2.0
pooch 1.8.0
pyarrow 15.0.0
pygrib 2.1.5
pyparsing 3.1.1
pyproj 3.6.1
python-dateutil 2.8.2
pytz 2024.1
requests 2.31.0
scipy 1.12.0
setuptools 69.1.0
six 1.16.0
traitlets 5.14.1
typing_extensions 4.9.0
tzdata 2024.1
urllib3 2.2.0
xarray 2024.1.1依存関係や関連性のあるPythonモジュールも一緒にインストールされました。PyGribやPandasもインストールされています。
MetPyで風向風速を求めるクラスを作る
今回は、前回の記事で作成した風向風速を求めるクラスを改造していきます。
MetPyを組み込んで改造したクラスの全体は以下のとおりです。
import pygrib
import metpy.calc as mpcalc
from metpy.units import units
from datetime import timedelta
class MetPhysModule :
def __init__(self, reflat, reflon, grib2) :
# 抽出する地点の緯度経度と抽出範囲を設定する
self.req_lat1 = reflat - 0.025
self.req_lat2 = reflat + 0.025
self.req_lon1 = reflon - 0.03125
self.req_lon2 = reflon + 0.03125
# GRIB2ファイル名を設定する
self.grib2dat = grib2
# UTCから日本時間に変換するための時刻差を設定する
self.diff_time = timedelta(hours=9)
def calcWindData(self) :
# GPVファイル名を指定してオープンする
gpv_file = pygrib.open(self.grib2dat)
# 地上風のコンポーネントベクトルを取り出す
u_comp = gpv_file.select(name="10 metre U wind component")
v_comp = gpv_file.select(name="10 metre V wind component")
# 予想時刻と抽出した要素データを格納するリストを初期化する
valid_dates = []
wind_dir = []
wind_spd = []
# 風向と風速をコンポーネントベクトルから計算する
for uwind, vwind in zip(u_comp, v_comp) :
# 予想時刻をリストに抽出する
valid_dates.append(uwind.validDate + self.diff_time)
u_pt = uwind.data(
lat1 = self.req_lat1, lat2 = self.req_lat2,
lon1 = self.req_lon1, lon2 = self.req_lon2
)[0][0][0]
v_pt = vwind.data(
lat1 = self.req_lat1, lat2 = self.req_lat2,
lon1 = self.req_lon1, lon2 = self.req_lon2
)[0][0][0]
# MetPyで風向風速を算出する
wdir = mpcalc.wind_direction(
u_pt * units('m/s'), v_pt * units('m/s'), convention='from'
)
wspd = mpcalc.wind_speed(
u_pt * units('m/s'), v_pt * units('m/s')
)
wind_dir.append(wdir.magnitude)
wind_spd.append(wspd.magnitude)
return valid_dates, wind_dir, wind_spdMetPyで風向風速を処理する
風向風速など物理量計算を行うためのmetpy.calcをmpcalcとしてインポートしています。また、MetPyでは、計算対象の物理量に対して単位を設定しなければならないため、metpy.unitsモジュールもインポートします。
これにより、気象データの物理量計算や単位の設定ができるようになります。
import pygrib
import metpy.calc as mpcalc
from metpy.units import units
from datetime import timedeltaMetPyで風向を計算する場合、metpy.calcモジュールのwind_directionメソッドを使用します。wind_directionメソッドの引数には、U成分ベクトル、V成分ベクトルとconvention=’from’を与えます。
wind_directionメソッドに与えるU成分ベクトルとV成分ベクトルには、単位を設定しなければなりません。単位は、「u_pt * units(‘m/s’)」のように設定します。単位が設定された変数はQuantityオブジェクトとして、物理量の値と単位のセットで扱われるようになります。
convention=’from’はデフォルトであり、これを設定することで、風が吹いてくる方向を計算して出力させることができます。他にもconvention=’to’で風が流れていく方向を計算して出力させることも可能です。
# MetPyで風向風速を算出する
wdir = mpcalc.wind_direction(
u_pt * units('m/s'), v_pt * units('m/s'), convention='from'
)
wspd = mpcalc.wind_speed(
u_pt * units('m/s'), v_pt * units('m/s')
)metpy.calcモジュールで処理された結果についても、物理量の値と単位のセットになったQuantityオブジェクトとして出力されます。物理量の値のみ出力させる場合は、Quantityオブジェクトにmagunitudeを付けます。
wind_dir.append(wdir.magnitude)
wind_spd.append(wspd.magnitude)風向風速を出力させてみる
MetPyを使って改造したクラスをインポートして、風向風速を出力してみます。Pythonのメインスクリプトは以下のとおりです。
from MetPhysModule import MetPhysModule
import pandas as pd
# 東京駅の緯度経度を設定する
land_lat = 35.681236
land_lon = 139.767125
# 気象庁MSMのGRIB2ファイル名を指定する
msm_gpv ='Z__C_RJTD_20240205000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin'
# MetPhysModuleクラスをインスタンス化する
calcwind = MetPhysModule(land_lat, land_lon, msm_gpv)
# calcWindDataメソッドを実行して、予報時刻、風向および風速のリスト変数を取得する
vt, wdir, wspd = calcwind.calcWindData()
# 予報時刻、風向および風速のリスト変数をPandasのデータフレームに入力する
df = pd.DataFrame({"ValidTime":vt, "WindDirection":wdir, "WindSpeed":wspd})
print(df)スクリプトを実行すると、以下のように出力されます。
(MetPy) [reafnex@devsrv01 MetPy]$ python3.12 main.py
ValidTime WindDirection WindSpeed
0 2024-02-05 09:00:00 16.92320143567693 0.858167
1 2024-02-05 10:00:00 322.874198396332 0.869812
2 2024-02-05 11:00:00 9.584191496389948 1.009093
3 2024-02-05 12:00:00 12.073230919555371 1.632468
4 2024-02-05 13:00:00 28.555091962085285 0.687119
5 2024-02-05 14:00:00 356.8093392845967 1.757192
6 2024-02-05 15:00:00 354.0491390949364 1.916778
7 2024-02-05 16:00:00 352.3529425373692 2.367616
8 2024-02-05 17:00:00 356.4629784162943 2.586107
9 2024-02-05 18:00:00 356.308984674244 2.771376
10 2024-02-05 19:00:00 352.34369567846204 2.968715
11 2024-02-05 20:00:00 356.07497335005314 3.361757
12 2024-02-05 21:00:00 355.7850835443618 3.645686
13 2024-02-05 22:00:00 345.4577411456487 3.819696
14 2024-02-05 23:00:00 337.5201470494043 3.444165
15 2024-02-06 00:00:00 358.9113010421894 2.701475MetPyを使わずに自力で求めた風向風速は以下のとおりです。
ValidTime WindDirection WindSpeed
0 2024-02-05 09:00:00 16.923201 0.858167
1 2024-02-05 10:00:00 322.874198 0.869812
2 2024-02-05 11:00:00 9.584191 1.009093
3 2024-02-05 12:00:00 12.073231 1.632468
4 2024-02-05 13:00:00 28.555092 0.687119
5 2024-02-05 14:00:00 356.809339 1.757192
6 2024-02-05 15:00:00 354.049139 1.916778
7 2024-02-05 16:00:00 352.352943 2.367616
8 2024-02-05 17:00:00 356.462978 2.586107
9 2024-02-05 18:00:00 356.308985 2.771376
10 2024-02-05 19:00:00 352.343696 2.968715
11 2024-02-05 20:00:00 356.074973 3.361757
12 2024-02-05 21:00:00 355.785084 3.645686
13 2024-02-05 22:00:00 345.457741 3.819696
14 2024-02-05 23:00:00 337.520147 3.444165
15 2024-02-06 00:00:00 358.911301 2.701475小数点以下の桁数の違いはありますが、同じ結果が得られました。
まとめ
いかがでしたでしょうか。
今回は、PythonのPygribモジュールを利用して気象庁MSMのGRIB2から特定地点の風成分ベクトルを抽出し、MetPyモジュールを利用して風向と風速を求めるコードを紹介しました。
MetPyモジュールは、気象学者やデータサイエンティストが気象データを効果的に解析するためにとても便利なライブラリです。
metpy.calcモジュールを利用する場合、このモジュールに渡す物理量には単位を設定する必要があります。MetPyは、物理量の単位を厳密に扱うことで計算結果の確からしさを担保しています。
今後もMetPyモジュールを利用した気象物理量の求め方について紹介していこうとおもいます。
参考になれば幸いです。
気象データ利用のお悩みについてご相談ください
気象システムの設計製造及び運用保守にかかわる現役システムエンジニアが、気象予報士の視点で気象データを活用したビジネスアイディアの創出と課題解決のお手伝いをさせていただきます。
気象や波浪などの数値予報モデルを活用した気象予報システム開発や運用についても是非ご相談ください。
