登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
x
本帖最后由 river 于 2017-4-2 22:11 编辑
最近发了一个利用grads批量合并nc文件的帖子http://bbs.06climate.com/forum.php?mod=viewthread&tid=14459，因为那不是最好的方法，所以非常推荐了利用ctl批量描述nc的方法。 很多人也都问怎么批量描述，，很多人也都去自己尝试，尝试过程中出了各种问题。推荐兰溪给nc文件写ctl的帖子http://bbs.06climate.com/forum.php?mod=viewthread&tid=6008，尤其是使用sdfopen后出现“SDF file has no discernable X coordinate”的同学们，非常推荐这个帖子。肯定很多人已经看过这个帖子，也照着做过，还是出了问题。
其实论坛还有很多其他有关ctl批量描述nc的帖子，大家参考一下，对照自己出现的问题，应该大部分情况下都能解决。但是有一些人很着急，连一个nc文件的描述文件都写不对，就直接批量的，那肯定只会出更多的错误，还一时半会儿改不对。只会造成时间的浪费。 我不是来说教的，只是我觉得什么事儿都是循序渐进的，不要那么浮躁。俗话说磨刀不误砍柴工，其实真的是这么回事。就说批量描述nc文件，它也是在正确描述一个nc文件的基础上来修改的。只要描述一个的出来了，往上加一两个语句，用正确的替代格式替代了文件名，再改一下时间长度就出来了。那能正确描述一个nc文件的ctl就是你磨得很锋利的砍柴刀，有了它后面的柴就相当好砍了。
下面我就以ncep逐日再分析资料为例说一下nc文件的描述和批量描述。文件名连续，hgt.1948.nc
hgt.1949.nc hgt.1950.nc`````````
首先看一下它自带的ctl，如下图红色圈起来的部分，就不详细解释了，对grads有一定了解的人，看一下就知道什么意思了。
下面就可以根据这个自己编写一个ctl文件来描述这个nc文件了。
照着自带的写下来先，存成1948.ctl然后使用open命令打开，画图，发现出来的图和缺测值设置错误时差不多。那就想方设法的改缺测值，比如用在grads中使用q attr查看有一句missing_value 32766，修改缺测值，再画图，还是那样；q undef出来的是-999000000.000000，再改，再画还是那样。无论怎么改，都还是那样。
如果你这么反复折腾，最后还没发现问题，那就说明你没有好好利用论坛的搜索功能，也没有看到兰溪的帖子，也没看到黎大叶子的用Fortran批量为nc写ctl的帖子
http://bbs.06climate.com/forum.php?mod=viewthread&tid=7267。
但是不用着急，我看见了，其中最重要的是打开netcdf格式数据的描述文件是需要用xdfopen命令的，那就先要去看看xdfopen能打开的ctl需要怎么写http://grads.iges.org/grads/gadoc/gadocindex.html。
看完了，差不多就能明白了，有这么多前人的成果，那就照着修改呗，最终我修改出来了，图画的也很正常了。写出来的如下
- dset F:/ncep/daily/hgt.1948.nc
- title mean daily NMC Reanalysis （1948）
- undef -999
- xdef lon 144 linear 0 2.5
- ydef lat 73 linear -90 2.5
- zdef level 17 levels 1000 925 850 700 600 500 400 300 250 200 150 100 70 50 30 20 10
- tdef time 366 linear 00Z01Jan1948 1440mn
- vars 1
- hgt=>hgt 17 t,z,y,x mean Daily Geopotential height
- endvars
其实到这步已经基本成功了，仔细看看叠加起来的那两张图，几乎是一样的，只是南北的方向是反的。那就好办了，在ctl里加上options yrev，告诉grads南北要反向不就行了么。于是最终的ctl出来了，如下
- dset F:/ncep/daily/hgt.1948.nc
- title mean daily NMC Reanalysis （1948）
- options yrev
- * yrev表示y轴反向
- undef -999
- xdef lon 144 linear 0 2.5
- ydef lat 73 linear -90 2.5
- zdef level 17 levels 1000 925 850 700 600 500 400 300 250 200 150 100 70 50 30 20 10
- tdef time 366 linear 00Z01Jan1948 1440mn
- vars 1
- hgt=>hgt 17 t,z,y,x mean Daily Geopotential height
- endvars
再用xdfopen打开画图，这回就一模一样了啊。 说明成功了！
那下面的批量描述就太简单了，比如我要批量描述1948和1949两年的，算一下一个闰年一个平年，一共有时次366+365=731，那么就修改ctl吧
dset F:/ncep/daily/hgt.%y4.nc
title mean daily NMC Reanalysis
options template
options yrev
* yrev表示y轴反向
*undef -999
xdef lon 144 linear 0 2.5
ydef lat 73 linear -90 2.5
zdef level 17 levels 1000 925 850 700 600 500 400 300 250 200 150 100 70 50 30 20 10
tdef time 731 linear 00Z01Jan1948 1440mn
vars 1
hgt=>hgt 17 t,z,y,x mean Daily Geopotential height
endvars
就这么简单，只要用%y4代表四位年，把总的时次改成731，再增加一句options template就可以了。然后就可以利用xdfopen命令打开画图和原来的对比一下，一样一样的吧，可以批量描述了吧！
注意：有人可能注意到我把缺测那一项给注释掉了，其实这个是完全不需要的，去掉或者改成任何值，都不影响。这是问什么？不细说了，因为我没仔细研究（比较晚了，我要去睡觉了，以后再说，或者有人知道可以告诉我）
其他的资料都同理了，思路都差不多，变通一下就可以了。- 首先编写可以正确描述一个资料的ctl，并能正确出图。
- 修改该ctl里dest后的文件名，使用合适的替代格式替代（具体的格式看实用手册）
- 增加options template，这个是批量描述必须加上的一句
- 计算出所要批量描述的文件的总时次，修改ctl里的总时次
- 画图进行验证，如果不对，再根据具体情况做出相应的修改
- 大功告成，可以用了
- 需注意的一点就是options yrev，这是可选项，要根据资料实际情况来使用。
今天这个帖子有点儿罗嗦了，大家捡有用的看就行。其实说那么多，无非是想和一些人说我们这个年代的人已经是站在前人的肩膀上了，很多东西别人都已经有了经验，分享给你你就要好好利用，不要把时间都浪费在发帖和等回复上面。很多东西，只要自己有一定的基础，在一定合理的范围内就会找到那个答案的·····
| 
[复制链接]
窥视卡
雷达卡
发表于 2013-5-22 22:49:34
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
发表于 2013-5-22 23:44:10
发表于 2013-5-23 07:38:32
发表于 2013-5-23 09:30:57