提取GrADS插值结果（oacres/gr2stn/grads2ascii）

mofangbao

昨天看到@topmad的提取surfer的插值结果（点击查看） ，想起来GrADS也可以提取插值结果，用到的是@兰溪之水的grads2ascii脚本，该脚本的下载地址如下：
http://bbs.06climate.com/forum.php?mod=viewthread&tid=5271

在使用此贴的方法之前，请确定你已经会使用GrADS进行站点数据作图，已经能够理解gr2stn的用法，如果还没有，请看以下帖子的介绍：
http://bbs.06climate.com/forum.php?mod=viewthread&tid=4903
http://bbs.06climate.com/forum.php?mod=viewthread&tid=1780

有了上面的3个知识，提取插值结果就是把3个方法结合起来即可，非常简单！

一、提取站点数据插值为格点数据的结果
本例子需要使用一个站点数据和一个格点背景数据，站点数据使用的是160站的降水数据，已经处理为GrADS可以支持的二进制形式，并配置了ctl文件，生成了map文件，格点背景为70E-140E,15N-55N的2.5*2.5的数据。

基本步骤：
1、打开格点背景数据和站点数据
2、调用插值函数插值到格点
3、调用输出脚本进行输出

代码如下：

1. # every time before you run the script,please del the text files created before
   
2. 'reinit'
   
3. 'open E:\program\example\interp\grid.ctl'
   
4. 'open E:\program\example\interp\rain.ctl'
   
5. 'define gridData = oacres(g,r.2)'
   
6. #export grid data
   
7. 'grads2ascii gridData E:\program\example\interp\grid_rain.txt %g 71'
   
8. 
   

# every time before you run the script,please del the text files created before<br /> 'reinit'<br /> 'open E:\program\example\interp\grid.ctl'<br /> 'open E:\program\example\interp\rain.ctl'<br /> 'define gridData = oacres(g,r.2)'<br /> #export grid data<br /> 'grads2ascii gridData E:\program\example\interp\grid_rain.txt %g 71'<br /> <br />

关键的就是define一个格点变量，然后使用外部脚本输出到文件即可。

输出的结果预览：


处理输出的数据：
这个数据的处理可以使用fortran编程实现，方法再此不做描述，自己输出的时候控制好每行的列数，尽量和你的径向格点数保持一致（比如代码中的71），这样看起来会很直观，处理也方便。

二、提取格点数据插值为站点数据的结果
这个例子中需要使用一个格点数据和一个站点背景数据，格点数据使用的是ncep再分析资料的地面气压数据pres.mon.mean.nc，站点背景使用的是中国的160站数据（直接使用上文中的rain数据即可，也可以自己另外制作一个，方法和站点作图的步骤相同）

基本步骤：
1、打开需要插值的海温数据，设置好经纬度
2、定义一个格点变量存放当前的格点数据
3、关闭nc数据，重置GrADS设置
4、打开站点背景数据，调用插值函数，并输出到文件

代码如下：

1. # every time before you run the script,please del the text files created before
   
2. 'reinit'
   
3. 'sdfopen E:\program\example\interp\pres.mon.mean.nc'
   
4. 'set lon 70 140'
   
5. 'set lat 15 55'
   
6. 'define temp=pres'
   
7. 'close 1'
   
8. 'reset'
   
9. 'open E:\program\example\interp\rain.ctl'
   
10. #now interp the grid data to station data and export
    
11. 'grads2ascii gr2stn(temp,r) E:\program\example\interp\station_rain.txt'
    
12. ;
    
13. 
    
14. 
    

# every time before you run the script,please del the text files created before<br /> 'reinit'<br /> 'sdfopen E:\program\example\interp\pres.mon.mean.nc'<br /> 'set lon 70 140'<br /> 'set lat 15 55'<br /> 'define temp=pres'<br /> 'close 1'<br /> 'reset'<br /> 'open E:\program\example\interp\rain.ctl'<br /> #now interp the grid data to station data and export<br /> 'grads2ascii gr2stn(temp,r) E:\program\example\interp\station_rain.txt'<br /> ;<br /> <br /> <br />

这里没有同时打开nc数据和站点背景，测试过的朋友会发现，如果同时打开这两个数据，需要做的工作很多，很容易会出现第二个打开的文件数据值为空，导致输出失败，所以采取一个个的打开，如果批量做的话，可以在循环里面进行打开和关闭，虽然这样效率会低一点，但是经过我自己的测试，还是很快的。

输出结果预览：


处理输出的数据：
这里的输出数据每两行为一个站点的数据，这两行中的第一行是站号，经纬度等信息，第二行才是数据信息，知道这个再编程处理就很方便了，大家可以按照自己的需求编写代码
最后要提醒注意的就是，兰溪的那个脚本中，文件输出采用的是append的方式，因此再运行第二遍之前，请先删除原来的输出数据，否则会直接加到源数据的后面，并不会覆盖旧数据。

whmooncool

这个不错，留着~

qxtlyf

好东西，谢谢了

lyroot@sina.com

楼主  顶一个！

lyroot@sina.com

楼主  顶一个！

lz007700

脚本很给力，我现在正好在处理格点插值为站点的问题，对于只有单层单一时刻的数据，比如地形数据，用来画平面上任一条斜线上地形的变化，是不是就要先生成这条斜线上各点的经纬度信息，当成站点（只是这些站点正好在一条直线上）后，再把地形格点值插到站点上，画出剖面变化（纵轴为高度，横轴为经纬度）。这么做好麻烦啊，楼主有没有更简单点儿的方法呢？类似于画斜剖面的方法？我试了试collect函数不行。

mofangbao

lz007700 发表于 2012-5-16 21:53

                               
登录/注册后可看大图

脚本很给力，我现在正好在处理格点插值为站点的问题，对于只有单层单一时刻的数据，比如地形数据，用来画平 ...

你的数据能打个包传上来吗

zytfay

谢谢楼主！

竹下妖

lz007700 发表于 2012-5-16 21:53

                               
登录/注册后可看大图

脚本很给力，我现在正好在处理格点插值为站点的问题，对于只有单层单一时刻的数据，比如地形数据，用来画平 ...

这个问题解决了么？我也在画地形的斜线剖面呢，collect只适用于3维数据。

luoziwuhui

非常不错的东西，对我帮助很大，感谢楼主大大的无私分享！O(∩_∩)O哈哈哈~楼主威武.谢谢！lz太给力了。辛苦了