气象编程 | Python绘制全球涡度、散度和垂直速度

 计算部分
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs

# 读取netCDF文件
data1 = xr.open_dataset('uwnd.20200812.nc')
data2 = xr.open_dataset('vwnd.20200812.nc')
data3 = xr.open_dataset('omega.20200812.1000hPa.nc')

LAT = data1['lat'].values
LON = data1['lon'].values

ORGANnC1 = data1['uwnd']
ORGANnC2 = data2['vwnd']
ORGANnC3 = data3['omega']

# 提取 1000 hpa 层级的数据
u1000 = ORGANnC1.sel(time="2020-08-12",level='1e+03').values
v1000 = ORGANnC2.sel(time="2020-08-12",level='1e+03').values
# 提取 850 hpa 层级的数据
u850 = ORGANnC1.sel(time="2020-08-12",level='850').values
v850 = ORGANnC2.sel(time="2020-08-12",level='850').values


#提取100 hpa 层级的垂直速度
omega = ORGANnC3.sel(time="2020-08-12",level='1e+03').values

'''
散度dv,涡度rv; 
'''

# 计算涡度和散度
def calculate_vorticity_divergence(u_field, v_field):

 计算涡度和散度。

    参数:
    u_field -- U风分量数组
    v_field -- V风分量数组
    a -- 地球半径
    dL -- 网格间距

    返回:
    rv -- 涡度数组
    dv -- 散度数组
    """
    rv = np.zeros_like(u_field)
    dv = np.zeros_like(v_field)
    # 常数和计算参数
    a = 6371110.0
    pi = 3.1415926
    L = 2.0 * pi * a / 360.0
    dL = 2.5 * L
    for m in range(1, u_field.shape[0] - 1):
        for n in range(1, u_field.shape[1] - 1):
            vv = v_field[m, n+1] - v_field[m, n-1]
            uu = u_field[m-1, n] - u_field[m+1, n]
            fi = (60 - m * 2.5) * 2 * np.pi / 360
            dx = dL * np.cos(fi)
            rv[m, n] = 0.5 * vv / dx - 0.5 * uu / dL + (u_field[m, n] / a) * np.tan(fi)

            vvv = v_field[m-1, n] - v_field[m+1, n]
            uuu = u_field[m, n+1] - u_field[m, n-1]
            dv[m, n] = 0.5 * vvv / dL + 0.5 * uuu / dx - (v_field[m, n] / a) * np.tan(fi)

    return rv, dv   





    

#调用函数,计算涡度、散度
rv1000, dv1000 = calculate_vorticity_divergence(u1000, v1000)
rv850, dv850 = calculate_vorticity_divergence(u850, v850)

# 计算不同气压层的垂直速度
omega850 = omega + (dv1000 + dv850) * 150 * 0.5 * 100

绘图部分
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import cartopy.mpl.ticker as cticker
import cmaps
import geopandas as gpd
import matplotlib.colors as colors
from matplotlib.path import Path
from cartopy.mpl.patch import geos_to_path
import matplotlib as mpl
from matplotlib.colors import TwoSlopeNorm
from matplotlib.ticker import ScalarFormatter

#中文
plt.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
plt.rcParams['axes.unicode_minus'] = False    # 修复保存图像是负号'-'显示为方块的问题

shp = gpd.read_file("china.shp", encoding='utf-8')["geometry"]

# 设定坐标轴及标题的函数
def setup_axes(ax, extent=[20, 160, -30, 60]):
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    ax.set_xticks(range(30, 161, 30), crs=ccrs.PlateCarree())
    ax.xaxis.set_major_formatter(cticker.LongitudeFormatter())
    ax.set_yticks(range(-30, 61, 30), crs=ccrs.PlateCarree())
    ax.yaxis.set_major_formatter(cticker.LatitudeFormatter())
    ax.tick_params(axis='x')
    ax.tick_params(axis='y')
    ax.set_xlabel('')
    ax.set_ylabel('')

# 计算最小值和最大值
a = rv850
min_cvalue = -a.max()
max_cvalue = a.max()
# 生成14个等间隔的点（产生12个间隔）
cvalues = np.linspace(min_cvalue, max_cvalue, 13)
# 寻找最接近0的两个cvalues值
zero_index = np.searchsorted(cvalues, 0)
lower_bound = cvalues[max(0, zero_index-1)]
upper_bound = cvalues[min(len(cvalues)-1, zero_index)]

# 绘图
fig = plt.figure(figsize=(20, 9))
gs = gridspec.GridSpec(1, 1)
ax = plt.subplot(gs[0, 0], projection=ccrs.PlateCarree())

# 绘制
contourf_obj = ax.contourf(LON, LAT, a,cmap=cmaps.BlueRed, levels=cvalues , transform=ccrs.PlateCarree())
setup_axes(ax)

#添加海岸线和中国地图
ax.coastlines()
ax.add_geometries(shp, crs=ccrs.PlateCarree(), facecolor="none")

ax.set_title('850hPa 涡度', fontsize=20)

# 添加颜色条，传入可绘制对象作为参数
cbar_ax = fig.add_axes([0.21, 0.03, 0.6, 0.05])
cbar = plt.colorbar(contourf_obj, cax=cbar_ax, orientation='horizontal')
# 创建并设置ScalarFormatter
formatter = ScalarFormatter(useMathText=True)
formatter.set_powerlimits((-1, 1))  # 设置科学记数法的范围，根据需要调整
cbar.ax.yaxis.set_major_formatter(formatter)  # 应用到颜色条的y轴

# 保存图片并显示
plt.savefig('G:/dataarray/tianzhen/涡度850.jpg', bbox_inches='tight', dpi=1000)
plt.show()