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

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 11236|回复: 11

ncl和grads联合对付非SDF规则的nc文件

[复制链接]

新浪微博达人勋

发表于 2013-10-21 21:12:55 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 游子 于 2013-10-21 21:22 编辑

ncl和grads联合对付非SDF规则的nc文件       


今天在处理数据时,出现这样一个似曾相识的问题:等值线巨大无比。
        使用CESM跑出来的数据,Grads竟然读不了,原因是输出数据不是SDF规格的数据。
        SDF在家园上已经有很多介绍了,不再赘述。所以不符合SDF规定的nc数据,使用sdfopen是打不开的。
        不过没关系了,可以使用IDL、MATLAB或者NCL来读取这类数据。把这类数据先转成二进制数据,然后再处理。以NCL为例子。
        比如,output.nc是实验跑出来的数据,由于不符合SDF规格,所以你用grads是打不开的,及时加上描述文件(Control文件)使用xdfopen仍然打不开。所以我们先用ncl转换成二进制文件。代码么,就贴在下面吧。
………………………………………………………………
do t=0,dimsizes(time)
fbindirwrite(hgt(i,:,:,:))
fbindirwrite(ta(i,:,:,:))
fbindirwrite(ts(i,:,:))
end do
………………………………………………………………

配上一个CTL文件output.ctl

比如我的一组实验数据是

dset output.grd
undef 9.96921e+36
title data_waccm
*options direct
xdef 144 linear 0.000 2.5   
ydef  96 levels -90  -88.1052631578947  -86.2105263157895  -84.3157894736842  
    -82.4210526315789  -80.5263157894737  -78.6315789473684  
    -76.7368421052632  -74.8421052631579  -72.9473684210526  
    -71.0526315789474  -69.1578947368421  -67.2631578947368  
    -65.3684210526316  -63.4736842105263  -61.5789473684211  
    -59.6842105263158  -57.7894736842105  -55.8947368421053  -54  
    -52.1052631578947  -50.2105263157895  -48.3157894736842  
    -46.421052631579  -44.5263157894737  -42.6315789473684  
    -40.7368421052632  -38.8421052631579  -36.9473684210526  
    -35.0526315789474  -33.1578947368421  -31.2631578947368  
    -29.3684210526316  -27.4736842105263  -25.5789473684211  
    -23.6842105263158  -21.7894736842105  -19.8947368421053  -18  
    -16.1052631578947  -14.2105263157895  -12.3157894736842  
    -10.4210526315789  -8.52631578947368  -6.63157894736842  
    -4.73684210526316  -2.8421052631579  -0.94736842105263  0.94736842105263  
    2.84210526315789  4.73684210526315  6.63157894736841  8.52631578947368  
    10.4210526315789  12.3157894736842  14.2105263157895  16.1052631578947  
    18  19.8947368421053  21.7894736842105  23.6842105263158  
    25.578947368421  27.4736842105263  29.3684210526316  31.2631578947368  
    33.1578947368421  35.0526315789474  36.9473684210526  38.8421052631579  
    40.7368421052632  42.6315789473684  44.5263157894737  46.4210526315789  
    48.3157894736842  50.2105263157895  52.1052631578947  54  
    55.8947368421053  57.7894736842105  59.6842105263158  61.578947368421  
    63.4736842105263  65.3684210526316  67.2631578947368  69.1578947368421  
    71.0526315789474  72.9473684210526  74.8421052631579  76.7368421052632  
    78.6315789473684  80.5263157894737  82.4210526315789  84.3157894736842  
    86.2105263157895  88.1052631578947  90 ;  
zdef 49 levels 1000  925  850  700  600  500  400  300  250  200  150  100  70  50  30  20
          10  7  5  3  2  1  0.7  0.5  0.3  0.2  0.1  0.07  0.05  0.03  0.02
          0.01  0.007  0.005  0.003  0.002  0.001  0.0007  0.0005  0.0003  0.0002
          0.0001  0.00007  0.00005  0.00003  0.00002  0.00001  0.000007  0.000005
tdef 211 linear jan0001 1mo
vars 8
z 49 -999 geopotential height
t 49 -999 air temperature
ts 0 -999 surface temperature
endvars
;

注意了,模式的输出没有缺测值,在计算的过程中,关闭了外插,出现了缺测。
默认值是9.96921e+36
这就是问题所在。在下面还要说到。注意了!!!

