Python计算涡度、散度

本文提供了两种涡度和散度的计算方法，一个是自定义的计算函数，一个是导入metpy库进行计算，最后对结果进行可视化。

涡度、散度概念

涡度表示流体的旋转程度，在北半球，气旋的气流逆时针旋转，涡度为正。南半球，气旋的气流顺时针旋转，涡度为负。在二维上，涡度表示为v风在x方向上的变率减去u风在y方向上的变率，即水平速度场的旋转率。

散度描述流体的收缩或扩张，是流体团在xyz三个方向上的法形变率之和，又称为体形变率，当散度为0时流体为不可压缩流体。散度在二维上表示为：

在实际的二维平面计算当中，由于地球自转，需要考虑行星涡度项，即2倍的地球自转角速度。

涡度：

散度： 

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
def div_vor_cal(u_field,v_field,pre):#输入uv和数据精度
    vor = np.empty([u_field.shape[0],u_field.shape[1]])#涡度
    div = np.empty([v_field.shape[0],v_field.shape[1]])#散度
    # 常数和计算参数
    R = 6371000  # 地球半径，单位：米
    Len = 2.0 * np.pi * R / 360.0#周长
    dy = pre * Len#经度方向格点大小
    lat_shape=u_field.shape[0]-1
    lon_shape = u_field.shape[1]-1
    lat=u_field.latitude
    lat_ct_start=lat[1]
    total_sum=0
    for i in range(1, lat_shape):
        for j in range(1,lon_shape):
            # ------------中间拆分------------
            # 考虑经纬网格大小随纬度变化
            fi = np.radians(lat_ct_start - i * pre)#弧度制高度角 or (lat_ct_start - i * pre) * 2 * np.pi / 360
            dx =  dy * np.cos(fi)
            #计算涡度
            vvy = v_field[i, j+1] - v_field[i, j-1]
            uux = u_field[i-1, j] - u_field[i+1, j]
            vor[i, j] = 0.5*(( vvy / dx) - (uux / dy)) + (u_field[i, j] / R) * np.tan(fi)#行星涡度项
            #计算散度
            vvx = v_field[i-1, j] - v_field[i+1, j]
            uuy = u_field[i, j+1] - u_field[i, j-1]
            div[i, j] = 0.5 * ( (vvx / dy) +  (uuy / dx) )- (v_field[i, j] / R) * np.tan(fi)
            total_sum+=1
        print("计算进度：" ,np.round(total_sum*100/(lat_shape*lon_shape),1),"%")

    return vor, div

由于自定义循环的计算效率问题，以上只考虑了中央拆分方法计算涡度和散度，若需要更精细的拆分完整的计算函数将在文末给出。

涡度、散度绘图

接下来就是绘图部分。这里使用了ERA5的逐日UTC 00:00时的再分析数据，精度为0.25*0.25，挑选了2024年8月17日850hpa上的数据进行可视化

df=xr.open_dataset(r"D:\Pangu-Weather-ReadyToGo\Test-Data\Z-Q-T-U-V-13levs-2024.08.00.nc")
date_sel='2024-08-17'
#(valid_time, pressure_level, latitude, longitude)
u=df.u.loc[date_sel,850,48:-1,99:150]#850hpa
v=df.v.loc[date_sel,850,48:-1,99:150]#m/s
vorticity,divergence=div_vor_cal(u,v,0.25)
lon,lat=u.longitude,u.latitude

# 为了更好地表示涡度，可以乘以一个比例因子（如1e5）
vorticity = vorticity * 1e5
divergence=divergence*1e5
# 创建画布和子图
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 5), subplot_kw={'projection': ccrs.PlateCarree()})
# 设置绘图范围
extent = [99, 145.1, 0, 45.1]
# 左侧绘制ws10
ax1 = axes[0]
ax1.set_extent(extent, crs=ccrs.PlateCarree())
ax1.coastlines('50m', linewidth=1.2)
ax1.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', edgecolor='k', facecolor='None'),zorder=2,alpha=0.5)
gl = ax1.gridlines(draw_labels=True, lw=1, color='k', alpha=0, ls='--')
gl.top_labels = False
gl.right_labels = False
ax1.set_title('850hPa Divergence', fontsize=12, pad=8, loc='center')
ax1.set_title(date_sel, fontsize=8, pad=8, loc='right')
cair1 = ax1.contourf(lon, lat, divergence, extend='both', levels=np.linspace(-10, 30, 41),
                     cmap='jet', transform=ccrs.PlateCarree())
