EC 数据驱动的颠簸指数计算全解析(Python 实现)

# ===================== 路径配置 =====================EC_FILE = "/202603121200-upper-xxxxxxxxx.nc"OUTPUT_DIR = "/"GIF_PATH = "turbulence_all_levels.gif"# =======================================================import numpy as npimport netCDF4import matplotlibmatplotlib.use("AGG")import matplotlib.pyplot as pltimport osfrom PIL import Imageimport warningswarnings.filterwarnings('ignore')# ===================== 配置 =====================EARTH_RADIUS = 6378100.0os.makedirs(OUTPUT_DIR, exist_ok=True)# 中国范围裁切LON_MIN, LON_MAX = 70, 140LAT_MIN, LAT_MAX = 15, 55# ===================== 读取数据 =====================with netCDF4.Dataset(EC_FILE) as ds:    lon = ds.variables['longitude'][:]    lat = ds.variables['latitude'][:]    isobaric = ds.variables['isobaricInhPa'][:]    t = ds.variables['t'][0]    u = ds.variables['u'][0]    v = ds.variables['v'][0]    z = ds.variables['z'][0]# ===================== 核心算法：颠簸指数 =====================def compute_lp_fast(t, u, v, z, lat):    lat_rad = np.radians(lat)    dx = 0.25 * np.pi / 180 * EARTH_RADIUS * np.cos(lat_rad)    dy = 0.25 * np.pi / 180 * EARTH_RADIUS    dx = dx[:, None]    dz = z + 1e-6    dtx_full = np.zeros_like(t)    dty_full = np.zeros_like(t)    du_full = np.zeros_like(t)    dv_full = np.zeros_like(t)    dtx_full[:, :, 1:-1] = 0.5 * (t[:, :, 2:] - t[:, :, :-2])    dty_full[:, 1:-1, :] = 0.5 * (t[:, 2:, :] - t[:, :-2, :])    du_full[:, :, 1:-1] = 0.5 * (u[:, :, 2:] - u[:, :, :-2])    dv_full[:, 1:-1, :] = 0.5 * (v[:, 2:, :] - v[:, :-2, :])    dvdz = np.sqrt((du_full / dz * 100) ** 2 + (dv_full / dz * 100) ** 2)    dtdn = np.sqrt((dtx_full / dx * 1e5) ** 2 + (dty_full / dy * 1e5) ** 2)    dvdn = np.sqrt((du_full / dx * 1e5) ** 2 + (dv_full / dy * 1e5) ** 2)    li = 7.268 * dvdz + 0.718 * dtdn + 0.133 * dvdn - 2.52    lp = 1 / (1 + np.exp(-0.59 * li))    lp[:, 0, :] = lp[:, -1, :] = lp[:, :, 0] = lp[:, :, -1] = 0    return lpprint("🔢 计算颠簸指数...")lp = compute_lp_fast(t, u, v, z, lat)print("✅ 计算完成！开始极速绘图...")# ===================== 绘图配置 =====================lon2d, lat2d = np.meshgrid(lon, lat)levels = np.linspace(0, 1, 21)cmap = "jet"images = []# ===================== 批量绘图 =====================for idx, p in enumerate(isobaric):    fig, ax = plt.subplots(figsize=(9, 6), dpi=100)    # 绘制颠簸指数    cf = ax.contourf(lon2d, lat2d, lp[idx], levels=levels, cmap=cmap, extend="max")    # 只裁切中国区域，不画任何边界    ax.set_xlim(LON_MIN, LON_MAX)    ax.set_ylim(LAT_MIN, LAT_MAX)    ax.set_title(f"EC Turbulence Index | {p:.0f} hPa", fontsize=13)    # 底部色标，每张图都有    cax = fig.add_axes([0.15, 0.05, 0.7, 0.025])    fig.colorbar(cf, cax=cax, orientation="horizontal", label="Turbulence Index (LP)")    # 保存    img_path = f"{OUTPUT_DIR}/lp_{int(p):03d}.png"    plt.savefig(img_path, bbox_inches="tight")    plt.close()    images.append(img_path)    print(f"✅ {int(p):3d} hPa 已保存")# ===================== 生成 GIF =====================print("\n🎞️ 生成 GIF...")frames = [Image.open(f) for f in images]frames[0].save(    GIF_PATH,    save_all=True,    append_images=frames[1:],    duration=500,    loop=0)