MATLAB
鉴于气象圈中使用matlab的比较多,先说一下使用matlab如何读取 grib2 格式数据。
使用matlab读取grib2格式数据需要用到 nctoolbox 工具箱[注1]
下载好后,解压,然后将路径添加到matlab搜索路径中即可。演示使用 MATLAB 版本为 2016a
开始读取数据之前,首先要执行以下语句:
setup_nctoolbox();
加载成功后即可进行数据读取
filename = 'E:\MATLAB\data\fnl_20160623_06_00.grib2';
data = ncdataset(filename);
lon = squeeze(data.data('lon'));
lat = squeeze(data.data('lat'));
t = squeeze(data.data('Temperature_surface'));
读取数据要知道相应的变量名,变量名的查看可以通过以下方式查看:
data.variables
比如下图,黑色框表示如何获取变量名信息,红色框除了包含变量名信息之外还包含了每个变量的其他信息,比如单位等,绿色框是变量名(没有全部列出)
数据读取好后就是画图了,下面展示一下绘制的全球温度分布及skewt图(仅展示绘图效果)
全球温度分布
使用FNL数据绘制某一经纬度点的 T-lnP 图
NCL
fnl再分析数据是转WRF模式经常用到的数据,而转模式就不得不提NCL这一绘图工具了。演示使用 NCL 6.3.0
NCL读取grib2格式数据同样非常方便。
fin = "fnl_20160623_06_00.grib2"
f = addfile(fin,"r")
读取数据之前,也需要确定变量名,使用 print(f) 可以查看变量名及其维度等信息,然后就可以读取数据并绘图了。
下面读取了压力,经纬度数据,其它变量同样可以按照以下方式读取
p = f->lv_ISBL0
lat = f->lat_0
lon = f->lon_0
绘图和常规绘图操作类似。
Python
python读取grib2格式数据主要有两种方式,1) 使用 pygrib 读取 2) 使用PyNio
由于 PyNio 的读取方式和 NCL 非常相似,这里主要说一下使用 pygrib 如何读取。
如果你没有 pygrib的话,需要先安装之,关于安装这里按下不表。
注:演示使用 Linux , Ipython2.7, pygrib 2.0.2 ,matplotlib 1.5.3
import pygrib
data = pygrib.open('fnl_20160623_06_00.grib2')
# 获取变量信息
variables = []
for var in data:
variables.append(var)
获取变量后,可查看含有哪些变量。注意:有些变量名相同,只是属于不同层的数据,所以要注意选择。
variables.sort
Out[7]: <function sort>
variables
Out[8]:
[1:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 100000 Pa:fcst time 0 hrs:from 200712061200,
2:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 97500 Pa:fcst time 0 hrs:from 200712061200,
3:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 95000 Pa:fcst time 0 hrs:from 200712061200,
4:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 92500 Pa:fcst time 0 hrs:from 200712061200,
5:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 90000 Pa:fcst time 0 hrs:from 200712061200,
6:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 85000 Pa:fcst time 0 hrs:from 200712061200,
7:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 80000 Pa:fcst time 0 hrs:from 200712061200,
8:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 75000 Pa:fcst time 0 hrs:from 200712061200,
9:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 70000 Pa:fcst time 0 hrs:from 200712061200,
10:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 65000 Pa:fcst time 0 hrs:from 200712061200,
11:Geopotential Height:gpm (instant):regular_ll:isobaricInhPa:level 60000 Pa:fcst time 0 hrs:from 200712061200,
many variables
由上面输出信息可以看出,其中包含了大量关于变量的信息。下面读取CAPE相关信息,通过 name 参数指定要读取的变量名,从而获取数据:
cape = data(name="Convective available potential energy")[0]
cape.data
Out[10]: <function data>
cape.data()
Out[11]:
(array([[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 0., 0., ..., 0., 0., 0.],
...,
[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 0., 0., ..., 0., 0., 0.],
[ 0., 0., 0., ..., 0., 0., 0.]]),
array([[ 90., 90., 90., ..., 90., 90., 90.],
[ 89., 89., 89., ..., 89., 89., 89.],
[ 88., 88., 88., ..., 88., 88., 88.],
...,
[-88., -88., -88., ..., -88., -88., -88.],
[-89., -89., -89., ..., -89., -89., -89.],
[-90., -90., -90., ..., -90., -90., -90.]]),
array([[ 0., 1., 2., ..., 357., 358., 359.],
[ 0., 1., 2., ..., 357., 358., 359.],
[ 0., 1., 2., ..., 357., 358., 359.],
...,
[ 0., 1., 2., ..., 357., 358., 359.],
[ 0., 1., 2., ..., 357., 358., 359.],
[ 0., 1., 2., ..., 357., 358., 359.]]))
capes = cape.data()[0]
capes.shape
Out[13]: (181, 360)
lats, lons = cape.latlons() # 获取经纬度信息
lats.shape
Out[15]: (181, 360)
绘制CAPE全球分布图(仅绘制简单的分布图)
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ax.contourf(lons, lats, capes)
pygrib 除了可以读取grib格式数据之外,还可以在 grib1 和 grib2 之间互相转换。
总结
读取数据的方式多种多样,只要能实现要求即可。
注1:https://github.com/nctoolbox/nctoolbox