请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 160047|回复: 299

[分享资料] GrADS中格点插值到站点(gr2stn)的详细方法

  [复制链接]

新浪微博达人勋

0
早起挑战累计收入
发表于 2011-8-11 19:36:52 | 显示全部楼层 |阅读模式

登录后查看更多精彩内容~

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

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

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

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

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

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

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

我使用的格点数据如下,
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'
;
         需要注意的就是上面维数的设置,因为输出的插值后的变量经纬度已经是固定的,因此必须要这样设置,当然你也可以不设置为1u,v就是ctl里面的变量,做出的图形如下:


gr2stn_lonlat

gr2stn_lonlat



下面是第二种方法:
使用的数据:
1、由于手头没有和我的站点文件匹配的格点数据(上面的范围太小,做的话需要生成合适的站点文件,后面会提到制作方法),因此使用站点数据先做了一些处理,用插值出来的格点数据进行反插,大致流程如下:
将由气候中心160站的20115月份降水量得来的降水量距平百分率数据转换为GrADS可以识别的二进制格式,并生成相应的.map文件,配出相应的grid文件和ctl文件(具体过程可以在本版找到相应帖子)。(这些是我由于数据的限制才做的,我的站点数据既是站点背景文件,又是数据文件,如果你有现成的格点数据是最好了)
2、站点的背景文件。
和站点插值到格点类似,也需要提供一个站点的背景文件,我就是用的160站作为背景站点资料,好了,到现在我已经有一个20115月的降水量距平百分率的站点二进制文件,一个grid.ctl文件,一个grid.grd文件,一个站点文件的mapctl,如果你是用的现成的格点文件来插值为站点,则grid.ctlgrid.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

gr2stn_sta



看来长江流域还是很多负距平的,再看一下它和源数据的误差(这其中的误差包括oacresgr2stn两次总共所造成的误差,如果直接提供的是格点文件,就是只有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_sta


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


gr2stn__comp_shaded

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 贡献)

评分

参与人数 13威望 +2 金钱 +37 贡献 +6 收起 理由
liulinws + 2 很给力!
zjtong + 2 + 2 很给力!
咏咏 + 2 赞一个!
979877673 + 1 很给力!
yrovl + 2 很给力!
爱智 + 2
饭团TT + 2 很给力!
jane1124 + 4 很给力!
qxtlyf + 2 很给力!
meehooqq + 2
灰太狼 + 2 很给力!
言深深 + 4 + 2 赞一个!
传说中的谁 + 2 + 10 + 2 这个不错的说

查看全部评分

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2011-8-11 19:52:21 | 显示全部楼层
要这么多贡献啊~~

点评

你又不缺贡献,主要是最近发现了一些让我非常郁闷的事情,不得已而为之,好吧,减掉一点了  发表于 2011-8-11 19:54
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2011-8-11 19:55:15 | 显示全部楼层
我要努力贡献。。。。。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2011-8-11 20:11:26 | 显示全部楼层
哦啊 希望在我攒够贡献前 帖子要永生啊

点评

贡献要分享材料才有得哦,或者,偷偷的告诉你,每天签到前两名都有5个贡献赠送,哈哈  发表于 2011-8-11 20:59
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2011-8-11 23:32:52 | 显示全部楼层
看来我要更加努力做贡献了

点评

加油加油哈 呵呵  发表于 2011-8-12 07:56
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2011-8-11 23:57:29 | 显示全部楼层
想说什么来着……这个是好东西,也许很快就用得着了。感谢清风的无私奉献!
PS:好东西卖钱也是我的馊主意,哈哈

点评

恩 给愿意互相帮助的人  发表于 2011-8-12 07:56
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2011-8-12 09:45:02 | 显示全部楼层

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

点评

大部分还是门槛很小的,原创的东西希望能有人一起分享啦  发表于 2011-8-12 10:07
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2011-8-12 10:42:43 | 显示全部楼层
好东西,求下载

点评

分享资料吧,可以得到贡献啦,同时我们还会给你加贡献的哦  发表于 2011-8-12 11:05
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2011-8-12 11:32:04 | 显示全部楼层
怎么还是下载不了啊
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2011-8-12 11:53:58 | 显示全部楼层
好东西,慢慢学用!
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表