MeteoInfoLab脚本示例：纬向平均的偏差

MeteoInfo

这个例子演示如何计算水平格点数据的纬向平均偏差，参考了NCL网站上的例子：http://www.ncl.ucar.edu/Applications/dev.shtml 。mean函数可以做不同维的平均。程序中ts是水平二维数组（经纬度格点场，64*128），zm = mean(ts, axis=1) 获得纬向平均一维数组（64），ts和zm的维数不同，直接相减必须符合数组的broadcasting规则，也就是从后往前，两个数组的维相同或者其中一个为1。显然目前是不符合的，可以对zm reshape，获得二维数组（64*1），这样ts和zm数组就符合broadcasting规则了，可以直接相减。也可以用broadcast_to函数将zm扩展为和ts一样的shape，然后再相减。

关于数组broadcasting的内容可以参考这里：https://docs.scipy.org/doc/numpy/user/basics.broadcasting.html 。

1. fn = 'D:/Temp/nc/83.nc'f = addfile(fn)
   
2. ts = f['TS'][0,:,:]
   
3. 
   
4. #calculate deviation from zonal mean
   
5. zm1 = mean(ts, axis=1)
   
6. zm = zm1.reshape(len(zm1), 1)
   
7. #zm = broadcast_to(zm, (ts.shape[0], ts.shape[1]))
   
8. nts = ts - zm
   
9. 
   
10. #Plot
    
11. #Plot origin data
    
12. lworld = shaperead('D:/temp/map/country1.shp')
    
13. axesm(position=[0.05, 0.55, 0.8, 0.4])
    
14. geoshow(lworld, edgecolor='b')
    
15. levs = arange(240, 301, 5)
    
16. layer = contourm(ts, levs, color='k')
    
17. ylim(-90, 90)
    
18. clabel(layer)
    
19. title('Origin data')
    
20. 
    
21. #Plot zonal mean
    
22. axes(position=[0.8, 0.55, 0.15, 0.4])
    
23. plot(zm1, ts.dimvalue(0))
    
24. yaxis(tickvisible=False)
    
25. xlim(200, 320)
    
26. ylim(-90, 90)
    
27. 
    
28. #Plot deviation data
    
29. axesm(position=[0.05, 0.05, 0.8, 0.4])
    
30. geoshow(lworld.clone(), edecolor='k')
    
31. levs = arange(-35, 36, 5)
    
32. layer = contourfm(nts, levs, edgecolor='gray')
    
33. ylim(-90, 90)
    
34. colorbar(layer)
    
35. title('Deviation from zonal average')

fn = 'D:/Temp/nc/83.nc'f = addfile(fn)<br /> ts = f['TS'][0,:,:]<br /> <br /> #calculate deviation from zonal mean<br /> zm1 = mean(ts, axis=1)<br /> zm = zm1.reshape(len(zm1), 1)<br /> #zm = broadcast_to(zm, (ts.shape[0], ts.shape[1]))<br /> nts = ts - zm<br /> <br /> #Plot<br /> #Plot origin data<br /> lworld = shaperead('D:/temp/map/country1.shp')<br /> axesm(position=[0.05, 0.55, 0.8, 0.4])<br /> geoshow(lworld, edgecolor='b')<br /> levs = arange(240, 301, 5)<br /> layer = contourm(ts, levs, color='k')<br /> ylim(-90, 90)<br /> clabel(layer)<br /> title('Origin data')<br /> <br /> #Plot zonal mean<br /> axes(position=[0.8, 0.55, 0.15, 0.4])<br /> plot(zm1, ts.dimvalue(0))<br /> yaxis(tickvisible=False)<br /> xlim(200, 320)<br /> ylim(-90, 90)<br /> <br /> #Plot deviation data<br /> axesm(position=[0.05, 0.05, 0.8, 0.4])<br /> geoshow(lworld.clone(), edecolor='k')<br /> levs = arange(-35, 36, 5)<br /> layer = contourfm(nts, levs, edgecolor='gray')<br /> ylim(-90, 90)<br /> colorbar(layer)<br /> title('Deviation from zonal average')

weinihou

来给大神点赞~~~~膜拜膜拜~~~~~
{:eb502:}{:eb502:}