【MetPy】大気の安定度指標を鉛直プロファイルから計算する【気象データ】

import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.plots import SkewT
from metpy.units import units
import numpy as np

#------------------------------------------------------------------------------
# データを読み込む
fname = "20170705_09JST_47807.txt"
data = np.genfromtxt(fname=fname, unpack=True)

prss_pcl = 950 * units.hPa # 持ち上げる空気塊がある気圧

prss = data[0]*units.hPa # 気圧
tmpr = data[1]*units.degC # 気温
rhmd = data[2] # 相対湿度

#------------------------------------------------------------------------------
# --- 計算
# 気温と相対湿度から露点温度を計算する。
dewp = mpcalc.dewpoint_from_relative_humidity(tmpr, rhmd/100.0)

# 持ち上げる空気塊の気圧のインデックスを取得
i_prss_pcl = np.argmin(np.abs(prss_pcl.magnitude - prss.magnitude))

# 空気塊を持ち上げた場合の温度変化を計算する
prss_wlcl, tmpr_wlcl, dewp_wlcl, prof_wlcl = mpcalc.parcel_profile_with_lcl(prss[i_prss_pcl:], tmpr[i_prss_pcl:], dewp[i_prss_pcl:])

# LCL (Lifted Condensation Level: 持ち上げ凝結高度)
lcl_prss, lcl_tmpr = mpcalc.lcl(prss[i_prss_pcl], tmpr[i_prss_pcl], dewp[i_prss_pcl])

# LFC (Level of Free Convection: 自由対流高度) を計算する
lfc_prss, lfc_tmpr = mpcalc.lfc(prss[i_prss_pcl:], tmpr[i_prss_pcl:], dewp[i_prss_pcl:])

# EL (Equilibrium Level: 平衡高度) を計算する
el_prss, el_tmpr = mpcalc.el(prss, tmpr, dewp)

# CAPE (convective available potential energy: 対流有効位置エネルギー) と
# CIN (convective inhibition: 対流抑制) を計算する
cape, cin = mpcalc.surface_based_cape_cin(prss, tmpr, dewp)

# SSI (Shawalter Stability Index: ショワルター安定指数) を計算する
ssi = mpcalc.showalter_index(prss, tmpr, dewp)

# KI (K-index: K指数) を計算する
ki = mpcalc.k_index(prss, tmpr, dewp)

#------------------------------------------------------------------------------
# *** プロット
# --- プロット範囲の設定
xmax = 40 # 横軸（温度）の最大値
xmin = -90 # 横軸（温度）の最小値
ymin = 100 # 縦軸（気圧）の最小値
ymax = 1020 # 縦軸（気圧）の最大値

# --- 全般の設定
fig = plt.figure(figsize=(210/25.4, 294/25.4), dpi=100)
out_fig_path = 'emagram.png'

# --- 安定度の指数を書く。値に応じて安定度のコメント、色を変える
fig.text(0.05, 0.95, f"starting parcel pressure = {prss[i_prss_pcl].magnitude:.0f} hPa")
fig.text(0.05, 0.93, f"LCL = {lcl_prss.magnitude:.0f} hPa")
fig.text(0.05, 0.91, f"LFC = {lfc_prss.magnitude:.0f} hPa")
fig.text(0.05, 0.89, f"EL = {el_prss.magnitude:.0f} hPa")

fig.text(0.55, 0.95, f"CIN = {cin.magnitude:.0f} J/kg")
fig.text(0.55, 0.93, f"CAPE = {cape.magnitude:.0f} J/kg")
fig.text(0.55, 0.91, f"SSI = {ssi[0].magnitude:.0f}℃")
fig.text(0.55, 0.89, f"KI = {int(ki.magnitude):.0f}℃")

# --- 図の設定 
skew = SkewT(fig, rotation=0, aspect=150)

# --- 軸の設定
skew.ax.set_xlim(xmin, xmax)
skew.ax.set_ylim(ymax, ymin)

# --- ラベルの設定
skew.ax.set_xlabel('Temperature (℃)')
skew.ax.set_ylabel('Pressure (hPa)')

# --- 気温，露点温度，風の鉛直分布をプロットする
skew.plot(prss, tmpr, 'red')
skew.plot(prss, dewp, 'blue')
skew.plot(prss_wlcl, prof_wlcl, 'black')

# --- LCL、LFCおよびELをプロット
skew.plot(lcl_prss, lcl_tmpr, 'ko', markerfacecolor='black')
skew.plot(lfc_prss, lfc_tmpr, 'ko', markerfacecolor='black')

# --- 乾燥断熱線のプロット
t0 = np.arange(200, 360, 10)*units.K # プロットする乾燥断熱線のスタート温度

# 乾燥断熱線をプロット
skew.plot_dry_adiabats(t0=t0, lw=0.5, colors='red')

