如何从grib文件中提取温度、相对湿度廓线

整理BY junewu@ 2009
一、GRIB数据的获取
相关气象资料的获取可以参考如下网站：
/article/2009-07/133.html
现在用到的主要是来自ucar的NCEP再分析资料。

链接如下：
/datasets/ds083.2/
数据获取步骤如下：
1.打开上述网页，如图
2.在该页面的下部，有一个“Data Access”，单击：
3.进入如下页面，该页面上提供了两种访问方式，任选其一即可（这里选用第二种方式，即Internet Download：
4.后续页面还是询问获取方式，任选其一（这里选Static List）：
5.后续页面上列出了可获取的两种格式（都是6小时间隔的再分析资料，分Grib1和Grib2两种格式，这里我们使用Grib1格式）：
后续页面上选择2009年的数据：
后续页面上列出了整个2009年的再分析数据，找到所需日期和时间（每天4组，6小时间隔），下载即可：
二、GRIB数据的转换
必备工具：包括Ultra Edit、WinRAR、Grads(Win32版本)等；
下载下来的数据文件是tar格式的压缩文件，用winrar将其解压到一个文件夹，可以看到里面的grib文件：
将解压获得的grib文件（这里为fnl_20091004_00_00_c）拷贝到Grads安装目录的win32e 文件夹下，打开grads，输入如下命令行：
这时文件夹下就生成了20091004_00_00_c.ctl，可以用Ultra Edit打开，查看其内容。

该文件的内容主要是对数据字段的定义等。

该文件还不能直接使用，要先转换成Unix格式才能继续。

利用Ultra Edit，文件->转换->DOS到UNIX，转换ctl的编码方式，然后保存。

UltraEdit默认的“文件”菜单中可能没有这个转换项，需要自己配置一下，如图：
打开“配置”对话框以后，找到“工具栏/菜单”->“定制”->“定制菜单”：
在下面的对话框中，从右边的选项栏中找到相应的功能（这里要使用“转换”），将其拖到左边就行了（这里已经添加到左边了）：
在Grads中继续输入如下命令行：
然后在文件夹下就生成了相应的idx文件，就可以使用idx文件来打开grib文件的转换结果了。

注意：idx的文件名应该与ctl中的一致
默认情况下，经纬度为全球范围，1000mb层，各种数据字段的名称可以在相应的ctl文件中看到，比如本例中的温度廓线字段：
如果不做特别的设置，直接绘制该字段的话，得到的是1000mb层的全球温度分布图：
三、具体廓线参数的获得
廓线的话，是在1000mb到0mb的所有层，在某个经纬度位置上定义的。

以合肥地区为例（31N51，117E17），在Grads中这样设置：
得到如下的温度廓线图（31N ，117E ——对于一个格网内的位置，也许要插值来获得较精确的结果）：
概廓线数据共有26层，其分层方法参见相应的ctl 文件。

# 廓线数据的导出方法
按照手册，应该有两种方法：fwrite 或write 写出数据，前者输出的是二进制文件，需要经过转换，后者输出ASCII 文件（但是使用write 的时候总有问题，也许是我方法不对，如果你发现了为什么，请告诉我）。

输出方法如下：
这样获得的文件是二进制的，非直接可读。

但是观察输出，发现文件大小为108KB ，输出项为27，格式为低位序（Little Endian ），也就是说每一项为4个字节，可以将其假定为float ，读出观察，如下图（2009-10-04-00：00的温度廓线，26层，最后一行为一个极大值9.99e+20，作为标记）。

set gxout fwrite; set fwrite output.txt; display tmpprs; disable fwrite;
这样，从NCEP得到的数据就有了26层的坐标、值配对，如何从这26层的值内插获得
76年大气分层下对应的每层的值呢？？
我的方法是利用MA TLAB中的interp1()，将NCEP数据插值到76年大气的分层位置，
插值结果如图，蓝色线为NCEP数据提供的温度廓线，绿色线为76年标准大气中的温度廓
线，红色线为将NCEP数据插值到76年大气分层位置的廓线，可以看到实际大气廓线与模
式大气廓线差别还是很大的，通过插值，我们可以把实际大气廓线加入模式大气进行计算，
（图中76年大气的温度廓线在接近0mb时，对应温度从270K迅速减到180K，从而提高精度。

因为NCEP的数据只提供到10mb高度，所以这段仍保持原值，应该影响不大）。