GrADS中格点插值到站点(gr2stn)的详细方法

mofangbao

先推荐一下MeteoInfo中的实现方法哦，很不错啦
http://bbs.06climate.com/forum.php?mod=viewthread&tid=1791

========================================================
这两天有朋友在讨论如何在GrADS中把格点数据插值为站点数据，之前整理GrADS命令的时候见过这个命令：gr2stn，使用这个命令即可在GrADS中将格点数据插值为站点数据，经过参考官方手册，我使用该命令进行了一些测试，现在把过程总结出来和大家交流。

先介绍下gr2stn命令的两种格式，并会在下文中分别测试两种格式的效果：
1、gr2stn(grid_expr, lon, lat, <-n>)
2、gr2stn(grid_expr, stn_expr, <-n>)

grid_expr为原始的格点数据中的变量，lon/lat为经纬度，-n是在grads2.0a6版本之后新增的一个功能，表示不进行插值，而采用最相邻的格点值来替代，stn_expr表示一个站点数据文件中的变量，如果现在不能理解这个参数的意思，没有关系，后面的例子看了就明白了。GrADS中的格点到站点的插值方法为双线性内插。

第一种格式：使用单点的经纬度来进行插值，该方式每次只能插值出一个站点数据；
第二种格式：使用一个变量表达式来提供站点背景数据。

先使用第一种方法来进行插值测试：
（该方法比较苛刻，插值前需要固定为Z或者T方向变化的一维数据）
使用的数据：一个规范的格点数据文件，至少需要可以再固定经度之后，Z或T方向存在变化。

我使用的格点数据如下，
Ctl文件：

dset E:\grads\gr2stn\sx6.grd

undef -9.99E+33

title  PROJECT

xdef   9 linear   80.00  2.500

ydef   5 linear   20.000  2.500

zdef    2 levels   850  200

tdef 48 linear   JAN1982      1mo

vars   4

U  2 99 u wind (m/s)

V  2 99 v wind (m/s)

H  1 99 H850

TSFC 1 99 TSFC DATA

Endvars

Gs文件：

'reinit'

'open e:\grads\gr2stn\sx6.ctl'

'set x 1'

'set y 1'

*'set z 1'

'set t 1 48'

'd mag(gr2stn(u,88,26),gr2stn(v,88,26))'

'printim e:\grads\gr2stn\sx6.png white'

;
         需要注意的就是上面维数的设置，因为输出的插值后的变量经纬度已经是固定的，因此必须要这样设置，当然你也可以不设置为1。u,v就是ctl里面的变量,做出的图形如下：

gr2stn_lonlat

下面是第二种方法：
使用的数据：
1、由于手头没有和我的站点文件匹配的格点数据（上面的范围太小，做的话需要生成合适的站点文件，后面会提到制作方法），因此使用站点数据先做了一些处理，用插值出来的格点数据进行反插，大致流程如下：
将由气候中心160站的2011年5月份降水量得来的降水量距平百分率数据转换为GrADS可以识别的二进制格式，并生成相应的.map文件，配出相应的grid文件和ctl文件（具体过程可以在本版找到相应帖子）。（这些是我由于数据的限制才做的，我的站点数据既是站点背景文件，又是数据文件，如果你有现成的格点数据是最好了）
2、站点的背景文件。
和站点插值到格点类似，也需要提供一个站点的背景文件，我就是用的160站作为背景站点资料，好了，到现在我已经有一个2011年5月的降水量距平百分率的站点二进制文件，一个grid.ctl文件，一个grid.grd文件，一个站点文件的map和ctl，如果你是用的现成的格点文件来插值为站点，则grid.ctl和grid.Grd就是你的格点数据源文件了，请千万不要混淆。
作图的流程很简单
1、  打开你的格点数据文件和站点数据文件的ctl
2、  我先插值出格点文件（你直接用的格点文件就跳过这一步）
3、  插值并作图
其实和站点插格点的过程是一样的。
下面是这个gs文件：

'reinit'

*文中的第一步

'open E:\grads\gr2stn\grid.ctl'