# ラベルを350hPaに付す
for t in t0:
    p1 = 350*units.hPa
    t1 = mpcalc.dry_lapse(p1, t, 1000*units.hPa) # p1=350hPaまで乾燥断熱線に従って持ち上げたときの温度
    if t1 < xmin*units.degC:
        continue
    skew.ax.text(t1, p1, f'{t.magnitude}', fontsize=9,
                 horizontalalignment='center', color='red')

# --- 湿潤断熱線のプロット
t0 = np.arange(245, 320, 5)*units.K # プロットする湿潤断熱線のスタート温度

# 湿潤断熱線のプロット
skew.plot_moist_adiabats(t0=t0, lw=0.5, colors='green')

# ラベルを250hPaに付す
for t in t0:
    p1 = 250 * units.hPa
    t1 = mpcalc.moist_lapse(p1, t, 1000*units.hPa) # p1=250hPaまで湿潤断熱線に従って持ち上げたときの温度
    if t1 < xmin*units.degC:
        continue
    skew.ax.text(t1, p1, f'{t.magnitude}', fontsize=9,
                 horizontalalignment='center', color='green')

# --- 等飽和混合比線のプロット
# プロットする混合比の値，プロットする気圧の範囲を指定する
mixing_ratio = np.array([0.2, 0.4, 0.6, 1, 2, 4, 5, 10, 15, 20, 25, 30, 35]).reshape(-1, 1) * units('g/kg')
pressure = np.arange(1000, 10, -50)*units.hPa

# 等飽和混合比線のプロット
skew.plot_mixing_lines(mixing_ratio=mixing_ratio, pressure=pressure, lw=0.5, colors='blue')

# ラベルを110hPa付近に付す
for i, mr in enumerate(mixing_ratio.flatten()):
    p1 = 110 * units.hPa + i%2 * 10 * units.hPa
    dewpt = mpcalc.dewpoint(mpcalc.vapor_pressure(p1, mr))
    if mr.magnitude < 1.0:
        skew.ax.text(dewpt, p1, f'{mr.magnitude:.1f}', fontsize=9,
                    horizontalalignment='center', color="blue")
    else:
        skew.ax.text(dewpt, p1, f'{mr.magnitude:.0f}', fontsize=9,
                    horizontalalignment='center', color="blue")
        

# CINとCAPEに相当する部分をそれぞれ青と赤で塗りつぶす
skew.shade_cin(prss_wlcl, tmpr_wlcl, prof_wlcl)
skew.shade_cape(prss_wlcl, tmpr_wlcl, prof_wlcl)

#------------------------------------------------------------------------------
# --- 図の保存
plt.savefig(out_fig_path, transparent=False)

$ cd ~/Desktop/labcode/python/plot-emagram-and-calc-params

$ python plot-emagram-and-calc-params.py

import matplotlib.pyplot as plt
import metpy.calc as mpcalc
from metpy.plots import SkewT
from metpy.units import units
import numpy as np

#------------------------------------------------------------------------------
# データを読み込む
fname = "20170705_09JST_47807.txt"
data = np.genfromtxt(fname=fname, unpack=True)

prss_pcl = 950 * units.hPa # 持ち上げる空気塊がある気圧

prss = data[0]*units.hPa # 気圧
tmpr = data[1]*units.degC # 気温
rhmd = data[2] # 相対湿度

#------------------------------------------------------------------------------
# --- 計算
# 気温と相対湿度から露点温度を計算する。
dewp = mpcalc.dewpoint_from_relative_humidity(tmpr, rhmd/100.0)

# 持ち上げる空気塊の気圧のインデックスを取得
i_prss_pcl = np.argmin(np.abs(prss_pcl.magnitude - prss.magnitude))

# 空気塊を持ち上げた場合の温度変化を計算する
prss_wlcl, tmpr_wlcl, dewp_wlcl, prof_wlcl = mpcalc.parcel_profile_with_lcl(prss[i_prss_pcl:], tmpr[i_prss_pcl:], dewp[i_prss_pcl:])

# LCL (Lifted Condensation Level: 持ち上げ凝結高度)
lcl_prss, lcl_tmpr = mpcalc.lcl(prss[i_prss_pcl], tmpr[i_prss_pcl], dewp[i_prss_pcl])

# LFC (Level of Free Convection: 自由対流高度) を計算する
lfc_prss, lfc_tmpr = mpcalc.lfc(prss[i_prss_pcl:], tmpr[i_prss_pcl:], dewp[i_prss_pcl:])

# EL (Equilibrium Level: 平衡高度) を計算する
el_prss, el_tmpr = mpcalc.el(prss, tmpr, dewp)

# CAPE (convective available potential energy: 対流有効位置エネルギー) と
# CIN (convective inhibition: 対流抑制) を計算する
cape, cin = mpcalc.surface_based_cape_cin(prss, tmpr, dewp)

# SSI (Shawalter Stability Index: ショワルター安定指数) を計算する
ssi = mpcalc.showalter_index(prss, tmpr, dewp)

# KI (K-index: K指数) を計算する
ki = mpcalc.k_index(prss, tmpr, dewp)

