【探讨和分享】一个利用GRADS插值到站点的方法

meehooqq

例如：
有格点描述文件rain.ctl如下：
dset ^rain.grd
undef 9.999E+20
ydef 28 linear 15.000000 1.5
xdef 48 linear 69.000000 1.5
tdef 61 linear 12Z30apr2013 6hr
zdef 1 linear 1 1
vars 1
r 0 99 ** surface Total Precipitation Rate [kg/m^2/s]
ENDVARS

另外，有站点信息文件文本station.txt ,格式如下：
57432   108.4000        30.7664
57516   106.4611        29.5758
57520   107.0667        29.8333
57536   108.7833        29.5333
57409   105.8000        30.1833
。。。。。。

利用GRADS读取站点信息文件station.txt,再把数据写出，对应GS脚本如下：
'reinit'
'open rain.ctl'
'set fwrite out.dat'
'set gxout fwrite'
cnt=35  ！cnt表示要读取的站点总个数
i=1
while (i<=cnt)
aa=read('station.txt')
aa1=sublin(aa,2)
'query 'aa1
alon=subwrd(aa1,2)
alat=subwrd(aa1,3)
'query 'alon
'query 'alat
'set lon 'alon
'set lat 'alat
'set t 1 61'
'define tmp1=r'
'd tmp1'
i=i+1
endwhile
ret=close('station.txt')
'disable fwrite'
'reinit'
'quit'


这样得到的数据就是插值出来的站点数据了。另外，发现如果原格点距离较大，而站点相对比较密集，插值出来的站点上的值没有变化，于是就修改了下CTL描述文件，加入pdef描述，再将格点加密，修改后的CTL文件如下，这样提取出来的值就不一样了。

dset ^rain.grd
undef 9.999E+20
pdef 45 28 lcc 33.5 102.75 23 13.5 60 30 102.75 165000.0 165000.0
ydef 200 linear 15.000000 0.2025
xdef 330 linear 69.000000 0.2
tdef 61 linear 12Z30apr2013 6hr
zdef 1 linear 1 1
vars 1
r 0 99 ** surface Total Precipitation Rate [kg/m^2/s]
ENDVARS

由于没有比较这种方法和用站点插值函数得到的结果哪个精度较高，所以，还请大家一起探讨。

river

首先谢谢楼主分享这个格点插值到站点方法。但是我有一些不太明白，后面你说改变ctl描述，加入pdef 45 28 lcc 33.5 102.75 23 13.5 60 30 102.75 165000.0 165000.0，而且还加密了网格。你用的是同一个资料，那ctl不应该也是一样的么，这样改变之后，对于原始资料的描述就不准确了啊

mofangbao

楼主，其实你这是取临近值的方法吧  建议你把数据先在屏幕上显示出来，然后设置一个固定的经纬度 输出值 观察一下取的是哪个值 应该没进行插值吧 我也还没具体测试

liutaoz21

同river，不明白这个ctl为啥要这样改

meehooqq

mofangbao 发表于 2013-6-1 09:40

                               
登录/注册后可看大图

楼主，其实你这是取临近值的方法吧  建议你把数据先在屏幕上显示出来，然后设置一个固定的经纬度 输出值 观 ...

我测试了下，数值变了，不过后来我又发现问题，就是感觉写成PDEF出来的图和不用PDEF出来的图不能完全重合，所以我现在也不感确定我这个方法是否可行了。

meehooqq

river 发表于 2013-6-1 08:04

                               
登录/注册后可看大图

首先谢谢楼主分享这个格点插值到站点方法。但是我有一些不太明白，后面你说改变ctl描述，加入pdef 45 28 lc ...

因为最开始的时候是提取WRF结果的插值，模式网格比较细，所以就算站点比较密提出来的数值也都不一样，后来我在处理粗网格数据的时候，用同样的方法做插值，发现由于网格粗了，插值出来很多站点的值是一样的，于是就想到了原来WRF转GRADS的时候也是用的PDEF参数，PDEF参数里面的格点和模式原有格点数一致，但是XDEF和YDEF却要密很多。于是就想试着用加入一条PDEF的参数来做一下，发现修改后输出的结果不一样。但是现在我又发现一个问题，就是修改了PDEF后，好像出来的图和不用PDEF的图有些差别，不能完全重合。所以，我在想这种方法是不是还是不适合把粗格点的数插值到细格点里面去。所以，对于这个PDEF参数，估计还得查查具体意思，还是说这个参数设置出来只适合话LAMBERT的投影

meehooqq

liutaoz21 发表于 2013-6-1 10:32

                               
登录/注册后可看大图

同river，不明白这个ctl为啥要这样改

因为修改后插值出来的各个站点的值完全都不一样了。所以才想这种办法是否可行，拿出来和大家讨论下。毕竟加入了PDEF后不能完全保证图形和不加PDEF出来的图完全一致，这也是我有点纠结的地方。

river

meehooqq 发表于 2013-6-3 23:13

                               
登录/注册后可看大图

因为最开始的时候是提取WRF结果的插值，模式网格比较细，所以就算站点比较密提出来的数值也都不一样，后来 ...

哦，好曲折啊。但是有一点我要说的，ctl必须如实描述你的资料。如果想把格点资料的分辨率增加，那就必须通过改变模式的分辨率或者利用fortran等程序把你的资料加密，然后再配上ctl。单纯的改变ctl不会改变原有资料的分辨率，只能是画出错误的图，得到错误的结论

meehooqq

river 发表于 2013-6-4 12:07

                               
登录/注册后可看大图

哦，好曲折啊。但是有一点我要说的，ctl必须如实描述你的资料。如果想把格点资料的分辨率增加，那就必须通 ...

是啊，我也考虑到这一点，所以暂时不用这种方法来提取粗格点的资料。

流浪的心

grads lambert 投影用原始格点作图