【MetPy】大気の鉛直構造を知るためのエマグラムの描画方法【Python】

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 = "20150924_09JST_47646.txt"
data = np.genfromtxt(fname=fname, unpack=True)

prss = data[0]*units.hPa # 気圧
hght = data[1]*units.meters # ジオポテンシャル高度
tmpr = data[2]*units.degC # 気温
rhmd = data[3] # 相対湿度
wndsp = data[4]*units('m/s') # 風速
wnddr = data[5]*units.degrees # 風向

#------------------------------------------------------------------------------
# --- 計算
# 風速，風向から風の水平成分を計算する。
uwnd, vwnd = mpcalc.wind_components(wndsp, wnddr)

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

#------------------------------------------------------------------------------
# *** プロット
# --- プロット範囲の設定
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'
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_barbs(prss, uwnd, vwnd)

# --- 乾燥断熱線のプロット
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 mr in mixing_ratio.flatten():
    p1 = 110 * units.hPa
    dewpt = mpcalc.dewpoint(mpcalc.vapor_pressure(p1, mr))
    skew.ax.text(dewpt, p1, f'{mr.magnitude}', fontsize=9,
                 horizontalalignment='center')
    

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

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 = "20150924_09JST_47646.txt"
data = np.genfromtxt(fname=fname, unpack=True)

prss = data[0]*units.hPa # 気圧
hght = data[1]*units.meters # ジオポテンシャル高度
tmpr = data[2]*units.degC # 気温
rhmd = data[3] # 相対湿度
wndsp = data[4]*units('m/s') # 風速
wnddr = data[5]*units.degrees # 風向

#------------------------------------------------------------------------------
# --- 計算
# 風速，風向から風の水平成分を計算する。
uwnd, vwnd = mpcalc.wind_components(wndsp, wnddr)

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

#------------------------------------------------------------------------------
# *** プロット
# --- プロット範囲の設定
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'
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_barbs(prss, uwnd, vwnd)

# --- 乾燥断熱線のプロット
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 mr in mixing_ratio.flatten():
    p1 = 110 * units.hPa
    dewpt = mpcalc.dewpoint(mpcalc.vapor_pressure(p1, mr))
    skew.ax.text(dewpt, p1, f'{mr.magnitude}', fontsize=9,
                 horizontalalignment='center')


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