はじめに

気象関係の要素を計算する関数、例えば気温と相対湿度から露点温度を求める関数などを集めた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)など

気象要素を計算するPythonモジュールを公開します

はじめに

インストール

ライセンス

主な機能

使い方

風を扱う関数

UV_to_SpdDir(U, V)

SpdDir_to_UV(Wspd, Wdir)

Deg_to_Dir8(Wdir, dir_zero=None, numeric=False)

Deg_to_Dir16(val, dir_zero=None, numeric=False)

Cross_Wind(Wspd, Wdir, RWY)

Tail_Wind(Wspd, Wdir, RWY)

Head_Wind(Wspd, Wdir, RWY)

大気中の水分量を扱う関数

RH_to_Td(T, RH, formula="Bolton")

Td_to_RH(T, Td, formula="Bolton")

T_to_WVP(T, formula="Bolton")

WVP_to_T(es, formula="Bolton")

T_Td(T, Td)

Mixing_Ratio(Td, P, formula="Bolton")

Specific_Humidity(Td, P, formula="Bolton")

Absolute_Humidity(T, Td, formula="Bolton")

Virtual_Temperature(T, Td, P, formula="Bolton")

大気の熱力学や安定度に係る関数

Theta(T, P)

Tlcl(T, Td)

Theta_e(T, Td, P, formula="Bolton")

SSI(P0, P1, T0, T1, Td0, formula="Bolton", **kwargs)

K_Index(T850, Td850, T700, Td700, T500)

PRES_to_PRMSL(P1, T1, Z1)

Surface_Height(P0, P1, T1)

単位変換

MPS_to_KT(x)

KT_to_MPS(x)

M_to_FT(x)

FT_to_M(x)

degF_to_degC(x)

degC_to_degF(x)

今後のバージョンアップについて