'open E:\grads\gr2stn\result.ctl'

*文中的第二步，插值出格点（如果自己提供格点文件那么这一步省略）

'define a=oacres(g,result.2)'

*文中第三步，插值出站点文件并作图

'd gr2stn(a,result.2)'

'printim  E:\grads\gr2stn\sta_result.png white'
；
现在出来的就是一个插值出来的站点数据图啦，如下所示：

gr2stn_sta

看来长江流域还是很多负距平的，再看一下它和源数据的误差（这其中的误差包括oacres和gr2stn两次总共所造成的误差，如果直接提供的是格点文件，就是只有gr2stn的误差了）
把原来的gs改成如下

'reinit'

*文中的第二步

'open E:\grads\gr2stn\grid.ctl'

'open E:\grads\gr2stn\result.ctl'

*文中的第一步，插值出格点（如果自己提供格点文件那么这一步省略）

'define a=oacres(g,result.2)'

*文中第三步，插值出站点文件并作图

*'d gr2stn(a,result.2)'

*'printim  E:\grads\gr2stn\sta_result.png white'

*再求一下误差，插回去和原来的对比下，插值这么多次，必然有误差了

'd result.2-gr2stn(a,result.2)'

'set gxout shaded'

'cnbasemap b'

'cbar'

;
那个b就是算出来的误差啦，图形如下：

gr2stn_comp_sta

不分析了…
当然也可以作成阴影图的形式啦：

gr2stn__comp_shaded

恩，差不多就这样了，总结一下：
1、  使用给定lon/lat的形式适合时间序列或者沿垂直方向变化的插值，而使用站点背景文件的形式使用范围好像更大一点；
2、  给定的站点背景文件需要事先生成好map映射；
3、  还没找到如何可以提取出插值出来的数据，希望能有人接着测试下去，分享出方法（方法已有，请论坛搜索）。
4、  我还没有测试多个时间和高度序列的批处理情形，因此如果多个时间序列的话在插值的时候是不是要加入类似 gr2stn(grid_var,sta_var.2(t=’t’))的语句。


附：制作站点背景文件
1、  这个其实和站点插格点是一样的，首先需要确定你的数据的范围，然后提取出你数据范围内的站点信息，这一步可以通过MeteoInfo来完成（http://bbs.06climate.com/forum.php?mod=viewthread&tid=1119&extra=page%3D1），也可以自己手工提取；
2、  提取后使用fortran来将站点信息写入GrADS可以识别的二进制文件，写一个变量就行，变量值无所谓，因为用的时候只是把经纬度信息当做一个模板来使用，不涉及变量值。
3、  配置ctl文件，生成map映射文件。

上面说的很清楚啦，如果还想使用我使用的数据进行测试，需要一定论坛贡献才能下载，人人为我，我为人人嘛，欢迎分享文件或者经验，这样会增加贡献啦
gr2stn.zip (35.32 KB, 下载次数: 259, 售价: 5 贡献)

2011-8-12 10:09 上传

点击文件名下载附件
GrADS中的格点到站点插值(gr2stn)以及站点到格点插值(oacres)方法.doc (90 KB, 下载次数: 434, 售价: 2 贡献)

2011-8-11 19:45 上传

点击文件名下载附件

善人/jw

要这么多贡献啊~~

尽头的尽头

我要努力贡献。。。。。

/sun满天星

哦啊 希望在我攒够贡献前 帖子要永生啊

Mouse

看来我要更加努力做贡献了

传说中的谁

想说什么来着……这个是好东西，也许很快就用得着了。感谢清风的无私奉献！
PS：好东西卖钱也是我的馊主意，哈哈

Mouse

传说中的谁 发表于 2011-8-11 23:57

                               
登录/注册后可看大图

想说什么来着……这个是好东西，也许很快就用得着了。感谢清风的无私奉献！
PS：好东西卖钱也是我的馊主意 ...

你的馊主意还真不错，这下好多人要使劲赚钱才行了

zhb

好东西，求下载

zhb

怎么还是下载不了啊

喜上眉梢

好东西，慢慢学用！