【气象常用】wrf后处理-风场

效果图：

主要步骤：

1. 数据准备：wrf输出文件

2. 数据处理：获取wrf中500hPa风场

3. 图像绘制：绘制风场

详细代码：着急的直接拖到最后有完整代码

步骤一：导入库包及图片存储路径并设置中文字体为宋体，西文为新罗马（没有的库包要先下好奥）

###############################################################################
# 导入库包并设置字体
import numpy as np
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from scipy.interpolate import RegularGridInterpolator  
from netCDF4 import Dataset
import os
from cartopy.io.shapereader import Reader
from wrf import getvar, interplevel
from wrf import (getvar, interplevel, to_np, latlon_coords, get_cartopy,smooth2d,
                 cartopy_xlim, cartopy_ylim)
# 设置西文字体为新罗马字体，中文为宋体，字号为12
from matplotlib import rcParams
config = {
            "font.family": 'serif',
            "font.size": 12,
            "mathtext.fontset": 'stix',
            "font.serif": ['SimSun'],
          }
rcParams.update(config)
rcParams['axes.unicode_minus']=False

datapath = r'H:/00.csdn/01data/'
figpath = r'H:/00.csdn/02fig/'
shppath = r'H:/00.csdn/04shp/cn_shp/Province_9/Province_9.shp'

步骤二：读入数据并获取500hPa风场

###############################################################################
# 读入数据
f = Dataset(datapath + "wrfout_d02_2023-06-01_18_00_00.nc")

p = getvar(f, "pressure")
u = getvar(f, "ua", units="kt")
v = getvar(f, "va", units="kt")
lats, lons = np.array(to_np(latlon_coords(p)))

# 对wrf输出数据进行处理，提取500hPa风场
u_500 = interplevel(u, p, 500)
v_500 = interplevel(v, p, 500)

步骤三：绘制风场图像

###############################################################################
# 画图
fig = plt.figure(figsize=(15, 8))
def plot_cont(u_data, v_data, location, extent, show, title, titleposition):
    ax = fig.add_axes(location, projection=ccrs.PlateCarree())

    cq = ax.quiver(lons[::5, ::5], lats[::5, ::5],
                   u_data[::5, ::5], v_data[::5, ::5],
                   color='k',
                   scale=300, zorder=3, width=0.004,
                   headwidth=3, headlength=4.5, transform=ccrs.PlateCarree())
    
    # 添加quiverkey作为参考图例  
    qk = ax.quiverkey(cq, 119.3, 31.2, 20, '20 m/s',
                      labelpos='E', color='k', 
                      coordinates='data', 
                      zorder=200)  
    #添加中国国界
    ax.add_geometries(Reader(shppath).geometries(), ccrs.PlateCarree(),
                      facecolor='none', edgecolor='k', linewidth=0.7)  
    
    # 设置显示范围
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    
    # 设置是否显示经纬网格线及坐标刻度
    if show == True:
        # 画经纬度网格
        gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=1.2, 
                          color='k', alpha=0.5, linestyle='--')
        
        title = ax.set_title(title, fontsize=12)
        title.set_position(titleposition)  # 调整位置
        
        gl.top_labels = False  # 关闭顶端的经纬度标签
        gl.right_labels = False  # 关闭右侧的经纬度标签
        gl.xformatter = LONGITUDE_FORMATTER  # x轴设为经度的格式
        gl.yformatter = LATITUDE_FORMATTER  # y轴设为纬度的格式
    
        # 设置经纬度网格的间隔及坐标标签
        gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+1, 2))
        gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+1, 2))
        ax.set_xticks(list(range(extent[0], extent[1]+1, 2)), crs=ccrs.PlateCarree())
        ax.set_yticks(list(range(extent[2], extent[3]+1, 2)), crs=ccrs.PlateCarree())
        
        lon_formatter = LongitudeFormatter(zero_direction_label=True)  # zero_direction_label用来设置经度的0度加不加E和W
        lat_formatter = LatitudeFormatter()
        ax.xaxis.set_major_formatter(lon_formatter)
        ax.yaxis.set_major_formatter(lat_formatter)

cc1 = plot_cont(u_500, v_500, 
                [0.1, 0.1, 0.8, 0.8], [110, 120, 23, 31], 
                show=True, title='(a)500hPa', titleposition=[0.04, 1.05])

步骤四：存图

