grads处理多个ctl文件和nc文件

grads处理多个ctl文件和nc文件2011-10-1021:03:59|分类：grads学习|标签：|举报|字号大中小订阅下载LOFTER我的照片书|用grads处理多个相同格式的数据时若单个单个处理非常麻烦，当文件非常多的时候是单个处理是不实际的。下面介绍一种方法;第一步，在这种情况下可以重新写一个ctl描述文件，其文件变量都和已知的ctl相同，若原来的n文件只是时间不同，那么新描述文件的时间维数是所有原文件的时间的和。同样，若其他维数不同时也用同样的方法处理。第二步，在第一行之后添加一行：optionstemplate表示多个时间序列原始数据文件想用一个描述文件统一地描述。这些原数据的原文件名由dset定义的形势命名文件名。第三步，修改dset的文件名。原路径不变，把文件名用％表示。其中：％y2代表两位数年％y4代表四位数年％m1代表一位或者两位数的月％m2代表两位数月（用0补齐1位数）％mc3个字符月份的缩写％d11或2位天％d2两位天％h11或者2位时％h22位时例如：原文件其中之一的文件名为gdas2006050812f00，且所有文件只有天和时的变化那么新描述文件的文件名为：gdas200605%d2%h2f00另外如果源文件里有index项的话，需要修改其idx的文件名，假设改成fnl.idx。并用在dos下用gribmap函数生成一个新的idx文件。gribmap-e-ifnl.ctl（加绝对路径）openfnl.ctl就可以打开所有文件。**********************************************************************************************************************************若想要提取从1951-2006年56年nc文件中的某些数据，一个一个处理非常麻烦，这里介绍种较为简易的方法。例如想提取6-8月的位势高度资料。'reinit't5=1951*作文件名循环while(t5=2006)'setgxoutfwrite''setfwriteD:\sichuan\hgt1\'%t5%'.dat''sdfopene:\ncep1\hgt\hgt.'%t5%'.nc't3=t5-1950*判断是否为闰年if（t3=2|t3=6|t3=10|t3=14|t3=18|t3=22|t3=26|t3=30|t3=34|t3=38|t3=42|t3=46|t3=50|t3=54)to=153elseto=152endift4=to+91while(to=t4)'sett'tot1=1while(t1=12)'setz't1'setlon80140''setlat1555''dhgt't1=t1+1endwhileto=to+1endwhile*这里必须先观点上述运行的文件，grads最多同时可以打开20个文件左右。'reinit't5=t5+1endwhile'reinit'这样可以提取你想要的年数据，然后你大可运用fortran对数据进行随心所欲的处理。能否直接生成一个文件还正在探索中。=viewthread&tid=7310&extra=&page=1批量读取nc数据，用你的方法成功了，谢谢！！！直接配个批量描述的ctl就可以了有一批nc数据，一个月一个文件，现将文件名改为：197901.nc，197902.nc，依次类推，对二进制的数据知道写ctl文件来进行批处理运算，那么nc数据应该怎么做呢？试过了写ctl文件，sdfopen***\%y4%m2.nc，year=1978while(year=2011)month=01while(month=12)'sdfopen***\'year''month'.nc'...month=month+1endwhileyear=year+1endwhile实我也是糊里糊涂的解决了。。。ctl文件如下：dset^%y4%m2.ncundef1e+15optionstemplatetitleMERRAdatadtypeNetCDFydef144linear-901.25xdef288linear-1801.25zdef21levels1000975950925900875850825800775750725700650600550500450400350300tdef396linear00Z01JAN19791movars3qv21t,z,y,xSpecifichumidityu21t,z,y,xEastwardwindcomponentv21t,z,y,xNorthwardwindcomponentendvars然后open***.ctl就行了，之前的问题是打不开ctl文件，怎么改也不行，后来换了台机子就好了。。。所以我说我也不知道怎么回事，希望对你有帮助。我之前就这样做的，能打开ctl文件，但是d之后，都显示allundefinedvalues，我的ctl如下，麻烦帮我看看哪里错了？dset^%y4%m2.ncundef1e+15optionstemplatetitleMERRAdatadtypenetcdfxdef288linear-179.3751.25ydef144linear-89.3751.25zdef21levels1000975950925900875850825800775750725700650600550500450400350300tdef396linear00Z01JAN19791mnvars3qv21t,z,y,xSpecifichumidityu21t,z,y,xEastwardwindcomponentv21t,z,y,xNorthwardwindcomponentendvars可以的，但文件名一定要连续这个时间长度一定要和你的文件对应好最近发了一个利用grads批量合并nc文件的帖子=viewthread&tid=14459，因为那不是最好的方法，所以非常推荐了利用ctl批量描述nc的方法。很多人也都问怎么批量描述，，很多人也都去自己尝试，尝试过程中出了各种问题。推荐兰溪给nc文件写ctl的帖子=viewthread&tid=6008，尤其是使用sdfopen后出现“SDFfilehasnodiscernableXcoordinate”的同学们，非常推荐这个帖子。肯定很多人已经看过这个帖子，也照着做过，还是出了问题。其实论坛还有很多其他有关ctl批量描述nc的帖子，大家参考一下，对照自己出现的问题，应该大部分情况下都能解决。但是有一些人很着急，连一个nc文件的描述文件都写不对，就直接批量的，那肯定只会出更多的错误，还一时半会儿改不对。只会造成时间的浪费。我不是来说教的，只是我觉得什么事儿都是循序渐进的，不要那么浮躁。俗话说磨刀不误砍柴工，其实真的是这么回事。就说批量描述nc文件，它也是在正确描述一个nc文件的基础上来修改的。只要描述一个的出来了，往上加一两个语句，用正确的替代格式替代了文件名，再改一下时间长度就出来了。那能正确描述一个nc文件的ctl就是你磨得很锋利的砍柴刀，有了它后面的柴就相当好砍了。下面我就以ncep逐日再分析资料为例说一下nc文件的描述和批量描述。文件名连续，hgt.1948.nchgt.1949.nchgt.1950.nc`````````首先看一下它自带的ctl，如下图红色圈起来的部分，就不详细解释了，对grads有一定了解的人，看一下就知道什么下面就可以根据这个自己编写一个ctl文件来描述这个nc文件了。照着自带的写下来先，存成1948.ctl然后使用open命令打开，画图，发现出来的图和缺测值设置错误时差不多。那就想方设法的改缺测值，比如用在grads中使用qattr查看有一句missing_value32766，修改缺测值，再画图，还是那样；qundef出来的是-999000000.000000，再改，再画还是那样。无论怎么改，都还是那样。如果你这么反复折腾，最后还没发现问题，那就说明你没有好好利用论坛的搜索功能，也没有看到兰溪的帖子，也没看到黎大叶子的用Fortran批量为nc写ctl的帖子=viewthread&tid=7267。但是不用着急，我看见了，其中最重要的是打开netcdf格式数据的描述文件是需要用xdfopen命令的，那就先要去看看xdfopen能打开的ctl需要怎么写。看完了，差不多就能明白了，有这么多前人的成果，那就照着修改呗，最终我修改出来了，图画的也很正常了。写出来的如下1.dsetF:/ncep/daily/hgt.1948.nc2.titlemeandailyNMCReanalysis（1948）3.undef-9994.xdeflon144linear02.55.ydeflat73linear-902.56.zdeflevel17levels100092585070060050040030025020015010070503020107.tdeftime366linear00Z01Jan19481440mn8.vars19.hgt=hgt17t,z,y,xmeanDailyGeopotentialheight10.endvars复制代码光看着正常不行啊，需要和原始的图对比验证了才能确定是对的吧。所以在grads里面用sdfopen命令打开hgt.1948.nc画第一层第一个时次的图，再用xdfopen打开你编写的ctl，也画第一层第一个时次的图看看。我擦！！！！对不上啊，看来还是有问题啊。说明上面的ctl是有问题的，还得改进才行。其实到这步已经基本成功了，仔细看看叠加起来的那两张图，几乎是一样的，只是南北的方向是反的。那就好办了，在ctl里加上optionsyrev，告诉grads南北要反向不就行了么。于是最终的ctl出来了，如下1.dsetF:/ncep/daily/hgt.1948.nc2.titlemeandailyNMCReanalysis（1948）3.optionsyrev4.*yrev表示y轴反向5.undef-9996.xdeflon144linear02.57.ydeflat73linear-902.58.zdeflevel17levels100092585070060050040030025020015010070503020109.tdeftime366linear00Z01Jan19481440mn10.vars111.hgt=hgt17t,z,y,xmeanDailyGeopotentialheight12.endvars复制代码再用xdfopen打开画图，这回就一模一样了啊。说明成功了！那下面的批量描述就太简单了，比如我要批量描述1948和1949两年的，算一下一个闰年一个平年，一共有时次366+365=731，那么就修改ctl吧dsetF:/ncep/daily/hgt.%y4.nctitlemeandailyNMCReanalysisoptionstemplateoptionsyrev*yrev表示y轴反向*undef-999xdeflon144linear02.5ydeflat73linear-902.5zdeflevel17levels10009258507006005004003002502001501007050302010tdeftime731linear00Z01Jan19481440mnvars1hgt