WRF模式是数值天气预报和大气模拟系统,其开发目的就是用语研究和实际应用。运行WRF模式时,可以利用多种初始场数据来驱动,然后配置好选项之后便可以模拟天气过程(说的好像很简单的样子==)。
个例模拟结束之后怎么办呢,我们怎么知道模拟的效果究竟如何呢?既然模拟了,那么就会有数据输出,想要检验模拟效果如何,就要对模拟数据进行处理。这就是所谓的后处理过程,也就是对模拟结果的处理,提取想要的数据,并进行分析及可视化的过程。
在地学系统中,尤其是大气科学领域,对WRF模式后处理主要使用的是GrADS和NCL,而GrADS同FORTRAN一样,属于历史悠久系列产品之一。而NCL是由NCAR开发的脚本语言,近些年来用户颇多,其主要应用与地学领域,其中包含了大量的气象相关脚本,很多气象要素的计算函数,以及气象常用图形的绘制函数,官网也给出了大量的示例,在气象领域中应用广泛。既然涉及气象领域,那么自然就少不了对WRF模式的支持,NCL中有一个脚本库,专门应用于WRF模式的后处理[注1]。
Python是一门新兴编程语言,号称胶水语言,其和众多其它编程语言或工具之间都有接口,可以非常方便的在不同语言之间进行“交流”。其近几年在气象领域的发展正如火如荼,在国内也正变得炙手可热。
关于python不作过多介绍,这也不是重点。重点是使用python进行WRF后处理。
Python进行WRF模式后处理,主要使用三个库:matplotlib(python中最火的可视化库),netCDF4(处理nc文件),Basemap(处理地图投影)。
当然关于处理 nc 文件的库还有不少,这里主要以 netCDF4为主,cartopy库也可以处理地图投影(之前也有介绍过)。这里不作过多介绍,有兴趣的可以搜索一下。
下面以一个脚本讲一下后处理过程:
#!/usr/bin/env python
# -*- coding: utf-8 -*-
from matplotlib import cm,colors
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
import netCDF4 as nc
from wrfpost import sharexy
fip = "/run/media/storm/Flash/WRF/623/"
fin1 = "wrfout_d02_2016-06-23_06_00_00"
shpfn = '/home/storm/MATLAB/shp/bou2_4p'
ti = [0,6,12,18]
xs = 0; xe = -1
ys = 0; ye = -1
zs = 0; ze = -1
data = nc.Dataset(fip + fin1, "r")
truelat = data.TRUELAT1
truelat = data.TRUELAT2
stalon = data.CEN_LON
stalat = data.CEN_LAT
# 国内雷达图配色
collev = ["#FFFFFF","#98F5FF","#87CEEB","#00FF00",\
"#008B00","#FFFF00","#FFA500","#FF8C00",\
"#FF0000","#AF0000","#8B0000","#FF00FF",\
"#8A2BE2","#360CF9"]
# 文献中雷达反射率配色
# ["#9808F0", "#51A5CE", "#3455A9", "#7FC500", "#54BE09", "#39932C",\
# "#F7E731","#E6BC25", "#F78C29", "#E12C2D", "#A22429", "#4A1910",\
# "#B85BC1", "#935CB6"]
cmaps = colors.ListedColormap(collev,'indexed')
cm.register_cmap(name = 'dbzcmap',data = collev,lut = 128)
fn = "Arial"
fs = 10
ld = 0.
nrows = 2
ncols = 2
fig,axes = plt.subplots(nrows = nrows, ncols = ncols, subplot_kw= dict(aspect = 'auto'))
for i, ax in zip(np.arange(0, nrows*ncols), axes.flat):
dbz = data.variables["COMPOSITE_REFL_10CM"][ti[i],xs:xe,ys:ye]
lat = data.variables["XLAT"][ti[i], xs:xe, ys:ye]
lon = data.variables["XLONG"][ti[i], xs:xe, ys:ye]
map = Basemap(ax= ax, projection="lcc", llcrnrlon = lon[-1,0], llcrnrlat = lat[0,0],\
urcrnrlon = lon[0,-1], urcrnrlat=lat[-1,-1], lat_0 = stalat,\
lon_0 = stalon, resolution ="h")
x, y = map(lon, lat)
# 添加经度,纬度坐标,海岸线,国界线
labelsx, labelsy = sharexy(ax, nrows, ncols, i)
map.readshapefile(shpfn, 'bou2_4p', linewidth = 1, color = 'k', ax = ax)
map.drawmeridians(np.arange(int(lon[-1,0]), int(lon[0,-1])+1, 2), labels= labelsx, \
fontname = fn, fontsize= fs, linewidth = ld, ax = ax)
map.drawparallels(np.arange(int(lat[0,0]), int(lat[-1,-1])+1), labels= labelsy, \
fontname = fn, fontsize= fs, linewidth = ld, ax = ax)
con = map.contourf(x, y, dbz, np.arange(0,71,5), cmap= cmaps)
ax.set_adjustable('box-forced') # 防止绘图和坐标区域不一致
fig.subplots_adjust(top = 0.9, bottom = 0.1, left = 0.12, right = 0.77,
hspace = 0.05, wspace = 0.05)
fig.colorbar(con, ax=axes.ravel().tolist(), pad = 0.01)
#plt.savefig("panel.eps")
plt.show()
首先,需要加载使用到的库:
from matplotlib import cm,colors
import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
import numpy as np
import netCDF4 as nc
# wrfpost 为个人编写的库,其中 sharexy用于设置共享 x-y轴,之前的例子有提及
from wrfpost import sharexy
然后,读取文件并获取变量:
fip = r"/run/media/storm/Flash/WRF/"
fin1 = "wrfout_d02_2015-09-28_06_00_00"
shpfn = r'/home/storm/MATLAB/shp/bou2_4p'
# 设置每一个维度的索引
ti = [0,6,12,18]
xs = 0; xe = -1
ys = 0; ye = -1
zs = 0; ze = -1
# Dataset 用于读取文件
data = nc.Dataset(fip + fin1, "r")
# 获取投影相关属性
truelat = data.TRUELAT1
truelat = data.TRUELAT2
stalon = data.CEN_LON
stalat = data.CEN_LAT
读取相关变量:
# “ ” 内为变量名,变量名之后为要读取的数据部分,通过索引设置
dbz = data.variables["COMPOSITE_REFL_10CM"][ti[i],xs:xe,ys:ye]
lat = data.variables["XLAT"][ti[i], xs:xe, ys:ye]
lon = data.variables["XLONG"][ti[i], xs:xe, ys:ye]
设置 colormap 和其他属性
# 国内雷达图配色
collev = ["#FFFFFF","#98F5FF","#87CEEB","#00FF00",\
"#008B00","#FFFF00","#FFA500","#FF8C00",\
"#FF0000","#AF0000","#8B0000","#FF00FF",\
"#8A2BE2","#360CF9"]
cmaps = colors.ListedColormap(collev,'indexed')
cm.register_cmap(name = 'dbzcmap',data = collev,lut = 128)
fn = "Arial"
fs = 10
ld = 0.
nrows = 2
ncols = 2
绘图:
先创建figure
fig,axes = plt.subplots(nrows = nrows, ncols = ncols, subplot_kw= dict(aspect = 'auto'))
添加地图投影:
使用 Basemap 创建地图对象,并指定投影方式和地图范围
map = Basemap(ax= ax, projection="lcc", llcrnrlon = lon[-1,0], llcrnrlat = lat[0,0],\
urcrnrlon = lon[0,-1], urcrnrlat=lat[-1,-1], lat_0 = stalat,\
lon_0 = stalon, resolution ="h")
x, y = map(lon, lat) # 坐标转换
# 添加经度,纬度坐标,海岸线,国界线
labelsx, labelsy = sharexy(ax, nrows, ncols, i) # 设置共享 x-y 轴
# 添加边界线
map.readshapefile(shpfn, 'bou2_4p', linewidth = 1, color = 'k', ax = ax)
# 添加经度线和纬度线
map.drawmeridians(np.arange(int(lon[-1,0]), int(lon[0,-1])+1, 2), labels= labelsx, \
fontname = fn, fontsize= fs, linewidth = ld, ax = ax)
map.drawparallels(np.arange(int(lat[0,0]), int(lat[-1,-1])+1), labels= labelsy, \
fontname = fn, fontsize= fs, linewidth = ld, ax = ax)
# 绘图
con = map.contourf(x, y, dbz, np.arange(0,71,5), cmap= cmaps)
ax.set_adjustable('box-forced') # 防止绘图和坐标区域不一致
调整图形属性:
fig.subplots_adjust(top = 0.9, bottom = 0.1, left = 0.12, right = 0.77,
hspace = 0.05, wspace = 0.05)
fig.colorbar(con, ax=axes.ravel().tolist(), pad = 0.01)
#plt.savefig("panel.eps")
plt.show()
除了上述三个库之外,python中有一个库和NCL 中的 WRF相关库具有几乎相同的功能,就是 wrf-python。上述脚本中并未使用此库进行后处理。关于此库的介绍在后面 —>