WRFOUT风向变量逐时次作差绘图

用户11172986

发布于 2024-06-20 17:28:45

5040

前言

读者来信说：风电场分析需要看两个时次的风向差。同时从“wrfout中提取变量，然后用08:10的风向wdir【ncl函数wind_direction（u，v，0）】减去08:00时刻的风向，做上循环语句do，就会出现差一个数值对不上的情况。笔者对ncl不太熟悉。但是以上功能实现Python不需要循环。因为wrfout的变量是xarray格式，想必大家知道要用哪个函数了。没错，就是xarray.diff() 废话半天了，开始写代码吧。

导入库与读取变量

In [2]:

# 导入数据读取模块
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 [3]:

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([os.path.join(wrfout_path, f) for f in os.listdir(wrfout_path) if f.startswith(filename_prefix)])
wrf_list = [Dataset(f) for f in wrf_files] 
# 提取变量
uvmet10 = getvar(wrf_list, 'uvmet10_wspd_wdir', timeidx=ALL_TIMES, method='cat')[1]
uvmet10

Out[3]:

以上代码提取了wrfout数小时的uvmet10_wspd_wdir变量，实际上是离地10m的风速风向，取了第二个维度，即是只取了风向。

这个看数值就知道是角度了，哪有两百多米/秒的风.

关于风变量的讲解可以看看https://www-k12.atmos.washington.edu/~ovens/wrfwinds.html

简单绘图

In [4]:

# #计算逐小时风向变化
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[i].plot.contourf(ax=axes[row, col], 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 [5]:

wdir = getvar(wrf_list, 'wspd_wdir', timeidx=ALL_TIMES, method='cat')[1]
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[5]:

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)

    wdir100change[i].plot.contourf(ax=axes[row, col], 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()