GPVデータファイルからWBGTを求める
この記事は、SDGsの目標１３「気候変動に具体的な対策を」をテーマに、気象データの利用方法の普及と気象関連ソリューション開発を促進するためのナレッジを提供することを目的として作成しています。
今回は、PythonのMetPyを使って気象庁メソモデル（MSM）から湿球黒球温度（WBGT）を求める方法について解説していきます。
GPVデータから物理量を抽出する方法や風速、露点温度の求め方については、以下で詳しく解説していますので参考にしてください。
GPVデータファイルからWBGTを求める
今回の解説では、PythonのモジュールのPygribとMetPyを利用します。また、作成したコードは再利用できるようにすることを前提としクラスで作成していきます。
なぜWBGTを求めるのか
WBGT（Wet Bulb Globe Temperature）は、湿球黒球温度という暑さ指数であり、熱中症を予防することを目的として1954年にアメリカ海兵隊で訓練時の熱中症リスクを事前に判断するための指標提案された指標です。
WBGTは、気温、湿度と輻射熱が人体に与える影響が考慮されており、湿球温度（Tw）、黒球温度（Tg）および乾球温度（Td）を使って算出されます。
湿球温度（wet-bulb temperature）は、空気中の湿度を考慮に入れた温度の指標で、湿った布に包まれた温度計を通風環境に暴露して測定することで得られます。湿度が低いほど蒸発が盛んとなり、湿球温度は低くなります。
黒球温度（black globe temperature）は、黒色に塗装された直径約15cm程度の薄い銅板の球の中心に温度計を入れて観測される温度です。直射日光と通風環境に暴露した状態で球体の中の温度を観測しており、日なたにおける体感温度と相関関係があるとされています。
乾球温度（NDB:Natural Dry Bulb temperature）は、乾球温度計を使用して測定される温度のことであり、乾燥した状態の空気の温度を示します。一般に気温として扱われる温度のことです。
近年、温暖化により地球全体の平均温度が上昇していると言われており、SDGsの目標１３「気候変動に具体的な対策を」で温暖化の原因となる二酸化炭素などの温室効果ガスの削減が叫ばれています。
日本における近年の気候も、春と秋が短く夏が長くなったように感じられます。真夏日や酷暑日の出現日数が増加しており、熱中症リスクが増大しているため、WBGTを活用した熱中症対策の呼びかけがとても重要です。
WBGTは、通常は専用の測器を用いた観測によって得られるものですが、数日後ぐらいまでのWBGT予報についても需要があると考えられるため、数値予報結果を利用したWBGTの算出についても研究がなされています。
GPVからWBGTを求める方法
環境省によると、WBGTは以下の計算式で計算されます。
WBGT（屋外） ＝0.7 × 湿球温度（Tw） ＋ 0.2 × 黒球温度（Tg） ＋ 0.1 × 乾球温度（Ta）
ここで、GPVから直接得られる物理量は、乾球温度（Ta）だけであり、湿球温度（Tw）と黒球温度（Tg）はいくつかの物理量を組み合わせて算出する必要があります。
湿球温度（Tw）は、地上気圧、気温と露点温度を使ってMetpyで算出することができます。
一方、黒球温度（Tg）はMetpyでは求めることはできないため、以下の論文に示される黒球温度推定式で算出します。
Maki.Okada Masumi.Okada and Hiroyuki.Kusaka, Parameter adjustment and application to an extension area of Okada and Kusaka’s formula for the black globe temperature
ここで、Tgは黒球温度（℃）、Taは気温（℃）、So は全天日射量（W/m2）、U は風速（m/s）を表す．完全な物理式であれば、a〜d は気象要素から成る関数で表現されるはずである．但し、Okada and Kusaka(23)の式は、気象官署から得られない気象要素を全て経験定数 a〜d として押し込めることで黒球温度を推定する．
https://www.heat-island.jp/web_journal/download/13A003.pdf
黒球温度推定式の原型は以下のとおりで、定数a～dが不定のため、このままでは黒球温度を算出できません。
そこで彼らは、環境省の黒球温度観測データと気象官署の気温、全天日射量、風速の観測データを用いて、非線形関数の最小値を探索する滑降シンプレックス法により定数 a～d の値を推定し、以下の式を導出しています。
黒球温度推定式は、日射量（W/m2）、地上付近の風速（m/s）と地上付近の気温（℃）を使って求められることが理解できます。
WBGTを求めるクラスを作る
ここからは、ある地点のWBGTを気象庁MSMのGRIB2から求めるPythonのコードについて紹介していきます。
はじめにクラスコードの全体を以下に示します。
import pygrib
import metpy.calc as mpcalc
from metpy.units import units
from datetime import timedelta
import numpy as np
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 extractValidTime(self, physquantity) :
validtime = physquantity.validDate + self.diff_time
return validtime
def extractPointData(self, physquantity) :
pointdata = physquantity.data(
lat1 = self.req_lat1, lat2 = self.req_lat2,
lon1 = self.req_lon1, lon2 = self.req_lon2
)[0][0][0]
return pointdata
def calcWBGT(self) :
# GPVファイル名を指定してオープンする
gpv_file = pygrib.open(self.grib2dat)
# 日射量、気温、相対湿度、風速ベクトルを取り出す。
shortwave = gpv_file.select(name="Mean surface downward short-wave radiation flux")
temp = gpv_file.select(name="Temperature")
rhum = gpv_file.select(name="Relative humidity")
pres = gpv_file.select(name="Surface pressure")
u_comp = gpv_file.select(name="10 metre U wind component")
v_comp = gpv_file.select(name="10 metre V wind component")
# 予想時刻と抽出した要素データを格納するリストを初期化する
valid_dates = []
solar_list = []
temp_list = []
rh_list = []
dew_list = []
wspd_list = []
wetb_list = []
for sw, tk, rh, psfc, uwind, vwind in zip(shortwave, temp, rhum, pres, u_comp,v_comp) :
# 予想時刻をリストに抽出する
vt = self.extractValidTime(uwind)
# 日射量を抽出する
sw_pt = self.extractPointData(sw)
# 風速ベクトルを抽出してする
u_pt = self.extractPointData(uwind)
v_pt = self.extractPointData(vwind)
# 風速をMetpyで算出する
wspd = mpcalc.wind_speed(
u_pt * units('m/s'), v_pt * units('m/s')
)
# 湿球温度をMetpyで算出する
psfc_pt = self.extractPointData(psfc)
tc_pt = self.extractPointData(tk) - 273.15
rh_pt = self.extractPointData(rh)
dew_pt = mpcalc.dewpoint_from_relative_humidity(
tc_pt * units.degC,
rh_pt * units.percent
)
wetb_pt = mpcalc.wet_bulb_temperature(
psfc_pt * units.Pa,
tc_pt * units.degC,
dew_pt.magnitude * units.degC
)
# リスト変数に追加する
valid_dates.append(vt)
solar_list.append(sw_pt)
temp_list.append(tc_pt)
rh_list.append(rh_pt)
dew_list.append(dew_pt.magnitude)
wspd_list.append(wspd.magnitude)
wetb_list.append(wetb_pt.magnitude)
# 黒球温度を計算する
So = np.array(solar_list)
U = np.array(wspd_list)
Ta = np.array(temp_list)
Tg = ((So - 38.5) / ((So * 0.0217) + (U * 4.35) + 23.5)) + Ta
# WBGTを計算する
Tw = np.array(wetb_list)
WBGT_out = (0.7 * Tw) + (0.2 * Tg) + (0.1 * Ta)
return (valid_dates, WBGT_out.tolist(),
temp_list, dew_list, rh_list,
Tg.tolist(), Tw.tolist(),
solar_list, wspd_list
)露点温度、湿球温度および風速の物理量計算を行うため、metpy.calcをインポートしエイリアス（mpcalc）を設定しています。
MetPyでは、計算対象の物理量に対して単位を設定しなければならないため、metpy.unitsモジュールもインポートし気象物理量の単位を設定できるようします。
今回は、黒球温度を自力で計算するため、Numpyもインポートしています。
import pygrib
import metpy.calc as mpcalc
from metpy.units import units
from datetime import timedelta
import numpy as np次に、コンストラクタ変数の初期化を行います。
続いて予報時刻を抽出するメソッド（extractValidTime）とポイントのデータを抽出するメソッド（extractPointData）を作成します。これらのメソッドは、WBGTを計算するメソッド内で使用します。
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 extractValidTime(self, physquantity) :
validtime = physquantity.validDate + self.diff_time
return validtime
def extractPointData(self, physquantity) :
pointdata = physquantity.data(
lat1 = self.req_lat1, lat2 = self.req_lat2,
lon1 = self.req_lon1, lon2 = self.req_lon2
)[0][0][0]
return pointdataここからがWBGTを算出する中心的なメソッドです。
PygribでGPVをオープンして、各GRIBメッセージを日射量（shortwave）、地上気温（temp）、地上相対湿度（rhum）、地上気圧（pres）および地上の成分風速（u_comp、v_comp）のオブジェクト変数に格納します。
予想時刻と抽出した物理量データを格納するリスト変数を初期化しておきます。
def calcWBGT(self) :
# GPVファイル名を指定してオープンする
gpv_file = pygrib.open(self.grib2dat)
# 日射量、気温、相対湿度、風速ベクトルを取り出す。
shortwave = gpv_file.select(name="Mean surface downward short-wave radiation flux")
temp = gpv_file.select(name="Temperature")
rhum = gpv_file.select(name="Relative humidity")
pres = gpv_file.select(name="Surface pressure")
u_comp = gpv_file.select(name="10 metre U wind component")
v_comp = gpv_file.select(name="10 metre V wind component")
# 予想時刻と抽出した要素データを格納するリストを初期化する
valid_dates = []
solar_list = []
temp_list = []
rh_list = []
dew_list = []
wspd_list = []
wetb_list = []以下の部分は、各物理量のGRIBメッセージを格納したオブジェクト変数をループさせて、予報時刻、特定地点の日射量を抽出しています。
スカラー風速を求めるため、抽出した特定地点の成分風速をMetpyに与えて計算させます。
湿球温度を求めるため、はじめに抽出した特定地点の地上気温と相対湿度からMetpyで露点温度を計算します。次に求められた露点温度と抽出した特定地点の地上気圧および地上気温からMetpyで湿球温度を計算させています。
for sw, tk, rh, psfc, uwind, vwind in zip(shortwave, temp, rhum, pres, u_comp,v_comp) :
# 予想時刻をリストに抽出する
vt = self.extractValidTime(uwind)
# 日射量を抽出する
sw_pt = self.extractPointData(sw)
# 風速ベクトルを抽出してする
u_pt = self.extractPointData(uwind)
v_pt = self.extractPointData(vwind)
# 風速をMetpyで算出する
wspd = mpcalc.wind_speed(
u_pt * units('m/s'), v_pt * units('m/s')
)
# 湿球温度をMetpyで算出する
psfc_pt = self.extractPointData(psfc)
tc_pt = self.extractPointData(tk) - 273.15
rh_pt = self.extractPointData(rh)
dew_pt = mpcalc.dewpoint_from_relative_humidity(
tc_pt * units.degC,
rh_pt * units.percent
)
wetb_pt = mpcalc.wet_bulb_temperature(
psfc_pt * units.Pa,
tc_pt * units.degC,
dew_pt.magnitude * units.degC
)
# リスト変数に追加する
valid_dates.append(vt)
solar_list.append(sw_pt)
temp_list.append(tc_pt)
rh_list.append(rh_pt)
dew_list.append(dew_pt.magnitude)
wspd_list.append(wspd.magnitude)
wetb_list.append(wetb_pt.magnitude)最後に、黒球温度とWBGTを自力で計算します。
リスト変数に格納された特定地点の日射量などの物理量は、このままでは計算できないため、Numpyアレイに代入してから計算します。
黒球温度は、日射量（SO）、風速（U）および地上気温（Ta）を黒球温度推定式に代入して求めます。
屋外におけるWBGT（WBGT_out）は、湿球温度（Tw）黒球温度（Tg）および気温（Ta）のそれぞれに係数を掛けてたものの総和により求まります。
このクラスからは、予想時間、WBGT、地上気温、露点温度、相対湿度、黒球温度、湿球温度、日射量、地上風速が戻り値に設定されています。
なお、戻り値はすべてリスト変数で返したいので、WBGT、黒球温度、湿球温度はtolist()関数でNumpyアレイからリスト変数に変換しています。
# 黒球温度を計算する
So = np.array(solar_list)
U = np.array(wspd_list)
Ta = np.array(temp_list)
Tg = ((So - 38.5) / ((So * 0.0217) + (U * 4.35) + 23.5)) + Ta
# WBGTを計算する
Tw = np.array(wetb_list)
WBGT_out = (0.7 * Twb) + (0.2 * Tg) + (0.1 * Ta)
return (valid_dates, WBGT_out.tolist(),
temp_list, dew_list, rh_list,
Tg.tolist(), Tw.tolist(),
solar_list, wspd_list
)WBGTを計算してみる
作成したクラスをインポートして、WBGTを計算してみます。Pythonのメインスクリプトは以下のとおりです。
from MetPhysModule import MetPhysModule
import pandas as pd
# 東京駅の緯度経度を設定する
land_lat = 35.681236
land_lon = 139.767125
# 気象庁MSMのGRIB2ファイル名を指定する
msm_gpv ='Z__C_RJTD_20230720000000_MSM_GPV_Rjp_Lsurf_FH00-15_grib2.bin'
# MetPhysModuleクラスをインスタンス化する
calcphys = MetPhysModule(land_lat, land_lon, msm_gpv)
# WBGTを計算する
vt, wbgt, tempc, dewp, rhum, tg, tw, sw, wspd = calcphys.calcWBGT()
# 予報時刻、戻り値のリスト変数をPandasのデータフレームに入力する
df = pd.DataFrame({
"ValidTime":vt,
"WBGT":wbgt,
"Temp C":tempc,
"Dewp C":dewp,
"Rhum %":rhum,
"Globe C":tg,
"WetBulb C":tw,
"Solar W/m**2":sw,
"Wind m/s":wspd
})
pd.set_option('display.max_rows', None)
pd.set_option('display.max_columns',None)
print(df)コードを実行すると以下のような結果が出力されます。
(MetPy) [reafnex@devsrv01 MetPy]$ python3.12 main.py
ValidTime WBGT Temp C Dewp C Rhum % Globe C \
0 2023-07-20 09:00:00 28.204719 27.333521 23.901170 81.552700 40.728725
1 2023-07-20 10:00:00 29.551712 29.249536 24.068106 73.676754 44.488780
2 2023-07-20 11:00:00 30.218298 29.818384 24.615339 73.675938 45.704081
3 2023-07-20 12:00:00 30.596670 31.011743 24.213162 67.163734 47.045105
4 2023-07-20 13:00:00 30.655998 31.878809 23.856524 62.573605 47.087742
5 2023-07-20 14:00:00 29.846093 31.465326 24.391543 66.153625 42.265909
6 2023-07-20 15:00:00 28.795842 30.401056 24.341985 70.092070 38.516961
7 2023-07-20 16:00:00 27.473617 29.820123 23.533158 69.029978 34.662994
8 2023-07-20 17:00:00 25.833955 28.832086 22.308251 67.857021 30.757492
9 2023-07-20 18:00:00 24.255273 27.586572 21.292224 68.571854 26.998699
10 2023-07-20 19:00:00 23.761430 26.878900 21.272135 71.392384 25.572846
11 2023-07-20 20:00:00 23.653113 26.396204 21.556524 74.745361 25.034782
12 2023-07-20 21:00:00 23.521129 25.999078 21.721967 77.300816 24.536084
13 2023-07-20 22:00:00 23.371278 25.621515 21.801221 79.435814 24.129072
14 2023-07-20 23:00:00 23.254207 25.258295 21.911953 81.721951 23.789356
WetBulb C Solar W/m**2 Wind m/s
0 24.750889 581.136551 1.011287
1 25.327147 654.402016 0.624138
2 25.850919 713.884396 0.810105
3 25.837820 730.350897 0.874044
4 25.786527 672.188646 0.822773
5 26.066255 480.616901 1.610380
6 25.789062 344.372106 1.543703
7 25.084295 205.157766 1.485307
8 23.998926 102.255630 1.699741
9 22.995537 19.375000 1.979792
10 22.798530 0.000000 1.374277
11 22.866480 0.000000 1.098681
12 22.877149 0.000000 0.647333
13 22.833303 0.000000 0.527962
14 22.815008 0.000000 0.6228522023年の関東甲信地方の梅雨明けは7月22日で、7月20日は日本の南海上に梅雨前線が停滞しており、東京は曇りで最高気温が32.6℃でした。
WBGTの計算結果から、この日の日中は厳重警戒（WBGT=28～31）の範囲であり、熱中症の危険が高い状態であったことが考えられます。
まとめ
いかがでしたでしょうか。
今回は、PythonのPygribやMetPyモジュールを利用して、暑さ指標であるWBGT（湿球黒球温度）を計算するコードを紹介しました。
また、WBGTを求めるうえで重要となる黒球温度については、一般的に観測によって得られるものですが、数日間の期間における暑さ指数の予想を可能とするため、気象庁MSMの数値予報結果とM.Okadaらの研究によって考案された黒球温度推定式を用いて予測するためのコードをPythonで作成しました。
近年における数値予報の予測精度は非常に高精度化しているため、求められたWBGTも高い精度を有している可能性が考えられます。
今回紹介したWBGTの計算方法には風速が考慮されているため、人体影響と相関関係にある黒球温度の推定は非常によく実現できていると考えられます。
しかしながら、降水、散乱による日射量変化や地面等から受ける赤外放射の影響など細かなところは考慮されていないため、一定の誤差を持っていることを十分に理解して利用する必要があります。
観測と予報の両側面から、効率的かつ効果的なWBGTの運用を目指し、地球温暖化により増加し続けている熱中症リスクから人類を守っていことが、SDGsの目標１３「気候変動に具体的な対策を」に対する具体的なアクションとなると信じています。
参考になれば幸いです。
気象データ利用のお悩みについてご相談ください
気象システムの設計製造及び運用保守にかかわる現役システムエンジニアが、気象予報士の視点で気象データを活用したビジネスアイディアの創出と課題解決のお手伝いをさせていただきます。
気象や波浪などの数値予報モデルを活用した気象予報システム開発や運用についても是非ご相談ください。
