我是气象绘图新手，十分感谢这款出色的工具！我正在制作供法国滑翔机使用的单日SkewT动画气象图，已拼凑出初步代码且能运行，但不确定质量如何。我已获取QNH气压，希望在对数刻度的气压轴上同步标注1000m、1500m、2000m、2500m、3000m及更高海拔线，却始终不知如何实现，恳请提供技术思路。

以下是我的测试代码：

#!/usr/bin/python3

import matplotlib.pyplot as plt
import xarray as xr
import metpy.calc as mpcalc
from metpy.plots import SkewT
from metpy.units import units
import pandas as pd
from metpy.plots import add_metpy_logo
from siphon.catalog import TDSCatalog
import numpy as np
from datetime import datetime, timedelta
import requests
import imageio.v2 as imageio
import os

## Parameters:
DaysAhead     = 4
PausePerHour  = 5000 #milliseconds
ParisPostCode = '75001'
StartTime     = 10
EndTime       = 15 ## one hour more than last time to compute

# --- 1. GEOLOCATION ---
def get_coords_france(query):
    try:
        url = f"https://api-adresse.data.gouv.fr/search/?q={query}&limit=1"
        response = requests.get(url).json()
        if response['features']:
            lon, lat = response['features'][0]['geometry']['coordinates']
            return lat, lon
    except Exception: pass
    return None, None

user_loc = input("Enter City or Postcode: ")
LAT, LON = get_coords_france(user_loc)
if LAT is None:
    print(f'{user_loc} not found, defaulting to {ParisPostCode}')
    user_loc = ParisPostCode

LAT, LON = get_coords_france(user_loc)

# --- 2. SETUP DATA CONNECTION ---
target_date = datetime.utcnow() + timedelta(days=DaysAhead)
cat_url = 'https://thredds.ucar.edu/thredds/catalog/grib/NCEP/GFS/Global_0p25deg/catalog.xml'
catalog = TDSCatalog(cat_url)
ds_name = [n for n in catalog.datasets.keys() if 'Best' in n][0]
ncss = catalog.datasets[ds_name].subset()

filenames = []

print(f"Generating animation for {user_loc}... ({LAT},{LON}) on {target_date.strftime('%Y %m %d')} {StartTime}:00 - {EndTime-1}:00")

for hour in range(StartTime,EndTime):
    target_time = datetime(target_date.year, target_date.month, target_date.day, hour, 0)
    
    try:
        # --- QUERY 1: 3D SOUNDING DATA ---
        q3d = ncss.query().lonlat_point(LON, LAT).time(target_time)
        q3d.accept('netcdf4').variables('Temperature_isobaric', 'Relative_humidity_isobaric', 
                                      'u-component_of_wind_isobaric', 'v-component_of_wind_isobaric')
        ds3d = xr.open_dataset(xr.backends.NetCDF4DataStore(ncss.get_data(q3d))).metpy.parse_cf()

        # --- QUERY 2: 2D SURFACE DATA ---
        q2d = ncss.query().lonlat_point(LON, LAT).time(target_time)
        q2d.accept('netcdf4').variables('Pressure_reduced_to_MSL_msl')
        ds2d = xr.open_dataset(xr.backends.NetCDF4DataStore(ncss.get_data(q2d))).metpy.parse_cf()

        # Extract QNH
        qnh_hpa = float(ds2d['Pressure_reduced_to_MSL_msl']) / 100.0

        # --- EXTRACTION (The emma13.py logic) ---
        T_var = ds3d['Temperature_isobaric'].squeeze()
        p_raw = T_var.metpy.vertical.values.flatten()
        p = (p_raw * units.hPa) if np.max(p_raw) < 2000 else (p_raw * units.Pa)
        T = T_var.values.flatten() * units.K
        rh = ds3d['Relative_humidity_isobaric'].squeeze().values.flatten() * units.percent
        u = ds3d['u-component_of_wind_isobaric'].squeeze().values.flatten() * units('m/s')
        v = ds3d['v-component_of_wind_isobaric'].squeeze().values.flatten() * units('m/s')

        # --- MONOTONIC FILTERING (The emma13.py logic) ---
        p_hpa = p.to('hPa').magnitude
        idx = np.argsort(p_hpa)[::-1]
        p_s, T_s, rh_s, u_s, v_s = p_hpa[idx], T[idx], rh[idx], u[idx], v[idx]
        
        keep_idx = [0]
        for i in range(1, len(p_s)):
            if p_s[i] < p_s[keep_idx[-1]]: keep_idx.append(i)
        
        p_f = p_s[keep_idx] * units.hPa
        t_f = T_s[keep_idx].to('degC')
        rh_f = rh_s[keep_idx]
        u_f = u_s[keep_idx]
        v_f = v_s[keep_idx]

        # --- CALCULATIONS ---
        rh_f[rh_f < 0.01 * units.percent] = 0.01 * units.percent
        td_f = mpcalc.dewpoint_from_relative_humidity(t_f, rh_f)
        prof = mpcalc.parcel_profile(p_f, t_f[0], td_f[0]).to('degC')
        cape, cin = mpcalc.cape_cin(p_f, t_f, td_f, prof)

        # --- PLOTTING ---
        fig = plt.figure(figsize=(10, 10))
        skew = SkewT(fig, rotation=45)
        
        # Bold Isotherms
        for val in [0,10,20,30]:
            skew.ax.axvline(val, color='blue', linestyle='-', linewidth=1.5, alpha=0.3)

        skew.plot(p_f, t_f, 'r', linewidth=2, label='Temperature')
        skew.plot(p_f, td_f, 'g', linewidth=2, label='Point de Rosé')
        skew.plot(p_f, prof, 'k--', alpha=0.5, label='Parcel')
        
        skew.shade_cape(p_f, t_f, prof)
        skew.shade_cin(p_f, t_f, prof, td_f)
        
        # Barb fix: ensure all arrays are same size via the [::2] slice
        skew.plot_barbs(p_f[::2], u_f[::2], v_f[::2])

        skew.plot_dry_adiabats(alpha=0.15)
        skew.plot_moist_adiabats(alpha=0.15)
        skew.plot_mixing_lines(alpha=0.15)

        # Add the MetPy logo!
        fig = plt.gcf()
        add_metpy_logo(fig, 150, 120, size='small');

        skew.ax.set_ylim(1050, 100)
        skew.ax.set_xlim(-40, 40)
        plt.title(f"{user_loc.upper()} | {target_time.strftime('%Y %m %d %H:%M')} UTC\nQNH: {qnh_hpa:.1f} hPa | CAPE/CIN: {cape.m:.0f}/{cin.m:.0f}")
        plt.legend(loc='upper right')
        
        fname = f"frame_{hour:02d}.png"
        plt.savefig(fname)
        filenames.append(fname)
        plt.close(fig)
        print(f"Processed {hour:02d}:00 ...")

    except Exception as e:
        print(f"Error at {hour:02d}:00: {e}")

