grads将多个nc文件处理并写入一个grb文件，出图明显错误，说明grb里面的数据有错误...

LemonTree007

将三个时间连续的月平均pdsi数据nc文件处理成年平均数据并写入grb文件中。
gs 如下：

1. 'reinit'
   
2. 'set gxout fwrite'
   
3. 'set fwrite d:/matlabht/pdsi/pdsi1901-2014.grd'
   
4. 'sdfopen d:/matlabht/remap05x05.pdsi.mon.mean.selfcalibrated_190001-194912.nc'
   
5. 'sdfopen d:/matlabht/remap05x05.pdsi.mon.mean.selfcalibrated_195001-199912.nc'
   
6. 'sdfopen d:/matlabht/remap05x05.pdsi.mon.mean.selfcalibrated_200001-201412.nc'
   
7. 
   
8. 'set dfile 1'
   
9. tt=13
   
10. while(tt<600)
    
11.     'set lon 70 140'
    
12.     'set lat 10 60'
    
13.     'set lev 0'
    
14.     'set t 'tt''
    
15.     'd ave(pdsi,t='tt',t='tt'+11)'
    
16.     tt=tt+12
    
17. endwhile
    
18. 
    
19. 'set dfile 2'
    
20. tt=1
    
21. while(tt<600)
    
22.     'set lon 70 140'
    
23.     'set lat 10 60'
    
24.     'set lev 0'
    
25.     if(tt=1)
    
26.         'set t 'tt''
    
27.         'd ave(pdsi,t='tt',t='tt'+11)'
    
28.         tt=tt+12
    
29.     else if(tt>12)
    
30.         'set t 'tt''
    
31.         'd ave(pdsi,t='tt',t='tt'+11)'
    
32.         tt=tt+12
    
33.     endif
    
34. endwhile
    
35. 
    
36. 'set dfile 3'
    
37. tt=1
    
38. while(tt<180)
    
39.     'set lon 70 140'
    
40.     'set lat 10 60'
    
41.     'set lev 0'
    
42.     if(tt=1)
    
43.         'set t 'tt''
    
44.         'd ave(pdsi,t='tt',t='tt'+11)'
    
45.         tt=tt+12
    
46.     else if(tt>12)
    
47.         'set t 'tt''
    
48.         'd ave(pdsi,t='tt',t='tt'+11)'
    
49.         tt=tt+12
    
50.     endif
    
51. endwhile
    
52. 'disable fwrite'
    
53. ;
    

复制代码

运行之后成功写生成grb文件，ctl文件如下：
dset d:\matlabht\pdsi\pdsi1901-2014.grd
undef -99999
title  1991-2014 annual average sc_PDSI
xdef 141 linear 70 0.5
ydef 101 linear 10 0.5
zdef 1   linear 0 1
tdef 114 linear 00Z01JAN1901 1yr
vars  1
pdsi  0 99  Annual Self-calibrated Palmer Drought Severity Index
endvars

打开grads，open .ctl文件顺利，说明描述文件是正确的。
编写如下脚本用ctl文件画图，counterline数值明显有误，分析应该是grb生成的数据有误，正常的数据范围应该是-4—4左右。

1. 'reinit'
   
2. 'open D:\matlabht\pdsi\pdsi1901-2014.ctl'
   
3. 'set mpdset cnbasemap'
   
4. 'set lon 70 140'
   
5. 'set lat 10 60'
   
6. 'set lev 0'
   
7. 'set grads off'
   
8. 'set gxout shaded'
   
9. 'set grid off'
   
10. 'cnbasemap ave(pdsi,t=1,t=114)'
    
11. 'cbarn'
    
12. 'run southsea.gs'
    
13. 'printim D:\matlabht\tp\pdsi000.png white'
    
14. pull dummy
    
15. 'c'
    

复制代码

LemonTree007

错误已经找到，是ctl的问题，undef -99999改为-9.99e+8

旺旺小小酥

好帖子，学习到了！

LemonTree007

旺旺小小酥 发表于 2018-4-26 12:42
好帖子，学习到了！

哈哈，对我来说也是学习的过程，是在逛别人帖子的时候发现的。