grads批量合并有规律的二进制文件（如NC逐日资料等）和批量提取特定时间的特定要素

river

grads批量合并有规律的二进制文件（如NC逐日资料等）和批量提取特定时间的特定要素

最近毕业季啊，开始各种论文写作啊，然后突然发现自己各种不会啊，有木有？！再然后论坛里各种奇葩问题，有木有？！很多热心坛子看帖帮忙回答的头晕目眩，有木有？！算了，不吐槽了。


本来不打算发这个帖子的，但是最近发现还是有很多人问这样的问题，比如求助：关于批处理nc数据的问题；比如请问如何将多个单时次NC资料合成一个多时次的数据?还有很多形如“求助求助”、“小白求助”之类的让人不明就里的帖子也都问到了相关问题。


不想发这个帖子，不是因为问这个问题的人少；也不是因为我自己忙没时间；更不会是因为觉得这种问题不值得回答。而是因为这种方法并不是处理这种问题最好的方法，但看到还是很多人问啊问的，那就先发上来让大家解解燃眉之急吧。


一 ，给出一个例子，使用set fwrite 命令将多个文件名连续，内部时间连续的NC逐日资料合并成一个二进制文件，文件名为hgt.1948.nc  hgt.1949.nc  hgt.1950.nc·······（对于ncep资料的熟悉的人一看就明白了，就不详细说明文件格式什么的了，如果您还不是很明白，那就搜索一下ncep资料的介绍看一下）并配上ctl文件画图。


首先给出的是合并多个文件的gs

1. 'reinit'
   
2. 'set gxout fwrite'
   
3. 'set fwrite f:/1948-1952.grd'
   
4. yy=1948
   
5. while(yy<=1952)
   
6. 'sdfopen f:/ncep/daily/hgt.'yy'.nc'
   
7. 'set x 1 144'
   
8. 'set y 1 73'
   
9. if ((math_mod(yy,4)=0&math_mod(yy,100)!=0)|math_mod(yy,400)=0)
   
10. tt=1
    
11. while(tt<=366)
    
12. 'set t 'tt''
    
13. zz=1
    
14. while(zz<=17)
    
15. 'set z 'zz''
    
16. 'd hgt'
    
17. zz=zz+1
    
18. endwhile
    
19. tt=tt+1
    
20. endwhile
    
21. else
    
22. tt=1
    
23. while(tt<=365)
    
24. 'set t 'tt''
    
25. zz=1
    
26. while(zz<=17)
    
27. 'set z 'zz''
    
28. 'd hgt'
    
29. zz=zz+1
    
30. endwhile
    
31. tt=tt+1
    
32. endwhile
    
33. endif
    
34. 'close 1'
    
35. yy=yy+1
    
36. endwhile
    
37. 'disable fwrite'
    
38. ;

'reinit'<br /> 'set gxout fwrite'<br /> 'set fwrite f:/1948-1952.grd'<br /> yy=1948<br /> while(yy<=1952)<br /> 'sdfopen f:/ncep/daily/hgt.'yy'.nc'<br /> 'set x 1 144'<br /> 'set y 1 73'<br /> if ((math_mod(yy,4)=0&math_mod(yy,100)!=0)|math_mod(yy,400)=0)<br /> tt=1<br /> while(tt<=366)<br /> 'set t 'tt''<br /> zz=1<br /> while(zz<=17)<br /> 'set z 'zz''<br /> 'd hgt'<br /> zz=zz+1<br /> endwhile<br /> tt=tt+1<br /> endwhile<br /> else<br /> tt=1<br /> while(tt<=365)<br /> 'set t 'tt''<br /> zz=1<br /> while(zz<=17)<br /> 'set z 'zz''<br /> 'd hgt'<br /> zz=zz+1<br /> endwhile<br /> tt=tt+1<br /> endwhile<br /> endif<br /> 'close 1'<br /> yy=yy+1<br /> endwhile<br /> 'disable fwrite'<br /> ;

其中需要注意的是'set gxout fwrite'  'set fwrite f:/1948-1952.grd'的位置，要合并成一个文件所以这个命令最好放在循环外面； 'close 1'，是因为GrADS只能同时打开20个文件左右，如果要合并的资料很多话就需要打开一个转换一个关掉一个，否则提取出来的资料是前面资料的重复或者直接报错。


提取完成，生成1948-1952.grd文件，配如下ctl：

1. dset F:/1948-1952.grd
   
2. title mean daily NMC Reanalysis （1948-1952）
   
3. undef -999000000.000000
   
4. xdef 144 linear 0 2.5
   
5. ydef  73 linear -90 2.5
   
6. zdef 17 levels 1000 925 850 700 600 500 400 300 250 200 150 100 70 50 30 20 10
   
7. tdef 1827 linear 00Z01Jan1948 1440mn
   
8. vars 1
   
9. hgt 17 99  mean Daily Geopotential height
   
10. endvars

复制代码

对不对啊这样？画图对比一下就知道了，当然是一样的了，不然我也不拿出来了。


二，开始第二个问题，从多个nc文件中提取特定时间特定要素的资料并分别存放在不同的grd中，也给出例子，所使用的资料同上。
若想要提取从1948-1952年nc逐日资料文件中6-8月的位势高度资料

1. 'reinit'
   
2. 
   
3. yy=1948
   