#------------------------------------------------------------------------------
# *** プロット
# --- プロット範囲の設定
xmax = 40 # 横軸（温度）の最大値
xmin = -90 # 横軸（温度）の最小値
ymin = 100 # 縦軸（気圧）の最小値
ymax = 1020 # 縦軸（気圧）の最大値

# --- 全般の設定
fig = plt.figure(figsize=(210/25.4, 294/25.4), dpi=100)
out_fig_path = 'emagram.png'

# --- 安定度の指数を書く。値に応じて安定度のコメント、色を変える
fig.text(0.05, 0.95, f"starting parcel pressure = {prss[i_prss_pcl].magnitude:.0f} hPa")
fig.text(0.05, 0.93, f"LCL = {lcl_prss.magnitude:.0f} hPa")
fig.text(0.05, 0.91, f"LFC = {lfc_prss.magnitude:.0f} hPa")
fig.text(0.05, 0.89, f"EL = {el_prss.magnitude:.0f} hPa")

fig.text(0.55, 0.95, f"CIN = {cin.magnitude:.0f} J/kg")
fig.text(0.55, 0.93, f"CAPE = {cape.magnitude:.0f} J/kg")
fig.text(0.55, 0.91, f"SSI = {ssi[0].magnitude:.0f}℃")
fig.text(0.55, 0.89, f"KI = {int(ki.magnitude):.0f}℃")

# --- 図の設定 
skew = SkewT(fig, rotation=0, aspect=150)

# --- 軸の設定
skew.ax.set_xlim(xmin, xmax)
skew.ax.set_ylim(ymax, ymin)

# --- ラベルの設定
skew.ax.set_xlabel('Temperature (℃)')
skew.ax.set_ylabel('Pressure (hPa)')

# --- 気温，露点温度，風の鉛直分布をプロットする
skew.plot(prss, tmpr, 'red')
skew.plot(prss, dewp, 'blue')
skew.plot(prss_wlcl, prof_wlcl, 'black')

# --- LCL、LFCおよびELをプロット
skew.plot(lcl_prss, lcl_tmpr, 'ko', markerfacecolor='black')
skew.plot(lfc_prss, lfc_tmpr, 'ko', markerfacecolor='black')

# --- 乾燥断熱線のプロット
t0 = np.arange(200, 360, 10)*units.K # プロットする乾燥断熱線のスタート温度

# 乾燥断熱線をプロット
skew.plot_dry_adiabats(t0=t0, lw=0.5, colors='red')

# ラベルを350hPaに付す
for t in t0:
    p1 = 350*units.hPa
    t1 = mpcalc.dry_lapse(p1, t, 1000*units.hPa) # p1=350hPaまで乾燥断熱線に従って持ち上げたときの温度
    if t1 < xmin*units.degC:
        continue
    skew.ax.text(t1, p1, f'{t.magnitude}', fontsize=9,
                 horizontalalignment='center', color='red')

# --- 湿潤断熱線のプロット
t0 = np.arange(245, 320, 5)*units.K # プロットする湿潤断熱線のスタート温度

# 湿潤断熱線のプロット
skew.plot_moist_adiabats(t0=t0, lw=0.5, colors='green')

# ラベルを250hPaに付す
for t in t0:
    p1 = 250 * units.hPa
    t1 = mpcalc.moist_lapse(p1, t, 1000*units.hPa) # p1=250hPaまで湿潤断熱線に従って持ち上げたときの温度
    if t1 < xmin*units.degC:
        continue
    skew.ax.text(t1, p1, f'{t.magnitude}', fontsize=9,
                 horizontalalignment='center', color='green')

# --- 等飽和混合比線のプロット
# プロットする混合比の値，プロットする気圧の範囲を指定する
mixing_ratio = np.array([0.2, 0.4, 0.6, 1, 2, 4, 5, 10, 15, 20, 25, 30, 35]).reshape(-1, 1) * units('g/kg')
pressure = np.arange(1000, 10, -50)*units.hPa

# 等飽和混合比線のプロット
skew.plot_mixing_lines(mixing_ratio=mixing_ratio, pressure=pressure, lw=0.5, colors='blue')

# ラベルを110hPa付近に付す
for i, mr in enumerate(mixing_ratio.flatten()):
    p1 = 110 * units.hPa + i%2 * 10 * units.hPa
    dewpt = mpcalc.dewpoint(mpcalc.vapor_pressure(p1, mr))
    if mr.magnitude < 1.0:
        skew.ax.text(dewpt, p1, f'{mr.magnitude:.1f}', fontsize=9,
                    horizontalalignment='center', color="blue")
    else:
        skew.ax.text(dewpt, p1, f'{mr.magnitude:.0f}', fontsize=9,
                    horizontalalignment='center', color="blue")
        

# CINとCAPEに相当する部分をそれぞれ青と赤で塗りつぶす
skew.shade_cin(prss_wlcl, tmpr_wlcl, prof_wlcl)
skew.shade_cape(prss_wlcl, tmpr_wlcl, prof_wlcl)

#------------------------------------------------------------------------------
# --- 図の保存
plt.savefig(out_fig_path, transparent=False)