可视化 | 台风方位角平均的半径-气压剖面图

from scipy import interpolate
import metpy.calc as mpcalc  
from metpy.units import units  
import xarray as xr   
import numpy as np 
import matplotlib.pyplot as plt
import matplotlib as mpl
# mpl.rcParams['font.sans-serif'] = ['Times New Roman']#设置默认字体
mpl.rcParams['font.size'] = '14' # 设置字体大小

Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.14.1

f = '/home/mw/input/nc_sample3575/data_example.nc'
#CMA给出的台风中心
lon_t=128.9
lat_t=20.0

#nc数据读取
ds = xr.open_dataset(f)
lat = ds.latitude
lon = ds.longitude
levels = ds.level

 #azimuths是极坐标系中的角度，ranges是半径，可以根据自己需要设置
azimuths = np.linspace(0,360,73)*units.degree
ranges = np.linspace(0,1000,101)*1000*units.meter 
vt_am = np.zeros((len(levels),len(ranges)))
vr_am = np.zeros((len(levels),len(ranges)))

#利用850hPa的相对涡度寻找再分析资料中的台风中心
lon_range = lon[(lon>=lon_t-2) & (lon<=lon_t+2)]
lat_range = lat[(lat>=lat_t-2) & (lat<+lat_t+2)]
vo = ds.vo.sel(level=850, longitude=lon_range, latitude=lat_range)

for m in lat_range:
    for n in lon_range:
        if vo.sel(longitude=n, latitude=m) == np.max(vo):
            lat_v = m
            lon_v = n

#利用metpy库可以十分便捷的得到插值后的经纬度坐标
lon_a,lat_a = mpcalc.azimuth_range_to_lat_lon(azimuths,ranges,lon_v,lat_v)

#因为ERA5的数据分辨率是0.25°，为了保证插值后不产生NAN，边界上各扩大一个格点
lons = lon[(lon>=lon_a.min()-0.25) & (lon<=lon_a.max()+0.25)] 
lats = lat[(lat>=lat_a.min()-0.25) & (lat<=lat_a.max()+0.25)]
lon_s,lat_s = np.meshgrid(lons,lats)

#构造插值前后的格点矩阵
grid_in = np.concatenate([lon_s.reshape(-1,1), lat_s.reshape(-1,1)], axis=1)
grid_out = np.concatenate([lon_a.reshape(-1,1), lat_a.reshape(-1,1)], axis=1)

for j in range(0,len(levels)):
    u_in = ds['u'].sel(level= levels[j], longitude=lons, latitude=lats)
    v_in = ds['v'].sel(level= levels[j], longitude=lons, latitude=lats)
    u_out = interpolate.griddata(grid_in, np.array(u_in).flatten(), grid_out, method='cubic')
    v_out = interpolate.griddata(grid_in, np.array(v_in).flatten(), grid_out, method='cubic')
    u_out = u_out.reshape((len(azimuths),len(ranges)))
    v_out = v_out.reshape((len(azimuths),len(ranges)))

    #计算切向风、径向风
    vt = np.zeros((len(azimuths),len(ranges)))
    vr = np.zeros((len(azimuths),len(ranges)))

    for k in range(0,len(azimuths)):
        vt[k,:] = v_out[k,:]*np.cos(azimuths[k]*np.pi/180)-u_out[k,:]*np.sin(azimuths[k]*np.pi/180)
        vr[k,:] = u_out[k,:]*np.cos(azimuths[k]*np.pi/180)+v_out[k,:]*np.sin(azimuths[k]*np.pi/180)
    #计算方位角平均
    vt_am[j,:] = np.mean(vt,axis=0)
    vr_am[j,:] = np.mean(vr,axis=0)

#为了图形美观，做2次9点平滑
vt_am = mpcalc.smooth_n_point(vt_am,9,2) 
vr_am = mpcalc.smooth_n_point(vr_am,9,2) 

#画图
plt.figure(1, figsize=(13., 5.))

ax1 = plt.subplot(121) 
ax1.invert_yaxis() 
ax1.set_yscale('symlog')
ax1.set_yticks([1000,850,700,500,400,300,200,150,100,50])
ax1.set_yticklabels([1000,850,700,500,400,300,200,150,100,50])        
ax1.set_xticks(range(0,1200000,200000))
ax1.set_xticklabels(range(0,12,2))
ax1.set_ylabel('Pressure (hPa)')
ax1.set_xlabel('Radius (100km)')
ax1.set_title('tangential wind')
fig1 = ax1.contourf(ranges, levels, vt_am, 
                    np.arange(-10,12,2),
                    cmap='bwr',
                    extend='both')
plt.colorbar(fig1,orientation='vertical',shrink=0.75)

ax2 = plt.subplot(122) 
ax2.invert_yaxis() 
ax2.set_yscale('symlog')
ax2.set_yticks([1000,850,700,500,400,300,200,150,100,50])
ax2.set_yticklabels([1000,850,700,500,400,300,200,150,100,50])        
ax2.set_xticks(range(0,1200000,200000))
ax2.set_xticklabels(range(0,12,2))
ax2.set_ylabel('Pressure (hPa)')
ax2.set_xlabel('Radius (100km)')
ax2.set_title('radial wind')
fig2 = ax2.contourf(ranges, levels, vr_am, 
                    np.arange(-10,12,2),
                    cmap='bwr',
                    extend='both')
plt.colorbar(fig2,orientation='vertical',shrink=0.75)
    
plt.show()