4. while(yy<=1952)
   
5. 'set gxout fwrite'
   
6. 'set fwrite f:/'yy'.grd'
   
7. 'sdfopen f:/ncep/daily/hgt.'yy'.nc'
   
8. 
   
9. if ((math_mod(yy,4)=0&math_mod(yy,100)!=0)|math_mod(yy,400)=0)
   
10. tt=153
    
11. while(tt<=244)
    
12. 'set t 'tt''
    
13. zz=1
    
14. while(zz<=17)
    
15. 'set z 'zz''
    
16. 'set x 1 144'
    
17. 'set y 1 73'
    
18. 'd hgt'
    
19. zz=zz+1
    
20. endwhile
    
21. tt=tt+1
    
22. endwhile
    
23. else
    
24. tt=152
    
25. while(tt<=243)
    
26. 'set t 'tt''
    
27. zz=1
    
28. while(zz<=17)
    
29. 'set z 'zz''
    
30. 'set x 1 144'
    
31. 'set y 1 73'
    
32. 'd hgt'
    
33. zz=zz+1
    
34. endwhile
    
35. tt=tt+1
    
36. endwhile
    
37. endif
    
38. 'close 1'
    
39. yy=yy+1
    
40. endwhile
    
41. 'disable fwrite'
    
42. ;

'reinit'<br /> <br /> yy=1948<br /> while(yy<=1952)<br /> 'set gxout fwrite'<br /> 'set fwrite f:/'yy'.grd'<br /> 'sdfopen f:/ncep/daily/hgt.'yy'.nc'<br /> <br /> if ((math_mod(yy,4)=0&math_mod(yy,100)!=0)|math_mod(yy,400)=0)<br /> tt=153<br /> while(tt<=244)<br /> 'set t 'tt''<br /> zz=1<br /> while(zz<=17)<br /> 'set z 'zz''<br /> 'set x 1 144'<br /> 'set y 1 73'<br /> 'd hgt'<br /> zz=zz+1<br /> endwhile<br /> tt=tt+1<br /> endwhile<br /> else<br /> tt=152<br /> while(tt<=243)<br /> 'set t 'tt''<br /> zz=1<br /> while(zz<=17)<br /> 'set z 'zz''<br /> 'set x 1 144'<br /> 'set y 1 73'<br /> 'd hgt'<br /> zz=zz+1<br /> endwhile<br /> tt=tt+1<br /> endwhile<br /> endif<br /> 'close 1'<br /> yy=yy+1<br /> endwhile<br /> 'disable fwrite'<br /> ;

同样要注意'set gxout fwrite' 、'set fwrite f:/'yy'.grd'的位置

这个我没有具体的去验证，因为看网上有类似的帖子，而且这和上面的合并是类似的。有需要的同学自己去验证吧。

最后还要强调一点，这个方法并不是处理上述问题最理想的方法。文件多的时候运行会比较慢，而且gs这么长，用到的循环也比较多，很容易出错。其实最好的方法处理上述问题就是使用一个ctl批量描述这些二进制文件，这和合并资料是一样的，而且省时省力，运行又快。所以不推荐以上方法。

ctl批量描述的帖子可以在这里看 nc资料的描述文件以及批量描述nc资料的方法。当然还有很多其他坛友发的相关帖子，大家可以自己在论坛搜索


就这么多，如果发现问题，还请不吝赐教！

river

jeffggg 发表于 2013-5-21 22:42

                               
登录/注册后可看大图

一个ctl描述2个NC   画出图有问题 大神帮忙看看

缺测值有问题吧，应该是32766，改下试试。使用fwrite命令提取出来的nc资料的缺测值才是-9.99E+33。

曾记否

顶一个，楼主辛苦了。建议楼主把gs分一下行

river

曾记否 发表于 2013-5-16 10:22

                               
登录/注册后可看大图

顶一个，楼主辛苦了。建议楼主把gs分一下行

速度很快啊，沙发被你抢了。添加代码那个命令太让人郁闷，还得自己排版。下回直接贴上去好了

kongfeng0824

顶一个！！！！！！！！！！

天道酬勤

楼主威武！

夏朗的芒果

LZ威武。

mofangbao

river 发表于 2013-5-16 10:35

                               
登录/注册后可看大图

速度很快啊，沙发被你抢了。添加代码那个命令太让人郁闷，还得自己排版。下回直接贴上去好了

你是不是用了GsEditor之类的编辑器啊，那个里面的文字都是带格式的，你可以复制出来，贴到记事本里面“过滤”一下，再复制出来，然后在插入代码的地方插入就行了。
这个帖子真不错，总结的挺到位啊！

兰溪之水

建议写ctl描述nc或二进制文件~

river

mofangbao 发表于 2013-5-16 12:06

                               
登录/注册后可看大图

你是不是用了GsEditor之类的编辑器啊，那个里面的文字都是带格式的，你可以复制出来，贴到记事本里面“过 ...

我直接从写字板里复制出来的，插入代码的时候估计粘贴成纯文本应该就好了。

river

兰溪之水 发表于 2013-5-16 12:15

                               
登录/注册后可看大图

建议写ctl描述nc或二进制文件~

呜呜好激动，兰溪哥哥第一次回我的贴啊
对的，我帖子里就说不推荐上面的方法，最好的还是ctl批量描述