cbar1 = fig.colorbar(cair1, ax=ax1, orientation="vertical", label='10$^-$$^5$ s$^-$$^1$',shrink=0.65)

# 右侧绘制vorticity
ax2 = axes[1]
ax2.set_extent(extent, crs=ccrs.PlateCarree())
ax2.coastlines('50m', linewidth=1.2)
ax2.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', edgecolor='k', facecolor='None'),zorder=2,alpha=0.5)
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=12, pad=8, loc='center')
ax2.set_title(date_sel, fontsize=8, pad=8, 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)
#调整间距，绘图
plt.tight_layout()
plt.show()

结果如下：

 对于涡度图可以很明显的看到日本东京附近有一个很明显的涡度极大值中心，这个就是2024年07号台风安比，在此时的中心风力达到15级，中心气压940hpa，最大风速48m/s，为强台风级别。

使用metpy库

对于自定义的涡度散度计算，由于python自身的计算效率问题以及循环、数据精度过高等原因，计算一次的时间成本过高，故采用已经由开发者封装好的metpy库来计算，时间成本将大大缩短。

代码如下：

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
df=xr.open_dataset(r"D:\Pangu-Weather-ReadyToGo\Test-Data\Z-Q-T-U-V-13levs-2024.08.00.nc")
#(valid_time, pressure_level, latitude, longitude)
u=df.u.loc['2024-08-16',850,48:-1,99:148]#850hpa
v=df.v.loc['2024-08-16',850,48:-1,99:148]#m/s
lon,lat=u.longitude,u.latitude
vorticity=mpcalc.vorticity(u,v)
vorticity=vorticity*1e5
divergence=mpcalc.divergence(u,v)
divergence=divergence*1e5
# 创建画布和子图
fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(10, 5), subplot_kw={'projection': ccrs.PlateCarree()})
# 设置绘图范围
extent = [99, 145.1, 0, 45.1]
# 左侧绘制divergence
ax1 = axes[0]
ax1.set_extent(extent, crs=ccrs.PlateCarree())
ax1.coastlines('50m', linewidth=1.2)
ax1.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', edgecolor='k', facecolor='None'),zorder=2,alpha=0.5)
gl = ax1.gridlines(draw_labels=True, lw=1, color='k', alpha=0, ls='--')
gl.top_labels = False
gl.right_labels = False
ax1.set_title('850hPa Divergence', fontsize=12, pad=8, loc='center')
ax1.set_title('2024-08-17', fontsize=8, pad=8, loc='right')
cair1 = ax1.contourf(lon, lat, divergence, extend='both', levels=np.linspace(-10, 30, 41),
                     cmap='jet', transform=ccrs.PlateCarree())
cbar1 = fig.colorbar(cair1, ax=ax1, orientation="vertical", label='10$^-$$^5$ s$^-$$^1$',shrink=0.65)

# 右侧绘制vorticity
ax2 = axes[1]
ax2.set_extent(extent, crs=ccrs.PlateCarree())
ax2.coastlines('50m', linewidth=1.2)
ax2.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', edgecolor='k', facecolor='None'),zorder=2,alpha=0.5)
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=12, pad=8, loc='center')
ax2.set_title('2024-08-17', fontsize=8, pad=8, 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)
#调整间距，绘图
plt.tight_layout()
plt.show()

以上使用metpy的计算结果与自定义函数的结果差别很小，故不展示绘图结果。

自定义函数的补充

