有关用指数时间序列画直方图的问题

nuistzhou

mofangbao 发表于 2012-5-18 17:24

                               
登录/注册后可看大图

做出来是这样的么，两个变量的距平？

嗯，原理差不多。不过我这里要求的是一个变量40年来在2个区域的距平差作为一个指数。因此有40个距平差值。我下午把40年的夏季值提到了一个grd里，再处理，那么做区域平均的时候，还需不需要在前面加上  set t 1 40

mofangbao

先把距平值提取出来到一个单独的文件中
然后使用这个文件中的两个变量相减（如果不需要距平值，就直接写出距平相减的值就行了）来作图

1. 'reinit'
   
2. 'sdfopen K:\ncep\pressure\uwnd.mon.mean.nc '
   
3. 'set gxout fwrite'
   
4. 'set fwrite K:\ncep\uwnd.anomaly.grd'
   
5. i=282
   
6. 'set lev 850'
   
7. 
   
8. 
   
9. 'define sumave1=aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=41,y=45)'
   
10. 'define sumave2=aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=47,y=51)'
    
11. 
    
12. i=i+12
    
13. 
    
14. #40年夏季的区域平均
    
15. while(i<=750)
    
16.     'sumave1=sumave1+aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=41,y=45)'
    
17.     'sumave2=sumave2+aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=47,y=51)'
    
18.     i=i+12
    
19. endwhile
    
20. *all year-summer aave
    
21. 'sumave1 = sumave1/40.0'
    
22. 'sumave2=sumave2/40.0'
    
23. 
    
24. #每年夏季的区域平均
    
25. i=282
    
26. 'set lev 850'
    
27. 'set x 1'
    
28. 'set y 1'
    
29. 'set z 1'
    
30. 'define summer1=aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=41,y=45)'
    
31. 'define summer2=aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=47,y=51)'
    
32. 'd summer1'
    
33. 'd summer2'
    
34. i=i+12
    
35. while(i<=750)
    
36. #提取距平数据
    
37.     'summer1=sumave1-aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=41,y=45)'
    
38.     'summer2=sumave2-aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=47,y=51)'
    
39.     'd summer1'
    
40.     'd summer2'
    
41.     i=i+12
    
42. endwhile
    
43. 'disable fwrite'
    
44. ;
    

'reinit'<br /> 'sdfopen K:\ncep\pressure\uwnd.mon.mean.nc '<br /> 'set gxout fwrite'<br /> 'set fwrite K:\ncep\uwnd.anomaly.grd'<br /> i=282<br /> 'set lev 850'<br /> <br /> <br /> 'define sumave1=aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=41,y=45)'<br /> 'define sumave2=aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=47,y=51)'<br /> <br /> i=i+12<br /> <br /> #40年夏季的区域平均<br /> while(i<=750)<br />     'sumave1=sumave1+aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=41,y=45)'<br />     'sumave2=sumave2+aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=47,y=51)'<br />     i=i+12<br /> endwhile<br /> *all year-summer aave<br /> 'sumave1 = sumave1/40.0'<br /> 'sumave2=sumave2/40.0'<br /> <br /> #每年夏季的区域平均<br /> i=282<br /> 'set lev 850'<br /> 'set x 1'<br /> 'set y 1'<br /> 'set z 1'<br /> 'define summer1=aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=41,y=45)'<br /> 'define summer2=aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=47,y=51)'<br /> 'd summer1'<br /> 'd summer2'<br /> i=i+12<br /> while(i<=750)<br /> #提取距平数据<br />     'summer1=sumave1-aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=41,y=45)'<br />     'summer2=sumave2-aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=47,y=51)'<br />     'd summer1'<br />     'd summer2'<br />     i=i+12<br /> endwhile<br /> 'disable fwrite'<br /> ;<br />

输出的文件对应的ctl

1. dset K:\ncep\uwnd.anomaly.grd
   
2. title testdata
   
3. undef -9.99e+33
   
4. xdef 1 linear 1 1
   
5. ydef 1 linear 1 1
   
6. zdef 1 linear 1 1
   
7. tdef 40 linear jan1958 1yr
   
8. vars 2
   
9. an1 0 99 an1 data
   
10. an2 0 99 an2 data
    
11. endvars
    

dset K:\ncep\uwnd.anomaly.grd<br /> title testdata<br /> undef -9.99e+33<br /> xdef 1 linear 1 1<br /> ydef 1 linear 1 1<br /> zdef 1 linear 1 1<br /> tdef 40 linear jan1958 1yr<br /> vars 2<br /> an1 0 99 an1 data<br /> an2 0 99 an2 data<br /> endvars<br />

剩下的自己搞定吧，最基础的东西了，出来的图上面贴过了，直接是两个距平的数据做的

mofangbao

上面CTL中的开始年份乱写的，没去仔细算，你自己改一下

nuistzhou

mofangbao 发表于 2012-5-18 19:12

                               
登录/注册后可看大图

上面CTL中的开始年份乱写的，没去仔细算，你自己改一下

清风，在你的基础上小改了一下：
'reinit'
'sdfopen e:\uwnd.mon.mean.nc '
'set gxout fwrite'
'set fwrite e:\uwnd.grd'
i=282

'define s1=0'
'define s2=0'
while(i<=750)
    'set lev 850'
    'define sumave1=aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=41,y=45)'  
    'define sumave2=aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=47,y=51)'
   
's1=s1+sumave1'
's2=s2+sumave2'
  i=i+12
endwhile

's1=s1/40.0'
's2=s2/40.0'

i=282
'set lev 850'
while(i<=750)
    'juping1=aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=41,y=45)-s1'
    'juping2=aave(ave(uwnd,t='i',t='i'+2),x=41,x=61,y=47,y=51)-s2'
    'd juping1-juping2'
    i=i+12
endwhile
'disable fwrite'
;