はじめに
気象関係の要素を計算する関数、例えば気温と相対湿度から露点温度を求める関数などを集めたPythonモジュール wxparams を作ったので、公開します。
（Weather Parametersの略）
Pythonを始める前はPerlで気象データを処理するスクリプトをよく書いていて、やはり気象関係でよく使う関数をPerlモジュールにまとめていました。コロナ禍の影響で家で過ごす時間が増えたので(笑)、この機にPythonに書き換えようと思ったのがきっかけです。
どうせ作るなら公開してしまおうということで、ソースコードはGitHub、コード例などの使い方はnote、マニュアルとして関数の説明などはこの記事にまとめました。
- Ver1.0 : 2020/04/24
- Ver1.1 : 2020/05/27
- 気温の単位：摂氏 <=> 華氏の変換を追加
- Ver1.2 : 2020/10/03
-
WVP_to_T関数を改修（ゼロの割り算が発生する場合に対応）
-
- Ver1.3 : 2020/11/20
- 比湿・絶対湿度・仮温度を求める関数を追加
- Ver1.4 : 2021/01/16
-
SSI関数を改修、海面気圧・地上の高度を求める関数を追加
-
- Ver1.5 : 2021/01/17
-
Theta_e関数で使われている$ R_{d}\ /\ C_{pd} $の値を修正
-
インストール
私のGitHubからpipインストールできます。
pip install git+https://github.com/Yoshiki443/weather_parameters
pipインストールに抵抗がある方は、GitHubのwxparams/wxparams.pyだけダウンロードして使っていただいても良いです。
ライセンス
MITライセンスとしています。
主な機能
- 風を扱う関数
- 風のUV成分 <=> 風向・風速の変換、8方位・16方位風向、cross windなど
- 大気中の水分量を扱う関数
- 相対湿度 <=> 露点温度の変換、飽和水蒸気圧、混合比など
- 大気の熱力学や安定度に係る関数
- 温位、相当温位、SSI、K-Indexなど
- 単位変換
- m/s <=> **knot** の変換, **meter** <=> feet の変換など
使い方
コード例を先にご紹介します。
import numpy as np
import wxparams as wx
temp = np.array([[0., 5.],[10., 20.]]) # 気温[C]
rh = np.array([[90., 50.], [70., 99.5]]) # 湿度[%]
td = wx.RH_to_Td(temp, rh) # 露点温度[C]
print(td)
# 出力
# [[-1.44330606 -4.56523582]
# [ 4.78251527 19.91913689]]
まずimportします。コード例にある省略名「wx」はweatherの略としてよく使われるものです。importしたら関数を呼び出します。上の例では、気温と相対湿度から露点温度を計算する場合を示しています。
ここで入力するデータは基本的に numpy.ndarray を想定しています。もしくは pandas.Series でも動作します。個々の気象要素を1つずつ入力した場合、動作する関数と動作しない関数があります。
入力データとして、例えばGPVのある気圧面のある気象要素（つまり2次元配列）を丸ごと入力して別の気象要素を計算したり、csvなど構造化データのある気象要素の列を入力して別の気象要素を計算したり、といった使い方を想定しています。
一部の処理で numpy.ndarray を前提としたものもあるので、作成者としては numpy.ndarray を推奨します。
以下、各関数の説明です。
風を扱う関数
UV_to_SpdDir(U, V)
風のUV成分から、風向・風速を計算します。
UV成分とは、U＝東西風・V＝南北風を表し、気象庁MSMなど数値予報GPVの風データは通常UV成分になっています。一般には東西風なら西風が正の値、南北風なら南風が正の値になっています。
風速の単位は m/s, knot など速度の単位なら何でもOKです。
風向は真北が0度ではなく、360度となります。風速が0だった場合に風向も0度となります。
引数(Parameters) :
- U : array_like
- 風のU成分（東西風成分）
- V : array_like
- 風のV成分（南北風成分）
戻り値(Returns) :
- Wspd : array_like
- 風速。単位は入力データと同じ
- Wdir : array_like
- 風向[degree]
SpdDir_to_UV(Wspd, Wdir)
UV_to_SpdDirの逆で、風向風速から風のUV成分を計算します。風速の単位は m/s, knot など速度の単位なら何でもOKです。
引数(Parameters) :
- Wspd : array_like
- 風速[m/s, kt, ...etc]
- Wdir : array_like
- 風向[degree]
戻り値(Returns) :
- U : array_like
- 風のU成分（東西風成分）
- V : array_like
- 風のV成分（南北風成分）
Deg_to_Dir8(Wdir, dir_zero=None, numeric=False)
風向を360度から8方位に変換します。
戻り値は「北・北東・東・南東・南・南西・西・北西」の8方位を表すアルファベット「N・NE・E・SE・S・SW・W・NW」に変換されます。
風向0の場合は、引数dir_zeroで指定した文字列に変換されます。デフォルトはNoneですが、例えば航空気象の場合に、風向が定まらない状態を表す「VRB」に変換する使い方ができます。
また引数numericをTrueにすると、数値「8・1・2・3・4・5・6・7」で出力します。このとき真北は0ではなく8となります。風向0の場合に出力0となります。
引数(Parameters) :
- Wdir : array_like
- 360度の風向[degree]
- dir_zero : str, optional
- 風向0の場合に変換する文字列を指定。デフォルトはNone
- numeric : bool, optional
- 8方位を表すアルファベットで出力されるか(False)、数値で出力させるか(True)を指定。デフォルトはFalse
戻り値(Returns) :
- Dir8 : array_like
- 8方位表記の風向
Deg_to_Dir16(val, dir_zero=None, numeric=False)
Deg_to_Dir8のように、風向を360度から16方位に変換します。
すなわち戻り値は「北・北北東・北東・東北東・東・東南東・南東・南南東・南・南南西・南西・西南西・西・西北西・北西・北北西」の16方位を表すアルファベット「N・NNE・NE・ENE・E・ESE・SE・SSE・S・SSW・SW・WSW・W・WNW・NW・NNW」に変換されます。
また引数numericをTrueにすると、数値「16・1・2・3・4・5・6・7・8・9・10・11・12・13・14・15」で出力します。風向0の扱いはDeg_to_Dir8と同じです。
引数(Parameters) :
- Wdir : array_like
- 360度の風向[degree]
- dir_zero : str, optional
- 風向0の場合に変換する文字列を指定。デフォルトはNone
- numeric : bool, optional
- 16方位を表すアルファベットで出力されるか(False)、数値で出力させるか(True)を指定。デフォルトはFalse
戻り値(Returns) :
- Dir16 : array_like
- 16方位表記の風向
Cross_Wind(Wspd, Wdir, RWY)
風速・風向・空港の滑走路の向きから、風の横風成分(Cross Wind)を計算します。
風速の単位は m/s, knot など速度の単位なら何でもOKです。風向の単位は degree です。
滑走路は2桁のRWY番号ではなく、360度で与えます。例えばRWY22であれば、220度を入力します。
引数(Parameters) :
- Wspd : array_like
- 風速[m/s, kt, ...etc]
- Wdir : array_like
- 風向[degree]
- RWY : array_like
- 滑走路の向き（360度方位）
戻り値(Returns) :
- crswind : array_like
- 風の横風成分の風速
Tail_Wind(Wspd, Wdir, RWY)
風速・風向・空港の滑走路の向きから、風の追い風成分(Tail Wind)を計算します。その他は Cross_Wind に同じです。
引数(Parameters) :
- Wspd : array_like
- 風速[m/s, kt, ...etc]
- Wdir : array_like
- 風向[degree]
- RWY : array_like
- 滑走路の向き（360度方位）
戻り値(Returns) :
- tailwind : array_like
- 風の追い風成分の風速
Head_Wind(Wspd, Wdir, RWY)
風速・風向・空港の滑走路の向きから、風の向かい風成分(Head Wind)を計算します。Tail_Windの正負を反対にした値に一致します。その他、Cross_Wind に同じです。
引数(Parameters) :
- Wspd : array_like
- 風速[m/s, kt, ...etc]
- Wdir : array_like
- 風向[degree]
- RWY : array_like
- 滑走路の向き（360度方位）
戻り値(Returns) :
- headwind : array_like
- 風の向かい風成分の風速
大気中の水分量を扱う関数
RH_to_Td(T, RH, formula="Bolton")
気温[C]と相対湿度[%]から、露点温度[C]を計算します。相対湿度は0%だと計算ができないので、0.1%が最小値となるように変換されます。
計算過程で飽和水蒸気圧[hPa]の計算をしますが、計算式が3種類あります。デフォルトはBoltonの式を使用しますが、オプションでTetensの式、WMO近似式も指定可能です。
詳しくはT_to_WVPを参照して下さい。
引数(Parameters) :
- T : array_like
- 気温[C]
- RH : array_like
- 相対湿度[%]
- formula : str, optional (default="Bolton")
- 飽和水蒸気圧[hPa]の計算式。"Bolton" の他に、"Tetens", "WMO"を指定可能
戻り値(Returns) :
- Td : array_like
- 露点温度[C]
Td_to_RH(T, Td, formula="Bolton")
RH_to_Tdの逆で、気温[C]と露点温度[C]から、相対湿度[%]を計算します。
その他はRH_to_Tdと同じです。
引数(Parameters) :
- T : array_like
- 気温[C]
- Td : array_like
- 露点温度[C]
- formula : str, optional (default="Bolton")
- 飽和水蒸気圧[hPa]の計算式。"Bolton" の他に、"Tetens", "WMO"を指定可能
戻り値(Returns) :
- RH : array_like
- 相対湿度[%]
T_to_WVP(T, formula="Bolton")
気温[C]から飽和水蒸気圧[hPa]を計算します。気温の代わりに露点温度[C]を入力すれば、水蒸気圧[hPa]が計算されます。
計算式は3種類あります。デフォルトはBoltonの式を使います。相当温位を計算する際に、気象庁が採用している計算式で実装しているためです。Theta_eも参照して下さい。
ただ私が確認した限り、実用上はどの式を用いても大差ないので、どれでも好みで良いと思います。個人的にはTetensの式になじみがあります。
Tetensの式 : $$ {\Large es = 6.1078 \times 10^{\left(\frac{7.5 T}{T + 237.3}\right)} } $$
WMO近似式 : $$ {\Large es = e^{\left(19.482 - \frac{4303.4}{T + 243.5}\right)} } $$
Boltonの式 : $$ {\Large es = 6.112 \times e^{\left(\frac{17.67 T}{T + 243.5}\right)} } $$
引数(Parameters) :
- T : array_like
- 気温[C]（もしくは露点温度[C]）
- formula : str, optional (default="Bolton")
- 飽和水蒸気圧[hPa]の計算式。"Bolton" の他に、"Tetens", "WMO"を指定可能
戻り値(Returns) :
- es : array_like
- 飽和水蒸気圧[hPa]（露点温度を入力した場合は水蒸気圧[hPa]）
WVP_to_T(es, formula="Bolton")
T_to_WVPの逆関数で、飽和水蒸気圧[hPa]から気温[C]を計算します。水蒸気圧[hPa]を入力すれば露点温度[C]が計算されます。
引数(Parameters) :
- es : array_like
- 飽和水蒸気圧[hPa]（もしくは水蒸気圧[hPa]）
- formula : str, optional (default="Bolton")
- 飽和水蒸気圧[hPa]の計算式。"Bolton" の他に、"Tetens", "WMO"を指定可能
戻り値(Returns) :
- T : array_like
- 気温[C]（水蒸気圧を入力した場合は露点温度[C]）
T_Td(T, Td)
気温[C]と露点温度[C]から、湿数[C]を計算します。単純に引き算を実行するだけの関数です。
引数(Parameters) :
- T : array_like
- 気温[C]
- Td : array_like
- 露点温度[C]
戻り値(Returns) :
- T-Td : array_like
- 湿数[C]
Mixing_Ratio(Td, P, formula="Bolton")
露点温度[C]と気圧[hPa]から、混合比[g/g]を計算します。
計算過程で飽和水蒸気圧を計算するので、計算式をformulaオプションで指定できます。オプションの種類はT_to_WVPを参照して下さい。
引数(Parameters) :
- Td : array_like
- 露点温度[C]
- P : array_like
- 気圧[hPa]
- formula : str, optional (default="Bolton")
- 飽和水蒸気圧[hPa]の計算式。"Bolton" の他に、"Tetens", "WMO"を指定可能
戻り値(Returns) :
- w : array_like
- 混合比[g/g]
Specific_Humidity(Td, P, formula="Bolton")
露点温度[C]と気圧[hPa]から、比湿[g/g]を計算します。
計算過程で飽和水蒸気圧を計算するので、計算式をformulaオプションで指定できます。オプションの種類はT_to_WVPを参照して下さい。
引数(Parameters) :
- Td : array_like
- 露点温度[C]
- P : array_like
- 気圧[hPa]
- formula : str, optional (default="Bolton")
- 飽和水蒸気圧[hPa]の計算式。"Bolton" の他に、"Tetens", "WMO"を指定可能
戻り値(Returns) :
- q : array_like
- 比湿[g/g]
Absolute_Humidity(T, Td, formula="Bolton")
気温[C]と露点温度[C]から、絶対湿度[$ g/m^3 $]を計算します。絶対湿度とは、要するに水蒸気密度のことです。
計算過程で飽和水蒸気圧を計算するので、計算式をformulaオプションで指定できます。オプションの種類はT_to_WVPを参照して下さい。
引数(Parameters) :
- T : array_like
- 気温[C]
- Td : array_like
- 露点温度[C]
- formula : str, optional (default="Bolton")
- 飽和水蒸気圧[hPa]の計算式。"Bolton" の他に、"Tetens", "WMO"を指定可能
戻り値(Returns) :
- $ \rho_w $ : array_like
- 絶対湿度[$ g/m^3 $]
Virtual_Temperature(T, Td, P, formula="Bolton")
気温[C]、露点温度[C]、気圧[hPa]から、仮温度[C]を計算します。
計算過程で飽和水蒸気圧を計算するので、計算式をformulaオプションで指定できます。オプションの種類はT_to_WVPを参照して下さい。
引数(Parameters) :
- T : array_like
- 気温[C]
- Td : array_like
- 露点温度[C]
- P : array_like
- 気圧[hPa]
- formula : str, optional (default="Bolton")
- 飽和水蒸気圧[hPa]の計算式。"Bolton" の他に、"Tetens", "WMO"を指定可能
戻り値(Returns) :
- $ T_v $ : array_like
- 仮温度[C]
大気の熱力学や安定度に係る関数
Theta(T, P)
気温[C]と気圧[hPa]から、温位[K]を計算します。
引数(Parameters) :
- T : array_like
- 気温[C]
- P : array_like
- 気圧[hPa]
戻り値(Returns) :
- Theta : array_like
- 温位[K]
Tlcl(T, Td)
気温[K]と露点温度[K]から、持ち上げ凝結高度の気温[K]を計算します。気温と露点温度の単位は[K]で入力する点に注意して下さい。
相当温位の計算過程で使っています。計算式についてはTheta_eも参照して下さい。
$$ {\Large T_{LCL} = \frac{1}{\frac{1}{Td-56}+\frac{ln(T\ /\ Td)}{800}}+56 } $$
引数(Parameters) :
- T : array_like
- 気温[K]
- Td : array_like
- 露点温度[K]
戻り値(Returns) :
- Tlcl : array_like
- 持ち上げ凝結高度の気温[K]
Theta_e(T, Td, P, formula="Bolton")
気温[C]・露点温度[C]・気圧[hPa]から、相当温位[K]を計算します。
相当温位の計算式も複数あるようですが、ここでは気象庁が採用した算出方法に準じています。詳しくはリンク先PDFの最後のページを参照して下さい。
$$ {\large \theta_{e} = T\left(\frac{1000}{P-e}\right)^\frac{R_{d}}{C_{pd}}\left(\frac{T}{T_{LCL}}\right)^{0.28w}exp\left(\left(\frac{3036.0}{T_{LCL}}-1.78\right)w(1+0.448w)\right) } $$
なお上記PDFだとRd/Cpdの値が0.2854となっていますが、通常は0.2857だと思います。なぜ0.2854なのか、誤植なのか or 理由があるのかまでは把握できていません。実装では0.2857を採用しました。
オリジナルの論文 Bolton (1980) にて、$ R_{d}\ /\ C_{pd} = 0.2854 $ となっているとの情報を頂きましたので、この値に修正しました。自分で確認した限りでは、この修正により0.01〜0.1[K]程度のオーダーで計算値に変化があります。
計算過程で飽和水蒸気圧を計算するので、計算式をformulaオプションで指定できます。デフォルトは気象庁の相当温位計算式に倣って"Bolton"としています。オプションの種類はT_to_WVPを参照して下さい。
引数(Parameters) :
- T : array_like
- 気温[C]
- Td : array_like
- 露点温度[C]
- P : array_like
- 気圧[hPa]
- formula : str, optional (default="Bolton")
- 飽和水蒸気圧[hPa]の計算式。"Bolton" の他に、"Tetens", "WMO"を指定可能
戻り値(Returns) :
- Thetae : array_like
- 相当温位[K]
SSI(P0, P1, T0, T1, Td0, formula="Bolton", **kwargs)
大気不安定指数の代表格、SSIを計算します。
入力は持ち上げる空気の高度の気圧[hPa]・気温[C]・露点温度[C]と、持ち上げた先の高度の気圧[hPa]・気温[C]です。通常SSIは850-500hPaで計算しますが、任意の気圧面での計算が可能です。
SSIは、P0高度のパーセル（空気塊）を断熱的にP1高度まで持ち上げたときの気温をTlとしたとき、T1 - Tlで計算されます。この関数ではTlを以下のように求めています。
- 測高公式を用いて、P0-P1間の層厚を計算します。ここでキーワード引数
h0とh1でジオポテンシャル高度を与えると、h1 - h0で層厚を計算します（※） - P0高度のパーセルを乾燥断熱減率でP1高度まで持ち上げた場合の気温（Tl_dry）と、持ち上げ凝結高度の気温（Tlcl）を計算します
-
Tl_dry > Tlclの場合、パーセルはP1高度まで飽和しないとみなし、Tl_dryをTlとして使います - それ以外はパーセルがP1高度までに飽和するとみなし、探索的にTlを推定します（詳細はソースコードをご覧下さい）
※ 層厚の計算はジオポテンシャル高度を与えた方が正確ですが、実用上は測高公式での推定で十分だと考えています。与える引数が少なくて済みますし、層厚が必要になるのはパーセルが持ち上げた先まで飽和しない場合のみなので、SSIの実用上もあまり問題になりません。
計算過程で飽和水蒸気圧を計算するので、計算式をformulaオプションで指定できます。オプションの種類はT_to_WVPを参照して下さい。
引数(Parameters) :
- P0 : array_like
- 持ち上げる元の高度の気圧[hPa]
- P1 : array_like
- 持ち上げた先の高度の気圧[hPa]
- T0 : array_like
- 持ち上げる元の高度の気温[C]
- T1 : array_like
- 持ち上げた先の高度の気温[C]
- Td : array_like
- 持ち上げる元の高度の露点温度[C]
- formula : str, optional (default="Bolton")
- 飽和水蒸気圧[hPa]の計算式。"Bolton" の他に、"Tetens", "WMO"を指定可能
- h0 : array_like
- 持ち上げる元のジオポテンシャル高度[m]
- キーワード引数で与えるオプションです
- h1 : array_like
- 持ち上げた先のジオポテンシャル高度[m]
- キーワード引数で与えるオプションです
戻り値(Returns) :
- SSI : array_like
- SSI
K_Index(T850, Td850, T700, Td700, T500)
大気不安定指数の1つ、K-Indexを計算します。
入力は850hPaの気温[C]と露点温度[C]、700hPaの気温[C]と露点温度[C]、500hPaの気温[C]です。
引数(Parameters) :
- T850 : array_like
- 850hPaの気温[C]
- Td850 : array_like
- 850hPaの露点温度[C]
- T700 : array_like
- 700hPaの気温[C]
- Td700 : array_like
- 700hPaの露点温度[C]
- T500 : array_like
- 500hPaの気温[C]
戻り値(Returns) :
- Kindex : array_like
- K-Index
PRES_to_PRMSL(P1, T1, Z1)
測高公式を用いて、海面気圧[hPa]を計算します。
入力は地上気圧[hPa]、地上気温[C]、地上の標高[m]です。
$$ {\large P_0 = P_1 \left( 1 - \frac{0.0065 Z_1}{T_1 + 273.15 + 0.0065 Z_1} \right)^{-5.257} } $$
引数(Parameters) :
- P1 : array_like
- 地上気圧[hPa]
- T1 : array_like
- 地上気温[C]
- Z1 : array_like
- 地上の標高[m]
戻り値(Returns) :
- P0 : array_like
- 海面気圧[hPa]
Surface_Height(P0, P1, T1)
PRES_to_PRMSLに記載した測高公式をZ1について解き、地上の標高[m]を計算します。
入力は海面気圧[hPa]、地上気圧[hPa]、地上気温[C]です。
引数(Parameters) :
- P0 : array_like
- 海面気圧[hPa]
- P1 : array_like
- 地上気圧[hPa]
- T1 : array_like
- 地上気温[C]
戻り値(Returns) :
- Z1 : array_like
- 地上の標高[m]
単位変換
MPS_to_KT(x)
風速など、速度の単位変換です。m/s を knot に変換します。
引数(Parameters) :
- x : array_like
- 速度[m/s]
戻り値(Returns) :
- y : array_like
- 速度[kt]
KT_to_MPS(x)
風速など、速度の単位変換です。knot を m/s に変換します。
引数(Parameters) :
- x : array_like
- 速度[kt]
戻り値(Returns) :
- y : array_like
- 速度[m/s]
M_to_FT(x)
高度など、長さの単位変換です。meter を feet に変換します。
引数(Parameters) :
- x : array_like
- 長さ[m]
戻り値(Returns) :
- y : array_like
- 長さ[ft]
FT_to_M(x)
高度など、長さの単位変換です。feet を meter に変換します。
引数(Parameters) :
- x : array_like
- 長さ[ft]
戻り値(Returns) :
- y : array_like
- 長さ[m]
degF_to_degC(x)
気温の単位変換です。華氏[F]を摂氏[C]に変換します。
引数(Parameters) :
- x : array_like
- 気温[F]
戻り値(Returns) :
- y : array_like
- 気温[C]
degC_to_degF(x)
気温の単位変換です。摂氏[C]を華氏[F]に変換します。
引数(Parameters) :
- x : array_like
- 気温[C]
戻り値(Returns) :
- y : array_like
- 気温[F]
今後のバージョンアップについて
最初に述べたとおり、もともとPerlモジュールとして作ったものを、一通りPythonに書き換えましたので、いったん完成としました。
今後、下記のような気象要素も取り入れていきます。需要があるのかわかりませんが(笑)、気象の専門家のこだわりで作っていきたいと思います。
- LFC、CAPE、CIN、Lifted Index など
- 移流、渦度、風の鉛直シア、安定度など
- 湿球温度、湿球温位など
- アメリカで使われているスーパーセル指数(SCP)や竜巻指数(STP)など