WRFOUT风向变量逐时次作差绘图
前言
读者来信说：风电场分析需要看两个时次的风向差。同时从“wrfout中提取变量，然后用08:10的风向wdir【ncl函数wind_direction（u，v，0）】减去08:00时刻的风向，
做上循环语句do，就会出现差一个数值对不上的情况。
笔者对ncl不太熟悉。但是以上功能实现Python不需要循环。因为wrfout的变量是xarray格式，想必大家知道要用哪个函数了。
没错，就是xarray.diff()
废话半天了，开始写代码吧。
导入库与读取变量
In :
# 导入数据读取模块
import numpy as np
import pandas as pd
from netCDF4 import Dataset
import xarray as xr
# 导入可视化模块
import matplotlib.pyplot as plt
from matplotlib.cm import get_cmap
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import cartopy.io.shapereader as shpreader
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
from cartopy.mpl.gridliner import LATITUDE_FORMATTER, LONGITUDE_FORMATTER
import shapely.geometry as sgeom
import cmaps
Warning: ecCodes 2.21.0 or higher is recommended. You are running version 2.14.1
Matplotlib is building the font cache; this may take a moment.
In :import os
import numpy as np
import matplotlib.pyplot as plt
from netCDF4 import Dataset
from wrf import getvar, ALL_TIMES,interplevel
# 定义 WRF 文件夹路径和文件名前缀
wrfout_path = "/home/mw/input/typhoon9537/"
filename_prefix = "wrfout_d01_"
# 获取 WRF 文件列表，并按照文件名排序
wrf_files = sorted()
wrf_list =
# 提取变量
uvmet10 = getvar(wrf_list, 'uvmet10_wspd_wdir', timeidx=ALL_TIMES, method='cat')
uvmet10
Out:
以上代码提取了wrfout数小时的uvmet10_wspd_wdir变量，实际上是离地10m的风速风向，取了第二个维度，即是只取了风向。
这个看数值就知道是角度了，哪有两百多米/秒的风.
关于风变量的讲解可以看看https://www-k12.atmos.washington.edu/~ovens/wrfwinds.html
简单绘图
In :
# #计算逐小时风向变化
hourly_dirchange = uvmet10.diff(dim='Time')
# 绘制逐小时风向变化
fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(20, 18))
for i in range(12):
row, col = divmod(i, 3)
hourly_dirchange.plot.contourf(ax=axes, cmap=cmaps.radar, add_colorbar=True,levels=np.arange(0,360,10))
fig.suptitle('wind dir change from 2019-08-08 19:00:00 to 2019-08-09 06:00:00')
plt.show()
/opt/conda/lib/python3.7/site-packages/cmaps/cmaps.py:3517: UserWarning: Trying to register the cmap 'radar' which already exists.
matplotlib.cm.register_cmap(name=cname, cmap=cmap)
当然，风电场一般要看离地一百米左右的风，提取方法在风能密度讲过
In :
wdir = getvar(wrf_list, 'wspd_wdir', timeidx=ALL_TIMES, method='cat')
z = getvar(wrf_list, 'z', timeidx=ALL_TIMES, method='cat')
hgt = getvar(wrf_list,'HGT', timeidx=ALL_TIMES, method='cat')
gmp = z - hgt
wdir100= interplevel(wdir, gmp ,100)
wdir100
Out:
wdir100change = wdir100.diff(dim='Time')
fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(20, 18))
for i in range(12):
row, col = divmod(i, 3)
<span style="font-style: italic;"><span style="font-style: normal;"> wdir100change.plot.contourf(ax=axes, cmap=cmaps.radar, add_colorbar=True,levels=np.arange(0,360,10))
fig.suptitle('wind dir 100m change from 2019-08-08 19:00:00 to 2019-08-09 06:00:00')
plt.show()</span></span>
好像拿台风个例来画图有点极端，不过计算方法get了问题不大吧
在此可获取文件在线运行
文章来源于微信公众号：气ython风雨
我觉得还可以从这几个方面来进一步改进。
页:
[1]