Skew-T LogP plot by MetPy

はじめに

大気の鉛直構造を把握するため、日々、高層気象観測が行われている。これらは観測としてのみならず、数値予報モデルの初期値・解析値としても用いられる重要な観測データである。本稿では、その観測データを地点毎に把握するための Skew-T LogP plot について述べる。ここでは、MetPy で提供されているものを利用し、特に理想化した数値実験の初期値に用いる一点鉛直プロファイルの可視化について紹介する。

データ

とある時刻・地点のデータを利用する。
1行目は、地上の気圧 (hPa)・温位 (K)・水蒸気混合比 (g/kg)。
2行目以降は、高度 (m)・温位 (K)・水蒸気混合比 (g/kg)・東西風 (m/s)・南北風 (m/s)。
これを、input_sounding という名前で保存。

 1000.00  298.81   18.08
  277.04  299.83   16.13   -5.29    1.54
  504.73  299.94   16.02   -5.96    1.84
  736.79  300.06   16.00   -6.08    2.16
  973.81  301.25   15.48   -5.57    3.11
 1465.20  303.95   14.38   -2.57    5.48
 1982.15  306.88   12.90    1.13    5.17
 3104.82  312.11    9.51    5.20    4.45
 4369.76  318.33    6.41    5.66    4.46
 5823.29  325.50    4.00   11.84   -0.45
 7542.22  334.36    1.78   15.32    3.46
 9655.91  343.85    0.65   15.49   -5.92
12000.00  351.53    0.11   15.49   -5.92
21000.00  500.53    0.01   15.49   -5.92

ソースコード

上記のデータを、下記で可視化する。

import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
import matplotlib
import matplotlib.dates as mdates

import metpy.calc as mpcalc
from metpy.calc import temperature_from_potential_temperature, height_to_pressure_std, dewpoint_from_specific_humidity
from metpy.plots import SkewT
from metpy.units import units, pandas_dataframe_to_unit_arrays
import pint_xarray

# ヘッダーを取り除いた上でカラムを指定、pandas dataframe として格納。
df = pd.read_fwf('input_sounding', header=0, usecols=[0, 1, 2, 3, 4], names=['height', 'theta', 'qv', 'uwnd', 'vwnd'])

# pandas dataframe を xarray dataset に変換
ds = df.to_xarray()

# MetPy.units を用い、単位を指定。
ds['height'] = ds['height'] * units.meter
ds['theta'] = ds['theta'] * units.kelvin
ds['qv'] = ds['qv'] * units('g/kg')
ds['uwnd'] = ds['uwnd'] * units.meter_per_second
ds['vwnd'] = ds['vwnd'] * units.meter_per_second

# SkewT logP plot の描画に必要な変数に変換。
ds['pressure'] = height_to_pressure_std(ds['height']).pint.to('hPa')
ds['temperature'] = temperature_from_potential_temperature(ds['pressure'], ds['theta']).pint.to('degC')
ds['dewpoint'] = dewpoint_from_specific_humidity(ds['pressure'], ds['temperature'], ds['qv']).pint.to('degC')

# ある高度から持ち上げたパーセルの気温・露点温度のプロファイルを計算。
parcel_prof = mpcalc.parcel_profile(ds['pressure'], ds['temperature'][0], ds['dewpoint'][0]).pint.to('degC')

# 持ち上げ凝結高度の気圧・気温・露点温度を計算。
lcl_pressure, lcl_temperature = mpcalc.lcl(ds['pressure'][0], ds['temperature'][0], ds['dewpoint'][0])


# 描画
fig = plt.figure(figsize=(8, 6))
skew = SkewT(fig)
skew.ax.set_ylim(1050,100)
skew.ax.set_xlim(-30,35)
skew.ax.set_title("A sample sounding")
skew.plot(ds['pressure'], ds['temperature'], color="k", linestyle='-')
skew.plot(ds['pressure'], ds['dewpoint'], linestyle='-')
skew.shade_cape(ds['pressure'], ds['temperature'], parcel_prof)
skew.plot(ds['pressure'], parcel_prof, linestyle=':', linewidth=2)
skew.plot(lcl_pressure, lcl_temperature, 'ko', markerfacecolor='black')
skew.plot_moist_adiabats(linewidth=0.5)
fig.savefig('skewt.png', bbox_inches='tight')


# 各種指数・パラメーターの計算
p = ds['pressure']
T = ds['temperature']
Td = ds['dewpoint']
height = ds['height']
u = ds['uwnd']
v = ds['vwnd']

ctotals = mpcalc.cross_totals(p, T, Td)
kindex = mpcalc.k_index(p, T, Td)
showalter = mpcalc.showalter_index(p, T, Td)
total_totals = mpcalc.total_totals_index(p, T, Td)
vert_totals = mpcalc.vertical_totals(p, T)
lift_index = mpcalc.lifted_index(p, T, parcel_prof)
cape, cin = mpcalc.cape_cin(p, T, Td, parcel_prof)
lclp, lclt = mpcalc.lcl(p[0], T[0], Td[0])
lfcp, _ = mpcalc.lfc(p, T, Td)
el_pressure, _ = mpcalc.el(p, T, Td, parcel_prof)
ml_t, ml_td = mpcalc.mixed_layer(p, T, Td, depth=50 * units.hPa)
ml_p, _, _ = mpcalc.mixed_parcel(p, T, Td, depth=50 * units.hPa)
mlcape, mlcin = mpcalc.mixed_layer_cape_cin(p, T, parcel_prof, depth=50 * units.hPa)
mu_p, mu_t, mu_td, _ = mpcalc.most_unstable_parcel(p, T, Td, depth=50 * units.hPa)
mucape, mucin = mpcalc.most_unstable_cape_cin(p, T, Td, depth=50 * units.hPa)
(u_storm, v_storm), *_ = mpcalc.bunkers_storm_motion(p, u, v, height)
critical_angle = mpcalc.critical_angle(p, u, v, height, u_storm, v_storm)