如果此时你希望得到一个气候场,12个月比如,那么此时的脚本如下:climate.gs
'reinit'  
'open output.grd'
'set t 1'
levels='1000  925  850  700  600  500  400  300  250  200  150  100  70  50  30  20 10  7  5  3  2  1  0.7  0.5  0.3  0.2  0.1  0.07  0.05  0.03  0.02 0.01  0.007  0.005  0.003  0.002  0.001  0.0007  0.0005  0.0003  0.0002 0.0001  0.00007  0.00005  0.00003  0.00002  0.00001  0.000007  0.000005'
'! rm -rf output_climatology.grd'
'set fwrite output_climatology.grd'
'set gxout fwrite'
'set lev 1000 0.000005'
'set y 1 96'
'set x 1 144'
******* ====>> 12 months for climatology
tt=1
while(tt<=12)
say '================================'
say 'month='tt
******* ====>> 49 labels
**
i=1
while(i<=49)
say 'levels='subwrd(levels, i)
'set lev 'subwrd(levels, i)
'set y 1 96'
'set x 1 144'
'set t 1'
'zclm=ave(z,t='tt',t=211,12)'
'd zclm'
i=i+1
endwhile
**
i=1
while(i<=49)
say 'levels='subwrd(levels, i)
'set lev 'subwrd(levels, i)
'set y 1 96'
'set x 1 144'
'set t 1'
'tclm=ave(t,t='tt',t=211,12)'
'd tclm'
i=i+1
endwhile
**
i=1
while(i<=49)
say 'levels='subwrd(levels, i)
'set lev 'subwrd(levels, i)
'set y 1 96'
'set x 1 144'
'set t 1'
'uclm=ave(u,t='tt',t=211,12)'
'd uclm'
i=i+1
endwhile
**
say 'surface temperature'
'set y 1 96'
'set x 1 144'
'set t 1'
'tsclm=ave(ts,t='tt',t=211,12)'
'd tsclm'
tt=tt+1
endwhile
;

生成了一组新的气候数据,名字就是output_climatology.grd

我们在配上一个ctl文件即可output_climatology.ctl
其中的211是比较随意的时间长度,可以根据需要修改。
此时output_climatology.ctl应该与output.ctl基本一样,除了dset和tdef
dset output_climatology.grd
……
……
tdef 12 linear jan0001 1mo
……
……

如果此时你就绘图,会惊奇的发现线条是不对的,
WHY ????????

原来如此^^^^^^^^^^^^^^^^^^^^^^^
在grads中,如果我们输出数据时不设定undef的数值,那么默认值就是
-999000000,所以在output_cliamtology.ctl中应该是
undef -999000000

或者我们可以保留原来的缺测,只要在climate.gs脚本中加一行,一定要加在
display的前面:
set undef file 1
或者
set undef dfile
这是针对只打开一个文件的情况,如果你打开了两个文件,想把输出缺测设定为第二个文件的缺测值
set undef file 2
或者直接索性设定具体指
set undef value
这样就大功告成了。
grads的好处在于出图快,适合快速检验你的实验是否正确。
ncl的属性多,画图要花掉大量的时间来修图。
所以即使学会了ncl,grads也不必急于丢掉。
当然了,在这里班门弄斧了。




评分

参与人数 2金钱 +33 贡献 +13 收起 理由
mofangbao + 15 + 5 探索精神~
kongfeng0824 + 18 + 8

查看全部评分

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

新浪微博达人勋

发表于 2013-10-21 21:28:08 | 显示全部楼层
只是为了画图快看个结果的话
不如试试meteoinfo和Panoply 后者可以打开几乎一切NC文件并出图。。前者也能打开部分grads打不开的nc 画图自然是很轻松自在。
也省去了楼主NCL转换的复杂过程。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2013-10-21 21:39:00 | 显示全部楼层

再简单一点直接用更懒惰的软件 : ncBrowse.
转换的过程也是思考的过程喔,为了出图而出图就没劲了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-10-21 21:44:57 | 显示全部楼层
顶!有时候处理这些资料挺烦的
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-10-21 22:09:29 | 显示全部楼层
游子 发表于 2013-10-21 21:39
再简单一点直接用更懒惰的软件 : ncBrowse.
转换的过程也是思考的过程喔,为了出图而出图就没劲了

那个和panoply差不多懒惰啦=。= 楼主的方法自然是很赞的,我也用其类似的方法用NCL转换过。。。
楼主的目的也就是为了简单的出图看看实验结果,免去NCL繁琐的设置不是吗~那么用懒惰的方法出图看看也挺好的~
密码修改失败请联系微信:mofangbao

新浪微博达人勋

0
早起挑战累计收入
发表于 2013-10-22 08:55:09 | 显示全部楼层
谢谢楼主前来分享~发挥各种工具的优势~
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2013-10-22 08:57:12 | 显示全部楼层
mofangbao 发表于 2013-10-22 08:55
谢谢楼主前来分享~发挥各种工具的优势~

清风大哥过奖了。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-10-28 14:50:08 | 显示全部楼层
经验值得借鉴
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2013-12-11 22:17:53 | 显示全部楼层
我是菜鸟,正在奇怪为什么会打不开自己的.nc格式的文件呢,原来是这个原因啊,我的也是cesm出的运行结果~学习了~~技多不压身,能多了解一些软件也是好的~~谢谢楼主的分享,和楼上的分享!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2015-9-28 08:57:11 | 显示全部楼层
经验贴值得借鉴
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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