代码实战 | 如何在垂直剖面图上加对流层顶

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from metpy.interpolate import cross_section

import skyborn
# SKYBORN 对流层顶计算（向量化版本，处理 2D/3D/4D 网格）
trop_wmo = skyborn.troposphere.trop_wmo

# 加载 ERA5 数据
file_path = 'ERA5-2023-08_pl.nc'
ds = xr.open_dataset(file_path)

# 定义剖面起点和终点（纬度, 经度）
start_point = (45, 100.0)   # (lat, lon)
end_point   = (20, 130.0)   # (lat, lon)

# 取第一个时刻，转换为 MetPy 的 CF 标准格式
ds_cf = ds.isel(time=0).metpy.parse_cf()

# 使用 MetPy 的 cross_section 函数提取剖面数据
cross = cross_section(ds_cf, start_point, end_point)
print(cross)

<xarray.Dataset> Size: 210kB
Dimensions:    (level: 37, index: 100)
Coordinates:
  * level      (level) int32 148B 1 2 3 5 7 10 20 ... 875 900 925 950 975 1000
  * index      (index) int64 800B 0 1 2 3 4 5 6 7 8 ... 92 93 94 95 96 97 98 99
    longitude  (index) float64 800B 100.0 100.4 100.8 ... 129.5 129.8 130.0
    latitude   (index) float64 800B 45.0 44.79 44.59 44.38 ... 20.56 20.28 20.0
    time       datetime64[ns] 8B 2023-08-02
    metpy_crs  object 8B Projection: latitude_longitude
Data variables:
    z          (level, index) float64 30kB 4.752e+05 4.751e+05 ... 445.0 507.9
    r          (level, index) float64 30kB -0.0005729 -0.0005729 ... 87.7 89.32
    t          (level, index) float64 30kB 259.9 259.9 260.1 ... 299.6 299.3
    q          (level, index) float64 30kB 4.049e-06 4.049e-06 ... 0.01905 0.019
    u          (level, index) float64 30kB -39.16 -38.9 -38.8 ... 8.306 7.847
    v          (level, index) float64 30kB 2.141 2.308 2.252 ... 13.12 13.05
    w          (level, index) float64 30kB -0.00102 -0.001254 ... 0.01719
Attributes:
    Conventions:  CF-1.6
    history:      2024-01-09 16:38:26 GMT by grib_to_netcdf-2.25.1: /opt/ecmw...

# 气压层（升序：1 -> 1000 hPa）
pressure = cross['level'].values
assert np.all(np.diff(pressure) > 0), '气压必须升序排列'

# 一次性向量化计算：trop_wmo 处理 2D (level, index) 数据
# levdim=0 (气压层), xdim=1 (剖面点), ydim=-1 (无纬度维)
result = trop_wmo(cross['t'].values, pressure, xdim=1, ydim=-1, levdim=0)
trop_pressure = result['pressure']   # (index,) 对流层顶气压 [hPa]
trop_height   = result['height']      # (index,) 对流层顶高度 [m]

print('对流层顶气压范围: {:.1f} ~ {:.1f} hPa'.format(np.nanmin(trop_pressure), np.nanmax(trop_pressure)))
print('对流层顶高度范围: {:.0f} ~ {:.0f} m'.format(np.nanmin(trop_height), np.nanmax(trop_height)))

对流层顶气压范围: 95.3 ~ 118.7 hPa
对流层顶高度范围: 15094 ~ 16488 m

# 字体设置
plt.rcParams['font.family'] = 'serif'
plt.rcParams['font.serif'] = ['DejaVu Serif', 'Times New Roman']
plt.rcParams['axes.unicode_minus'] = False

fig, ax = plt.subplots(figsize=(12, 8), dpi=150)

lon = cross['longitude'].values
lat = cross['latitude'].values
lev = cross['level'].values

# 温度转为摄氏度
temp_c = cross['t'].values - 273.15   # (level, index)

# --- 温度填色 ---
cf = ax.contourf(lon, lev, temp_c, levels=np.arange(-80, 41, 2),
                  cmap='RdYlBu_r', extend='both')
cbar = fig.colorbar(cf, ax=ax, pad=0.02)
cbar.set_label('Temperature (degC)', fontsize=12)

# --- 对流层顶线 ---
ax.plot(lon, trop_pressure, color='crimson', lw=2.5, ls='-',
        label='Tropopause (WMO, skyborn)')
ax.legend(loc='lower left', fontsize=11, framealpha=0.9)

# --- 纵轴：对数气压，高值在下 ---
ax.set_yscale('log')
ax.set_ylim(1000, 10)            # 1000 hPa 在底部, 10 hPa 在顶部
ax.set_yticks([1000, 700, 500, 300, 200, 100, 50, 20, 10])
ax.set_yticklabels(['1000', '700', '500', '300', '200', '100', '50', '20', '10'])

# --- 顶部副轴：纬度 ---
ax2 = ax.secondary_xaxis('top')
ax2.set_xticks(lon[::10])
ax2.set_xticklabels([f'{v:.1f}'for v in lat[::10]])
ax2.set_xlabel('Latitude (degN)', fontsize=12)

# --- 标签与样式 ---
ax.set_title('ERA5 Temperature Cross-Section with WMO Tropopause\n'
             f'{cross.time.values}  |  {start_point} -> {end_point}', fontsize=14)
ax.set_xlabel('Longitude (degE)', fontsize=12)
ax.set_ylabel('Pressure (hPa)', fontsize=12)
ax.tick_params(axis='both', labelsize=10)
ax.grid(True, ls='--', alpha=0.5)

plt.tight_layout()
plt.savefig('temperature_cross_section_tropopause.png', dpi=200, bbox_inches='tight')
plt.show()
print('图已保存: temperature_cross_section_tropopause.png')

图已保存: temperature_cross_section_tropopause.png