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

# はじめに
気象関係の要素を計算する関数、例えば気温と相対湿度から露点温度を求める関数などを集めたPythonモジュール **```wxparams```** を作ったので、公開します。
（**Weather Parameters**の略）

Pythonを始める前はPerlで気象データを処理するスクリプトをよく書いていて、やはり気象関係でよく使う関数をPerlモジュールにまとめていました。コロナ禍の影響で家で過ごす時間が増えたので(笑)、この機にPythonに書き換えようと思ったのがきっかけです。

どうせ作るなら公開してしまおうということで、ソースコードは[GitHub](https://github.com/Yoshiki443/weather_parameters)、コード例などの使い方は[note](https://note.com/calmtree/n/n6692a736f6c1)、マニュアルとして関数の説明などはこの記事にまとめました。

- 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](https://github.com/Yoshiki443/weather_parameters)の```wxparams/wxparams.py```だけダウンロードして使っていただいても良いです。


# ライセンス
MITライセンスとしています。


# 主な機能
1. 風を扱う関数
    - 風のUV成分 <=> 風向・風速の変換、8方位・16方位風向、cross windなど
2. 大気中の水分量を扱う関数
    - 相対湿度 <=> 露点温度の変換、飽和水蒸気圧、混合比など
3. 大気の熱力学や安定度に係る関数
    - 温位、相当温位、SSI、K-Indexなど
4. 単位変換
    - **m/s** <=> **knot** の変換, **meter** <=> **feet** の変換など


# 使い方
コード例を先にご紹介します。

```python
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]を計算します。
相当温位の計算式も複数あるようですが、ここでは[気象庁が採用した算出方法](https://www.data.jma.go.jp/add/suishin/jyouhou/pdf/371.pdf)に準じています。詳しくはリンク先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)](https://doi.org/10.1175/1520-0493(1980)108%3C1046:TCOEPT%3E2.0.CO;2) にて、$ 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を以下のように求めています。

1. 測高公式を用いて、P0-P1間の層厚を計算します。ここでキーワード引数`h0`と`h1`でジオポテンシャル高度を与えると、`h1 - h0`で層厚を計算します（※）
2. P0高度のパーセルを乾燥断熱減率でP1高度まで持ち上げた場合の気温（Tl_dry）と、持ち上げ凝結高度の気温（Tlcl）を計算します
3. `Tl_dry > Tlcl`の場合、パーセルはP1高度まで飽和しないとみなし、Tl_dryをTlとして使います
4. それ以外はパーセルがP1高度までに飽和するとみなし、探索的にTlを推定します（詳細はソースコードをご覧下さい）

<font color="gray">※ 層厚の計算はジオポテンシャル高度を与えた方が正確ですが、実用上は測高公式での推定で十分だと考えています。与える引数が少なくて済みますし、層厚が必要になるのはパーセルが持ち上げた先まで飽和しない場合のみなので、SSIの実用上もあまり問題になりません。</font>

計算過程で飽和水蒸気圧を計算するので、計算式を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に書き換えましたので、いったん完成としました。

今後、下記のような気象要素も取り入れていきます。需要があるのかわかりませんが(笑)、気象の専門家のこだわりで作っていきたいと思います。

1. LFC、CAPE、CIN、Lifted Index など
2. 移流、渦度、風の鉛直シア、安定度など
2. 湿球温度、湿球温位など
3. アメリカで使われているスーパーセル指数(SCP)や竜巻指数(STP)など