对于自定义的计算，在之前还未考虑经纬网格边缘的计算，这里将使用四角拆分和四边拆分对自定义函数进行补充。代码如下：

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
def div_vor_cal(u_field,v_field,pre):#输入uv和数据精度
    vor = np.empty([u_field.shape[0],u_field.shape[1]])#涡度
    div = np.empty([v_field.shape[0],v_field.shape[1]])#散度
    # 常数和计算参数
    R = 6371000  # 地球半径，单位：米
    Len = 2.0 * np.pi * R / 360.0#周长
    dy = pre * Len#经度方向格点大小
    lat_shape=u_field.shape[0]-1
    lon_shape = u_field.shape[1]-1
    lat=u_field.latitude
    lat_ct_start=lat[1]
    total_sum=0
    for i in range(1, lat_shape):
        for j in range(1,lon_shape):
            # ------------中间拆分------------
            # 考虑经纬网格大小随纬度变化
            fi = np.radians(lat_ct_start - i * pre)#弧度制高度角 or (lat_ct_start - i * pre) * 2 * np.pi / 360
            dx =  dy * np.cos(fi)
            #计算涡度
            vvy = v_field[i, j+1] - v_field[i, j-1]
            uux = u_field[i-1, j] - u_field[i+1, j]
            vor[i, j] = 0.5*(( vvy / dx) - (uux / dy)) + (u_field[i, j] / R) * np.tan(fi)#行星涡度项
            #计算散度
            vvx = v_field[i-1, j] - v_field[i+1, j]
            uuy = u_field[i, j+1] - u_field[i, j-1]
            div[i, j] = 0.5 * ( (vvx / dy) +  (uuy / dx) )- (v_field[i, j] / R) * np.tan(fi)
            #考虑边缘的更细节的拆分
            # ------------四边拆分------------
            #计算涡度
            vv_edge1 = v_field[0, j+1] - v_field[0, j-1]
            uu_edge1  = u_field[0, j] - u_field[1, j]
            vor[0, j] = 0.5*(( vv_edge1  / dx) - (uu_edge1 / dy)) + (u_field[i, j] / R) * np.tan(fi)

            vv_edge2 = v_field[i, 1] - v_field[i, 0]
            uu_edge2  = u_field[i-1, 0] - u_field[i+1, 0]
            vor[i, 0] = 0.5*(( vv_edge2  / dx) - (uu_edge2 / dy)) + (u_field[i, j] / R) * np.tan(fi)

            vv_edge3 = v_field[i, lon_shape] - v_field[i, lon_shape-1]
            uu_edge3  = u_field[i-1, lon_shape] - u_field[i+1, lon_shape]
            vor[i, lon_shape] = 0.5*(( vv_edge3  / dx) - (uu_edge3 / dy)) + (u_field[i, j] / R) * np.tan(fi)

            vv_edge4 = v_field[lat_shape, j+1] - v_field[lat_shape, j-1]
            uu_edge4  = u_field[lat_shape-1, j] - u_field[lat_shape, j]
            vor[lat_shape, j] = 0.5*(( vv_edge4  / dx) - (uu_edge4 / dy)) + (u_field[i, j] / R) * np.tan(fi)
            #计算散度
            vvv_edge1 = v_field[0, j] - v_field[1, j]
            uuu_edge1 = u_field[0, j + 1] - u_field[0, j - 1]
            div[0, j] = 0.5 * ((vvv_edge1 / dy) + (uuu_edge1 / dx)) - (v_field[0, j] / R) * np.tan(fi)

            vvv_edge2 = v_field[i - 1, 0] - v_field[i + 1, 0]
            uuu_edge2 = u_field[i, 1] - u_field[i, 0]
            div[i, 0] = 0.5 * ((vvv_edge2 / dy) + (uuu_edge2 / dx)) - (v_field[i, 0] / R) * np.tan(fi)

            vvv_edge3 = v_field[i - 1, lon_shape] - v_field[i + 1, lon_shape]
            uuu_edge3 = u_field[i, lon_shape] - u_field[i, lon_shape - 1]
            div[i, lon_shape] = 0.5 * ((vvv_edge3 / dy) + (uuu_edge3 / dx)) - (v_field[i, lon_shape] / R) * np.tan(fi)

            vvv_edge4 = v_field[lat_shape - 1, j] - v_field[lat_shape, j]
            uuu_edge4 = u_field[lat_shape, j + 1] - u_field[lat_shape, j - 1]
            div[lat_shape, j] = 0.5 * ((vvv_edge4 / dy) + (uuu_edge4 / dx)) - (v_field[lat_shape, j] / R) * np.tan(fi)
            total_sum+=1
        print("计算进度：" ,np.round(total_sum*100/(lat_shape*lon_shape),1),"%")
    # 四个角点的散度和涡度计算
    fi1=np.radians(lat[0])
    dx1=dy * np.cos(fi1)
    fi2=np.radians(lat[len(lat)])
    dx2=dy * np.cos(fi2)
    # ------------左上角 (0, 0) ------------
    vvv_corner1 = v_field[0, 0] - v_field[1, 0]
    uuu_corner1 = u_field[0, 1] - u_field[0, 0]
    div[0, 0] = 0.5 * ((vvv_corner1 / dy) + (uuu_corner1 / dx1)) - (v_field[0, 0] / R) * np.tan(fi1)

    vv_corner1 = v_field[0, 1] - v_field[0, 0]
    uu_corner1 = u_field[0, 0] - u_field[1, 0]
    vor[0, 0] = 0.5 * ((vv_corner1 / dx1) - (uu_corner1 / dy)) + (u_field[0, 0] / R) * np.tan(fi1)
    # ------------右上角 (0, lon_shape) ------------
    vvv_corner2 = v_field[0, lon_shape] - v_field[1, lon_shape]
    uuu_corner2 = u_field[0, lon_shape] - u_field[0, lon_shape - 1]
    div[0, lon_shape] = 0.5 * ((vvv_corner2 / dy) + (uuu_corner2 / dx1)) - (v_field[0, lon_shape] / R) * np.tan(fi1)

    vv_corner2 = v_field[0, lon_shape] - v_field[0, lon_shape - 1]
    uu_corner2 = u_field[0, lon_shape] - u_field[1, lon_shape]
    vor[0, lon_shape] = 0.5 * ((vv_corner2 / dx1) - (uu_corner2 / dy)) + (u_field[0, lon_shape] / R) * np.tan(fi1)
    # ------------左下角 (lat_shape, 0) ------------
    vvv_corner3 = v_field[lat_shape, 0] - v_field[lat_shape - 1, 0]
    uuu_corner3 = u_field[lat_shape, 1] - u_field[lat_shape, 0]
    div[lat_shape, 0] = 0.5 * ((vvv_corner3 / dy) + (uuu_corner3 / dx2)) - (v_field[lat_shape, 0] / R) * np.tan(fi2)

    vv_corner3 = v_field[lat_shape, 1] - v_field[lat_shape, 0]
    uu_corner3 = u_field[lat_shape - 1, 0] - u_field[lat_shape, 0]
    vor[lat_shape, 0] = 0.5 * ((vv_corner3 / dx2) - (uu_corner3 / dy)) + (u_field[lat_shape, 0] / R) * np.tan(fi2)
    # ------------右下角 (lat_shape, lon_shape) ------------
    vvv_corner4 = v_field[lat_shape, lon_shape] - v_field[lat_shape - 1, lon_shape]
    uuu_corner4 = u_field[lat_shape, lon_shape] - u_field[lat_shape, lon_shape - 1]
    div[lat_shape, lon_shape] = 0.5 * ((vvv_corner4 / dy) + (uuu_corner4 / dx2)) - (
                v_field[lat_shape, lon_shape] / R) * np.tan(fi2)
    # 涡度
    vv_corner4 = v_field[lat_shape, lon_shape] - v_field[lat_shape, lon_shape - 1]
    uu_corner4 = u_field[lat_shape - 1, lon_shape] - u_field[lat_shape, lon_shape]
    vor[lat_shape, lon_shape] = 0.5 * ((vv_corner4 / dx2) - (uu_corner4 / dy)) + (
                u_field[lat_shape, lon_shape] / R) * np.tan(fi2)

    return vor, div