专栏首页气象杂货铺真・WRF模式后处理之Python版

真・WRF模式后处理之Python版

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。上述脚本中并未使用此库进行后处理。关于此库的介绍在后面 —>

本文分享自微信公众号 - 气象杂货铺(meteogs),作者:lightning

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2017-06-01

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 使用 Cartopy 和 netCDF4 可视化 WRF 模式数据

    对比使用 Basemap,gdal 和 Cartopy,netCDF4 读取 WRF 模式数据并绘图。

    bugsuse
  • 使用 Basemap 和 Cartopy 绘制子图实例

    平时绘制地图时,经常会将多个图放到同一个 figure 中,而这些图的地图范围通常是相同的,所以可以设置共享 x-y 轴。

    bugsuse
  • 建议收藏!Matplotlib常见组件设置整理

    继上一篇文章为大家介绍了plt和ax绘图的区别后,这篇文章结合我自己的一些使用经历,为大家整理了Matplotlib中比较常用的一些组件设置。

    bugsuse
  • 求100以内所有奇数的和,存于字变量X中。

    注释:在debug中显示的是十六进制,可以看到bx中为09c4 换算成10进制

    cuptobjut
  • testNG的高级用法 --DataProvider

    @DataProvider Method参数 数据提供者的第一个参数是java.lang.reflect.Method,TestNG传递这个将调用的测试方法。如...

    千往
  • 教程 | 5种快速易用的Python Matplotlib数据可视化方法

    选自towardsdatascience 作者:George Seif 机器之心编译 参与:刘晓坤、思源 数据可视化是数据科学家工作的重要部分。在项目的早期...

    机器之心
  • 有这5小段代码在手,轻松实现数据可视化(Python+Matplotlib)

    大数据文摘
  • 第24次文章:结构性模式(续)

    本周把结构性模式的剩余几种模式学完了,并且开始接触行为型模式。为便于对文章的分类,本周的文章就只把结构性模式全部扫尾处理。行为型模式放在下周的文章中一起介绍吧!

    鹏-程-万-里
  • 5个快速而简单的数据可视化方法和Python代码

    数据可视化是数据科学家工作的重要组成部分。在项目的早期阶段,你通常会进行探索性数据分析(EDA),以获得对数据的一些见解。创建可视化确实有助于使事情更清晰和更容...

    AI算法与图像处理
  • 这5小段代码轻松实现数据可视化(Python+Matplotlib)

    本文要讲的是Matplotlib,一个强大的Python可视化库。一共5小段代码,轻松实现散点图、折线图、直方图、柱状图、箱线图,每段代码只有10行,也是再简单...

    昱良

扫码关注云+社区

领取腾讯云代金券