台风路径判定-初始场的可视化

台风路径相关初始场

对于台风路径的追踪，ECMWF（European Centre for Medium-Range Weather Forecasts）在早些年就做出了相关规定。

给定初始时间和位置
给定准备时间（6 小时），在 445km 范围内找到当地平均海平面气压 (MSLP) 的最小值（小于1015hpa），并同时满足以下条件：

半径 278km 范围内有一个 850hPa 涡度最大值，绝对值大于 5×10^(-5)

半径 278km 范围内 850hPa 和 200hPa 之间存在厚度最大

当气旋位于陆地上时，半径 278km 范围内的最大 10m 风速大于 8m/s

当所处地形海拔高于1000m时停止跟踪。

[A] P. White, Newsletter No. 102 – Winter 2004/05, by ECMWF, 2005.

数据准备 

对于初始场数据，可以在欧洲气候数据存取中心（Climate Data Store, CDS）的官网进行下载，对于以上的判断条件，我们需要下载single level（地表）的海平面气压（mslp），10m风速（10u、10mv）以及pressure level（高度层）上的U、V、Z（geopotential）位势。下载的网址如下：

高度层：ERA5 hourly data on pressure levels from 1940 to present (copernicus.eu)

地表：ERA5 hourly data on single levels from 1940 to present (copernicus.eu) 

下载好之后就能进行使用。 

数据可视化

由于metpy的计算效率高，本文使用其对涡度进行计算，自编写的计算函数在上一篇当中有提到。

本文以2024年8月16 UTC 00:00 时刻的ERA5数据进行台风路径判断相关初始场的可视化。

import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import metpy.calc as mpcalc
#选择日期
date_sel='2024-08-16'
# 设置绘图范围
extent = [99, 145.1, 0, 45.1]
# 读取数据
df = xr.open_dataset(r"D:\Pangu-Weather-ReadyToGo\Test-Data\Z-Q-T-U-V-13levs-2024.08.00.nc")
# 导入第二个文件
df1 = xr.open_dataset(r"D:\Pangu-Weather-ReadyToGo\Test-Data\mslp-10mU-10mV-2mT-2024.08.00.nc")
u = df.u.loc[date_sel, 850, 48:-1, 99:148]  # 850hpa
v = df.v.loc[date_sel, 850, 48:-1, 99:148]  # m/s
lon, lat = u.longitude, u.latitude
#计算涡度
vorticity = mpcalc.vorticity(u, v)
vorticity = vorticity * 1e5
#计算厚度
z850 = df.z.loc[date_sel, 850, 48:-1, 99:148]  # 850 hPa 位势高度
z200 = df.z.loc[date_sel, 200, 48:-1, 99:148]  # 200 hPa 位势高度
thickness = z200 - z850
# 计算风速
u10 = df1.u10.loc[date_sel, 48:-1, 99:148]
v10 = df1.v10.loc[date_sel, 48:-1, 99:148]
ws10 = np.sqrt(u10**2 + v10**2)
#海平面气压
slp = df1.msl.loc[date_sel, 48:-1, 99:148]
slp = slp / 100  # hpa
# ---------创建画布----------
fig = plt.figure(figsize=(10, 8))
# 左上子图：绘制ws10
ax1 = fig.add_axes([0.05, 0.55, 0.4, 0.4], projection=ccrs.PlateCarree())
ax1.set_extent(extent, crs=ccrs.PlateCarree())
ax1.coastlines('50m', linewidth=1.2)
ax1.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', edgecolor='face', facecolor='lightgrey'))
gl = ax1.gridlines(draw_labels=True, lw=1, color='k', alpha=0, ls='--')
gl.top_labels = False
gl.right_labels = False
ax1.set_title('10m Wind', fontsize=10, pad=5, loc='center')
ax1.set_title(date_sel, fontsize=8, pad=5, loc='right')
cair1 = ax1.contourf(lon, lat, ws10, extend='both', levels=np.linspace(0, 25, 26),
                     cmap='viridis', transform=ccrs.PlateCarree())
cbar1 = fig.colorbar(cair1, ax=ax1, orientation="vertical", label='m/s', shrink=0.65)
# 右上子图：绘制vorticity
ax2 = fig.add_axes([0.5, 0.55, 0.4, 0.4], projection=ccrs.PlateCarree())
ax2.set_extent(extent, crs=ccrs.PlateCarree())
ax2.coastlines('50m', linewidth=1.2)
ax2.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', edgecolor='face', facecolor='lightgrey'))
gl = ax2.gridlines(draw_labels=True, lw=1, color='k', alpha=0, ls='--')
gl.top_labels = False
gl.right_labels = False
ax2.set_title('850hPa Vorticity', fontsize=10, pad=5, loc='center')
ax2.set_title(date_sel, fontsize=8, pad=5, loc='right')
cair2 = ax2.contourf(lon, lat, vorticity, extend='both', levels=np.linspace(-20, 60, 81),
                     cmap='jet', transform=ccrs.PlateCarree())
cbar2 = fig.colorbar(cair2, ax=ax2, orientation="vertical",
                     label='10$^-$$^5$ s$^-$$^1$', shrink=0.65)
# 左下子图：绘制厚度
ax3 = fig.add_axes([0.05, 0.08, 0.4, 0.4], projection=ccrs.PlateCarree())
ax3.set_extent(extent, crs=ccrs.PlateCarree())
ax3.coastlines('50m', linewidth=1.2)
ax3.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', edgecolor='face', facecolor='lightgrey'))
gl = ax3.gridlines(draw_labels=True, lw=1, color='k', alpha=0, ls='--')
gl.top_labels = False
gl.right_labels = False
ax3.set_title('Thickness 850-200hPa', fontsize=10, pad=5, loc='center')
ax3.set_title(date_sel, fontsize=8, pad=5, loc='right')
cair3 = ax3.contourf(lon, lat, thickness, 60, extend='both',
                     cmap='jet', transform=ccrs.PlateCarree())
cbar3 = fig.colorbar(cair3, ax=ax3, orientation="vertical", label='m', shrink=0.65)
# 右下子图：绘制海平面气压
ax4 = fig.add_axes([0.5, 0.08, 0.4, 0.4], projection=ccrs.PlateCarree())
ax4.set_extent(extent, crs=ccrs.PlateCarree())
ax4.coastlines('50m', linewidth=1.2)
ax4.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', edgecolor='face', facecolor='lightgrey'))
gl = ax4.gridlines(draw_labels=True, lw=1, color='k', alpha=0, ls='--')
gl.top_labels = False
gl.right_labels = False
ax4.set_title('Mean SLP', fontsize=10, pad=5, loc='center')
ax4.set_title(date_sel, fontsize=8, pad=5, loc='right')
cair4 = ax4.contourf(lon, lat, slp, extend='both', levels=np.arange(990, 1021, 1),
                     cmap='jet', transform=ccrs.PlateCarree())
cbar4 = fig.colorbar(cair4, ax=ax4, orientation="vertical",
                     label='hPa', shrink=0.65)
# 展示
plt.show()

结果展示

四张图中在140°E，30°N左右位置上的四个明显圆圈即为2024年07号台风安比，此时是强台风级别。