GPVデータファイルから自力で風向風速を求める
この記事は、SDGsの目標１３「気候変動に具体的な対策を」をテーマに、気象データの利用方法の普及と気象関連ソリューション開発を促進するためのナレッジを提供することを目的として作成しています。
今回は、Pythonを使って気象庁メソモデル（MSM）から地上の風向と風速を自力で求める方法について解説していきます。
数値予報モデルやGPVデータファイル、Pythonを使ったGPVデータ処理の基本的方法については、以下の記事で紹介しています。
GPVデータファイルから自力で風向風速を求める
気象庁のメソモデル（MSM)、全球モデル（GSM）や局地モデル（LFM）などの数値予報結果は、GRIB2というデータフォーマットで出力されています。
GRIB2には、各種数値予報モデルが計算した風、気温、湿度、気圧などの物理パラメーターが２次元配列で記録されており、１回分の数値予報結果は、地上、上層、予報時間などの単位の複数のGRIB2に分割されて提供されています。
また、風は一般的に風向と風速で表現されますが、GRIB2に記録されている風の物理パラメータは、東西成分と南北成分の２つのベクトルとして表現されています。そのため、２つの成分ベクトルを合成して風向と風速を求める必要があります。
風の成分ベクトルとは
数値予報モデルの中では、風は南北成分（V成分）と東西成分（U成分）の２つのベクトルで処理されています。
風の南北成分（V成分）と東西成分（U成分）をベクトル合成することで、「風の流れていく方向」と「風速」を求めることができます。また、ベクトル合成で求めた「風の流れていく方向」から「風の吹いてくる方向」に角度を変換することで、風向を求めることができます。
Pythonで風向風速を求めるクラスを作る
ここからは、ある地点の風向風速を気象庁MSMのGRIB2から求めるPythonのコードについて紹介していきます。
なお、Pythonはオブジェクト指向のコーディングが可能であるため、風向風速を求めるクラスを作成して、再利用できるようにコーディングしていきます。
クラスコードの全体は以下のとおりです。
import pygrib
from datetime import timedelta
import math
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]
if v_pt == 0 and u_pt == 0 :
wdir = 0.0
wspd = 0.0
else :
# 風向計算
angle_rad = math.atan2(v_pt, u_pt)
angle_deg = math.degrees(angle_rad)
wdir = ((270.0 - angle_deg) + 360.0) % 360.0
# 風速計算
wspd = math.sqrt(u_pt ** 2.0 + v_pt ** 2.0)
wind_dir.append(wdir)
wind_spd.append(wspd)
return valid_dates, wind_dir, wind_spdこのクラスはPythonのPygribモジュールをインポートしてGRIB2データを読み込みます。予報時刻はUTCで記録されているため、これをJSTに変換するtimedeltaモジュールと、ベクトル合成するためmathモジュールもインポートします。
クラスの_コンストラクタでは、風向風速を求める地点の緯度・経度と気象庁MSMのGRIB2ファイル名を受け取ってコンストラクタ変数を初期化しています。
import pygrib
from datetime import timedelta
import math
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)calcWindDataメソッドでは、PygribモジュールのopenメソッドにGRIB2ファイル名を指定してファイルをオープンしています。オープンしたGRIB2ファイルから、地上付近の風の南北成分（V成分）と東西成分（U成分）のGRIBメッセージをv_compおよびu_compというオブジェクト変数にそれぞれ抽出し格納しています。
また、予想時刻、風向、および風速を格納するためのリスト変数（valid_dates、wind_dir、wind_spd）を初期化しています。
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 = []以下の部分では、u_compのオブジェクト変数から予報時刻をvalid_datesのリスト変数に追加しています。
また、u_compおよびv_compのオブジェクト変数から特定の緯度経度に該当する東西成分ベクトルと南北成分ベクトルをu_ptおよびv_ptという変数にそれぞれ代入しています。
# 風向と風速をコンポーネントベクトルから計算する
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]以下の部分では、東西成分ベクトルと南北成分ベクトルを格納した変数に代入した値を使って、風向風速を計算しています。
東西成分ベクトルと南北成分ベクトルの値が「0」の場合、風向（wdir）と風速（wspd）に「0.0」を代入します。それ以外の場合、風向はatan2を使用して合成ベクトルのラジアンを求め、その後ラジアンを角度に変換します。
求められた角度（angle_deg）は「風の流れていく方向」であるため、「風の吹いてくる方向」に角度を変換して風向（wdir）変数に代入します。
また、風速は東西成分ベクトルと南北成分ベクトルのそれぞれの二乗和の平方根を求めて、風速（wspd）変数に代入します。算出した風向と風速は、wind_dirとwind_spdのリスト変数に追加します。
最後に、予報時刻、風向および風速を格納したリスト変数を戻り値としてリターンしています。
if v_pt == 0 and u_pt == 0 :
wdir = 0.0
wspd = 0.0
else :
# 風向計算
angle_rad = math.atan2(v_pt, u_pt)
angle_deg = math.degrees(angle_rad)
wdir = ((270.0 - angle_deg) + 360.0) % 360.0
# 風速計算
wspd = math.sqrt(u_pt ** 2.0 + v_pt ** 2.0)
wind_dir.append(wdir)
wind_spd.append(wspd)
return valid_dates, wind_dir, wind_spd風向風速を出力させてみる
作成したクラスをインポートして、風向風速を出力してみます。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)スクリプトを実行すると、以下のように出力されます。
(Pygrib) [reafnex@devsrv01 Pygrib]$ python3 main.py
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.7014752024年2月5日は、東京タワーが立つ関東南部でも積雪となった日です。気象庁MSMの数値予報結果によると、この日は2～3m/s前後の北寄りの風が予想されていたことになります。
まとめ
いかがでしたでしょうか。
今回は、PythonのPygribモジュールを使い、気象庁MSMのGRIB2から特定地点の風成分ベクトルを抽出し、風向と風速を自力で求めるコードを紹介しました。
実は、Pythonは気象物理計算に特化したMetPyモジュールが利用可能であり、もっと簡単に風成分ベクトルから風向や風速を求めることができます。
次回は、MetPyモジュールを利用した風向と風速の求め方について紹介しようとおもいます。
参考になれば幸いです。
気象データ利用のお悩みについてご相談ください
気象システムの設計製造及び運用保守にかかわる現役システムエンジニアが、気象予報士の視点で気象データを活用したビジネスアイディアの創出と課題解決のお手伝いをさせていただきます。
気象や波浪などの数値予報モデルを活用した気象予報システム開発や運用についても是非ご相談ください。