# 各種指数・パラメーターの値を確認
print(f'        CAPE: {cape:.2f}')
print(f'         CIN: {cin:.2f}')
print(f'LCL Pressure: {lclp:.2f}')
print(f'LFC Pressure: {lfcp:.2f}')
print(f' EL Pressure: {el_pressure:.2f}')
print()
print('Mixed Layer - Lowest 50-hPa')
print(f'     ML Temp: {ml_t:.2f}')
print(f'     ML Dewp: {ml_td:.2f}')
print(f'     ML CAPE: {mlcape:.2f}')
print(f'      ML CIN: {mlcin:.2f}')
print()
print('Most Unstable - Lowest 50-hPa')
print(f'     MU Temp: {mu_t:.2f}')
print(f'     MU Dewp: {mu_td:.2f}')
print(f' MU Pressure: {mu_p:.2f}')
print(f'     MU CAPE: {mucape:.2f}')
print(f'      MU CIN: {mucin:.2f}')
print()
print(f'   Lifted Index: {lift_index:.2f}')
print(f'        K-Index: {kindex:.2f}')
print(f'Showalter Index: {showalter:.2f}')
print(f'   Cross Totals: {ctotals:.2f}')
print(f'   Total Totals: {total_totals:.2f}')
print(f'Vertical Totals: {vert_totals:.2f}')
print()
print('Bunkers Storm Motion Vector')
print(f'  u_storm: {u_storm:.2f}')
print(f'  v_storm: {v_storm:.2f}')
print(f'Critical Angle: {critical_angle:.2f}')

出力結果

図が Skew-T logP plot。横軸は気温または露点温度 (degree C)、縦軸は気圧 (hPa)。黒実線で気温、青実線で露点温度、橙点線で持ち上げたパーセルの気温をそれぞれ示す、赤色の陰影で示される面積は、持ち上げたパーセルと周辺大気との気温差、すなわち対流有効位置エネルギー (CAPE) である。黒丸印は持ち上げ凝結高度 (LCL) を示す。
図の左下から右上に向かう灰色の線は等温線を、図の右下から左上に伸びる青色点線は等湿潤断熱線をそれぞれ示す。

図の 925 hPa から少なくとも 500 hPa までは、気温と露点温度の線が重なっている。これは、ほぼ相対湿度 100% であることを示す。加えて、800 hPa から 600 hPa までの層では、気温減率の傾きが湿潤断熱線の傾きよりも大きくなっており、湿潤絶対不安定 (Moist Absolutely Unstable Layer, Bryan and Fritsch 2000) の可能性がある。

下記は、一点鉛直プロファイルから得られる指数やパラメーターの出力である。一地点の鉛直プロファイルから、これだけ多様な診断パラメーターを算出することが可能となる。

        CAPE: 729.54 joule / kilogram
         CIN: -3.21 joule / kilogram
LCL Pressure: 926.10 hectopascal
LFC Pressure: 885.54 hectopascal
 EL Pressure: 260.51 hectopascal

Mixed Layer - Lowest 50-hPa
     ML Temp: 22.89 degree_Celsius
     ML Dewp: 20.64 degree_Celsius
     ML CAPE: 1997.43 joule / kilogram
      ML CIN: -0.15 joule / kilogram

Most Unstable - Lowest 50-hPa
     MU Temp: 20.55 degree_Celsius
     MU Dewp: 20.17 degree_Celsius
 MU Pressure: 927.77 hectopascal
     MU CAPE: 754.89 joule / kilogram
      MU CIN: -1.08 joule / kilogram

   Lifted Index: -3.48 delta_degree_Celsius
        K-Index: 41.32 degree_Celsius
Showalter Index: -3.25 delta_degree_Celsius
   Cross Totals: 24.40 delta_degree_Celsius
   Total Totals: 48.64 delta_degree_Celsius
Vertical Totals: 24.24 delta_degree_Celsius

Bunkers Storm Motion Vector
  u_storm: 0.25 meter_per_second
  v_storm: -4.63 meter_per_second
Critical Angle: 179.85 degree

指数・診断パラメーターは、発案された地域や特性に依存するため、計算結果をそのまま異なる地域に当てはめることは得策ではない。上記で調べたように、可視化することで指数が意図した計算となっているか、愚直に確認することが重要である (時に時間を要するが)。ここでは、計算方法の紹介に留める。

おわりに

大気の状態を把握するための基礎的な観測の一つである、高層気象観測データの可視化と関連する指数・パラメータの計算について紹介した。

参考文献

• May, R. M., Goebbert, K. H., Thielen, J. E., Leeman, J. R., Camron, M. D., Bruick, Z., Bruning, E. C., Manser, R. P., Arms, S. C., and Marsh, P. T., 2022: MetPy: A Meteorological Python Library for Data Analysis and Visualization. Bull. Amer. Meteor. Soc., 103, E2273-E2284, doi:10.1175/BAMS-D-21-0125.1.
• Bryan, G. H., and J. M. Fritsch, 2000: Moist Absolute Instability: The Sixth Static Stability State. Bull. Amer. Meteor. Soc., 81, 1207–1230, doi:10.1175/1520-0477(2000)081<1287:MAITSS>2.3.CO;2.

更新履歴

• 2024-01-01: 初稿

Miscellaneous — Jan 1, 2024