# --- 5. STITCH ---
print("Finalizing GIF...")

file_name = f"forecast_{user_loc}_{target_date.strftime('%Y_%m_%d')}.gif"

with imageio.get_writer(file_name, mode='I', duration=PausePerHour) as writer:
    for f in filenames:
        writer.append_data(imageio.imread(f))
        os.remove(f)

print(f"Success! Created {file_name}")

核心逻辑 #

要在对数气压轴上标注海拔线，需将海拔高度转换为对应的气压值，然后在Skew-T图上绘制水平参考线并添加标注。利用MetPy内置的气压-海拔转换工具即可实现，无需手动编写复杂公式。

具体步骤 #

1. 定义需要标注的海拔列表：根据需求设置要显示的海拔值（单位：米）。
2. 海拔转气压：使用mpcalc.height_to_pressure_std函数，将海拔转换为标准大气下的气压值（单位：hPa）。该函数基于国际标准大气模型，适合滑翔机飞行参考的通用场景。
3. 绘制水平参考线：对每个转换后的气压值，在Skew-T图的y轴（气压轴）上绘制水平虚线，并在图的右侧添加海拔标注，避免遮挡曲线。

修改后的绘图代码片段 #

在你的# --- PLOTTING ---区块中，添加以下代码（建议放在skew.ax.set_xlim(-40, 40)之后，plt.title(...)之前）：

# 添加海拔参考线
elevations = [1000, 1500, 2000, 2500, 3000, 3500, 4000]  # 可按需扩展
for elev in elevations:
    # 将海拔转换为标准大气下的气压
    pressure = mpcalc.height_to_pressure_std(elev * units.meter).to('hPa').magnitude
    # 绘制水平虚线
    skew.ax.axhline(pressure, color='gray', linestyle='--', alpha=0.7)
    # 在图右侧标注海拔
    skew.ax.text(skew.ax.get_xlim()[1] + 1, pressure, f'{elev}m', 
                va='center', ha='left', color='gray', fontsize=9)

额外说明 #

• 如果需要更精确的气压值（结合实际探空温度廓线），可使用mpcalc.height_to_pressure_agl函数，但需要传入地面温度和气压参数：

  # 示例：结合实际探空数据计算
  surface_pressure = p_f[0]  # 探空的地面气压
  surface_temp = t_f[0]      # 探空的地面温度
  pressure = mpcalc.height_to_pressure_agl(elev * units.meter, surface_pressure, surface_temp).to('hPa').magnitude
  

• 调整elevations列表可添加更高或更低的海拔线，修改color、linestyle参数可自定义参考线样式。

内容的提问来源于stack exchange，提问作者gratefulfrog