Grads处理GRIB格式文件的准备
LYanbing
2007-7-2初稿，2007-7-12修改1 开场说明
WHAT IS GRIB?
GRIB (GRIdded Binary) is an international, public, binary format for the efficient storage of meteorological/oceanographic variables. Typically, GRIB data consists of a sequence of 2-D (lon,lat) chunks of a (in most general sense) 4-D variable (e.g., u comp on the wind = f(lon,lat,level,time)). The sequence is commonly organized in files containing all variables at a particular time (i.e., 3-D (lon,lat,level) volume).
大气所的NCEP再分析资料使用这种格式。
这里针对6小时一次的1°×1°，26层数据来处理。
Grads中识别路径的方式基本为Unix的方式，即路径中用斜杠/，而不是反斜杠\，cmd中也支持这种方式，所以，使用Grads及其相关组件时，指定路径用斜杠/会很方便。
Grads中用!pwd可以看到当前目录，C：盘对应/cygdrive/c/，d：盘对应/cygdrive/d/。cmd中用pwd看到的也是如此，cygdrive是怎么来的？与cygwin程序有关，它能把Unix程序嫁接到windows下使用，它的目录系统以/cygdrive/为根目录。
如果ctl文件中，数据文件指定不是全路径，而是^，则可以在open命令中指定上述形式的全路径，例如：
ga-> open /cygdrive/d/data/l
它等效于：
ga-> open d:/data/l
2 生成描述文件
PCGrads软件的User’s Guide中介绍了GRIB及其处理方法，但不完全。实际上有两种方法：1）利用工具生成整个文件的描述文件.ctl，再利用工具生成映射文件.idx；2）利用工具解码文件中需要使用的部分记录，建立新的数据文件，然后人工建立描述文件.ctl。
为了使用方便，环境变量Path中增加Grads可执行文件所在目录，则在cmd中其他路径下亦可访问所有该目录下的工具。
2.1 方法1
使用工具生成GRIB数据文件的描述文件ctl，之后还要使用gribmap工具生成映射文件.idx。
1. 工具的获得。
查看相关信息：
v/products/wesley/grib2ctl.html
一般网上得到的是源码文件grib2ctl.pl。
ftp://v/wd51we/wgrib.scripts/grib2ctl.pl
grib2ctl.pl是perl语言代码，转换成exe格式需要使用perl2exe工具。转换的过程还需要借助perl编译器，如ActivePerl。
perl2exe工具可从以下网址查看和下载：
www.indigostar/perl2exe.htm
www.indigostar/download/p2x-8.60-Win32.zip
ActivePerl可从以下网址下载：
ftp.activestate/ActivePerl/Windows/5.8/ActivePerl-5.8.6.811-MSWin32-x86-1222 08.msi
ActivePerl下载后安装好，perl2exe下载后解压到某个目录，把grib2ctl.pl拷贝入该目录。在cmd命令行窗口执行：
grib2ctl.pl
perl2exe
于是生成。为了方便使用，把它拷贝到Grads目录下，这里是PCGrads\win32e\。
2. 工具的使用。
使用过程中要用到wgrib工具，Grads自带了但不是最新的，最新版本需要自己下载。
查看：
v/products/wesley/wgrib.html
可选择下载源码或已编译版本。
这里使用已编译版本，下载和其依赖的cygwin1.dll文件:
ftp://v/wd51we/wgrib/machines/
ftp://v/wd51we/wgrib/machines/Windows/cygwin1.dll
把这两个文件放置在可以访问到的地方，例如和放置在同一目录下，在cmd 命令行窗口，cd到数据文件所在路径（这里是d:/data/ncep/），然后执行： grib2006060100 &l
注意，这样生成的.ctl文件还不能直接使用，需要处理：
去掉文件结尾多余部分，即ENDV ARS之后的关于Perl2exe声明的部分（因为我们按照前文地址得到的Perl2exe是试用版）。
得到的ctl文件举例见附录1。
3．gribmap的使用
在cmd命令行窗口，cd到数据文件所在路径，然后执行：
gribmap –v –l &
则，该路径下生成映射文件grib2006060100.idx，b.txt中记录了映射的过程，打开它可以查看所有field匹配的情况，全是！！！！！MATCH才是完全匹配成功。（这里发现有的行不是！！！！！MATCH，而是..... NOOOO，问题已经解决，见2.3节）
如果不需要诊断，直接执行：
gribmap –l
保证数据文件、ctl文件、idx文件，三者齐全，则Grads下可以用命令打开：
ga->open d:/data/l
4．总结
实际使用过程就是：
d:\data\ncep> grib2006060100 &l
打开ctl文件，把结尾多余信息去掉
d:\data\ncep> gribmap –v –l &
根据gribmap的映射过程文件b.txt确定匹配情况，根据需要修改ctl文件（分析和修改过程见2.3节的例子），再次gribmap
2.2 方法2
利用工具解码文件中需要使用的部分记录，建立新的数据文件，然后人工建立描述文件.ctl。这种方法需要生成数据的备份，所以会占用较多的磁盘空间。
1．提取数据
工具除了像方法1中那样被grib2ctl工具调用外，还可以单独使用，用于选择性地解码GRIB文件中的记录。
要确定选择哪些记录，需要查看记录信息列表，用wgrib来生成，在数据文件路径下执行：
wgrib –v grib2006060100 &
该文本文件就包含了数据文件中所有记录的列表，包括记录号、变量名、变量描述等。
一般我们要解码的是多条记录，可利用批处理文件。
例如，这里我们要解码T、U、V三个变量，可以通过查看记录信息列表获知对应的记录号，然后编写批处理文件，见附录2。这个批处理带有两个参数， %1表示输入的ncep 资料文件名，%2表示输出的二进制数据文件名，使用举例：
ncep.bat grib2006060100 2006060100TUV.dat
2．建立描述文件
解码后的数据文件要配上ctl文件才能使用，这要手工建立。注意ctl中要有OPTIONS yrev，表示y轴数据是反向的，因为Grads默认y轴（纬度）从南到北，而解码出的数据（纬度）排列是从北到南。前面生成的数据文件对应的ctl内容应为：
DSET  ^2006060100TUV.dat
TITLE TUV data
UNDEF -9.99E33
OPTIONS yrev
XDEF 360 LINEAR 0    1
YDEF 181 LINEAR -90 1
ZDEF 26 LEVELS 1000 975 950 925 900 850 800 750 700 650 600 550 500 450 400 350 300
250 200 150 100 70 50 30 20 10
TDEF  1 LINEAR 01jun2006 06hr
V ARS    3
temp  26  99  temperature
u    26  99  u wind
v    26  99  v wind
ENDV ARS
2.3 解决idx中不匹配问题
原来生成映射过程文件b.txt中，发现有7行不是！！！！！MATCH，而是..... NOOOO，这
个错误已经解决！
下面举例说明解决过程（针对最后一个..... NOOOO）：
调查现场线索
1. gribmap –v –l &生成的b.txt中的出事现场
..... NOOOO:  285  23782648  3  1  0  27  100 500 （后面省略）
可以发现出错的地方对应field号是285。
2. 用wgrib –v grib2006060100 &得到的数据field描述文件中，285号field附近的
现场
284:23592070:D=2006060100:GPA:1000 mb:kpds=27,100,1000:anl:"Geopotential height anomaly [gpm]
285:23698038:D=2006060100:GPA:500 mb:kpds=27,100,500:anl:"Geopotential height anomaly [gpm]
可以发现285号的field表示的是500hpa的“Geopotential height anomaly”，另有284号field表示的是1000hpa的“Geopotential height anomaly”，两者简称都是“GPA”。
3. grib2ctl grib2006060100 &l生成的l中与
“Geopotential height anomaly”及“GPA”相关的现场
GPAprs 2 27,100,0 ** (profile) Geopotential height anomaly [gpm]
可以发现，284和285号的field被grib2ctl合成了一个2层的Var，即GPAprs。
分析
现场3中，如果按照GPAprs这样的2层Var定义方式，那么对应的两个气压层次应该是zdef 26 levels中最初的2层，即1000和975，1000对应现场2中284号field是对的，但975无
法对应285号field，因为这个field表示的层次是500，这就是问题所在了！
纠正错误
修改ctl文件，把现场3中的2层GPAprs变量改成2个0层的变量，kpds要对应现场2：GPA1000hpa 0 27,100,1000 ** 1000hpa Geopotential height anomaly [gpm]
GPA500hpa 0 27,100,500 ** 500hpa Geopotential height anomaly [gpm]
然后，还要记得把变量数Vars 95改成Vars 96，保存，这时新的ctl就更正确了一些，再用gribmap诊断发现原错误行变成了！！！！！MATCH。
其他的错误行可用类似的办法修正！
需要说明的是，如果你需要使用的变量都！！！！！MATCH成功，上文的这种修改就不是
必须的，如果追求完美，那自然要修改至清一！！！！！MATCH为止，这是你的自由！
附录3是这种完美的一种实现！
附录1
dset ^grib2006060100
index ^grib2006060100.idx
undef 9.999E+20
title grib2006060100
* produced by grib2ctl v0.9.12.5p39c
dtype grib 3
options yrev
ydef 181 linear -90.000000 1
xdef 360 linear 0.000000 1.000000
tdef 1 linear 00Z01jun2006 1mo
* z has 26 levels, for prs
zdef 26 levels
1000 975 950 925 900 850 800 750 700 650 600 550 500 450 400 350 300 250 200 150 100 70 50 30 20 10
vars 95
no4LFTXsfc 0 132,1,0 ** surface Best (4-layer) lifted index [K]
no5WA V A500mb 0 230,100,500 ** 500 mb 5-wave geopot. height anomaly [gpm]
no5WA VH500mb 0 222,100,500 ** 500 mb 5-wave geopotential height [gpm]unknown怎么处理
ABSVprs 26 41,100,0 ** (profile) Absolute vorticity [/s]
CAPEsfc 0 157,1,0 ** surface Convective Avail. Pot. Energy [J/kg]
CAPE180_0mb 0 157,116,46080 ** 180-0 mb above gnd Convective Avail. Pot. Energy [J/kg] CINsfc 0 156,1,0 ** surface Convective inhibition [J/kg]
CIN180_0mb 0 156,116,46080 ** 180-0 mb above gnd Convective inhibition [J/kg] CLWMRprs 21 15
3,100,0 ** (profile) Cloud water [kg/kg]
CWATclm 0 76,200,0 ** atmos column Cloud water [kg/m^2]
GPAprs 2 27,100,0 ** (profile) Geopotential height anomaly [gpm]
HGTsfc 0 7,1,0 ** surface Geopotential height [gpm]
HGTprs 26 7,100,0 ** (profile) Geopotential height [gpm]
HGTpv2 0 7,117,2000 ** pot vorticity = 2000 units level Geopotential height [gpm] HGTpvneg2 0 7,117,34768 ** pot vorticity = -2000 units level Geopotential height [gpm] HGThtfl 0 7,204,0 ** highest trop freezing level Geopotential height [gpm]
HGT0deg 0 7,4,0 ** 0C isotherm level Geopotential height [gpm]
HGTmwl 0 7,6,0 ** max wind level Geopotential height [gpm]
HGTtrp 0 7,7,0 ** tropopause Geopotential height [gpm]
HPBLsfc 0 221,1,0 ** surface Planetary boundary layer height [m]
ICECsfc  0 91,1,0  ** surface Ice concentration (ice=1;no ice=0) [fraction]
LANDsfc 0 81,1,0 ** surface Land cover (land=1;sea=0) [fraction]
LFTXsfc  0 131,1,0  ** surface Surface lifted index [K]
O3MRprs 6 154,100,0 ** (profile) Ozone mixing ratio [kg/kg]
POTsig995  0 13,107,9950 ** sigma=.995  Potential temp. [K]

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。