# ###############################################################################
# # 存图
plt.savefig(figpath+'213 wrf风场', dpi=600, bbox_inches = 'tight')
plt.show()

完整代码在这里：

###############################################################################
# 导入库包并设置字体
import numpy as np
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
from scipy.interpolate import RegularGridInterpolator  
from netCDF4 import Dataset
import os
from cartopy.io.shapereader import Reader
from wrf import getvar, interplevel
from wrf import (getvar, interplevel, to_np, latlon_coords, get_cartopy,smooth2d,
                 cartopy_xlim, cartopy_ylim)
# 设置西文字体为新罗马字体，中文为宋体，字号为12
from matplotlib import rcParams
config = {
            "font.family": 'serif',
            "font.size": 12,
            "mathtext.fontset": 'stix',
            "font.serif": ['SimSun'],
          }
rcParams.update(config)
rcParams['axes.unicode_minus']=False

datapath = r'H:/00.csdn/01data/'
figpath = r'H:/00.csdn/02fig/'
shppath = r'H:/00.csdn/04shp/cn_shp/Province_9/Province_9.shp'
###############################################################################
# 读入数据
f = Dataset(datapath + "wrfout_d02_2023-06-01_18_00_00.nc")

p = getvar(f, "pressure")
u = getvar(f, "ua", units="kt")
v = getvar(f, "va", units="kt")
lats, lons = np.array(to_np(latlon_coords(p)))

# 对wrf输出数据进行处理，提取500hPa风场
u_500 = interplevel(u, p, 500)
v_500 = interplevel(v, p, 500)
###############################################################################
# 画图
fig = plt.figure(figsize=(15, 8))
def plot_cont(u_data, v_data, location, extent, show, title, titleposition):
    ax = fig.add_axes(location, projection=ccrs.PlateCarree())

    cq = ax.quiver(lons[::5, ::5], lats[::5, ::5],
                   u_data[::5, ::5], v_data[::5, ::5],
                   color='k',
                   scale=300, zorder=3, width=0.004,
                   headwidth=3, headlength=4.5, transform=ccrs.PlateCarree())
    
    # 添加quiverkey作为参考图例  
    qk = ax.quiverkey(cq, 119.3, 31.2, 20, '20 m/s',
                      labelpos='E', color='k', 
                      coordinates='data', 
                      zorder=200)  
    #添加中国国界
    ax.add_geometries(Reader(shppath).geometries(), ccrs.PlateCarree(),
                      facecolor='none', edgecolor='k', linewidth=0.7)  
    
    # 设置显示范围
    ax.set_extent(extent, crs=ccrs.PlateCarree())
    
    # 设置是否显示经纬网格线及坐标刻度
    if show == True:
        # 画经纬度网格
        gl = ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=1.2, 
                          color='k', alpha=0.5, linestyle='--')
        
        title = ax.set_title(title, fontsize=12)
        title.set_position(titleposition)  # 调整位置
        
        gl.top_labels = False  # 关闭顶端的经纬度标签
        gl.right_labels = False  # 关闭右侧的经纬度标签
        gl.xformatter = LONGITUDE_FORMATTER  # x轴设为经度的格式
        gl.yformatter = LATITUDE_FORMATTER  # y轴设为纬度的格式
    
        # 设置经纬度网格的间隔及坐标标签
        gl.xlocator = mticker.FixedLocator(np.arange(extent[0], extent[1]+1, 2))
        gl.ylocator = mticker.FixedLocator(np.arange(extent[2], extent[3]+1, 2))
        ax.set_xticks(list(range(extent[0], extent[1]+1, 2)), crs=ccrs.PlateCarree())
        ax.set_yticks(list(range(extent[2], extent[3]+1, 2)), crs=ccrs.PlateCarree())
        
        lon_formatter = LongitudeFormatter(zero_direction_label=True)  # zero_direction_label用来设置经度的0度加不加E和W
        lat_formatter = LatitudeFormatter()
        ax.xaxis.set_major_formatter(lon_formatter)
        ax.yaxis.set_major_formatter(lat_formatter)

cc1 = plot_cont(u_500, v_500, 
                [0.1, 0.1, 0.8, 0.8], [110, 120, 23, 31], 
                show=True, title='(a)500hPa', titleposition=[0.04, 1.05])

# ###############################################################################
# # 存图
plt.savefig(figpath+'213 wrf风场', dpi=600, bbox_inches = 'tight')
plt.show()