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

游子

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也不必急于丢掉。
当然了，在这里班门弄斧了。

/猪头哼哼

只是为了画图快看个结果的话
不如试试meteoinfo和Panoply 后者可以打开几乎一切NC文件并出图。。前者也能打开部分grads打不开的nc 画图自然是很轻松自在。
也省去了楼主NCL转换的复杂过程。

游子

/猪头哼哼 发表于 2013-10-21 21:28

                               
登录/注册后可看大图

只是为了画图快看个结果的话
不如试试meteoinfo和Panoply 后者可以打开几乎一切NC文件并出图。。 ...

再简单一点直接用更懒惰的软件 ： ncBrowse.
转换的过程也是思考的过程喔，为了出图而出图就没劲了

ayzqs

顶！有时候处理这些资料挺烦的

/猪头哼哼

游子 发表于 2013-10-21 21:39

                               
登录/注册后可看大图

再简单一点直接用更懒惰的软件 ： ncBrowse.
转换的过程也是思考的过程喔，为了出图而出图就没劲了

那个和panoply差不多懒惰啦=。= 楼主的方法自然是很赞的，我也用其类似的方法用NCL转换过。。。
楼主的目的也就是为了简单的出图看看实验结果，免去NCL繁琐的设置不是吗~那么用懒惰的方法出图看看也挺好的~

mofangbao

谢谢楼主前来分享~发挥各种工具的优势~

游子

mofangbao 发表于 2013-10-22 08:55

                               
登录/注册后可看大图

谢谢楼主前来分享~发挥各种工具的优势~

清风大哥过奖了。

走走

经验值得借鉴

singapple

我是菜鸟，正在奇怪为什么会打不开自己的.nc格式的文件呢，原来是这个原因啊，我的也是cesm出的运行结果~学习了~~技多不压身，能多了解一些软件也是好的~~谢谢楼主的分享，和楼上的分享！

huangwei

经验贴值得借鉴