MeteoInfoLab脚本示例：计算不同区域平均值

MeteoInfo

这里用美国做例子，有一个美国区域的格点温度场数据（usgrid.data），需要计算出每个州（state）的平均温度。当然需要有一个包含各州行政区域的shape文件了（相关文件可以在此帖中下载：http://bbs.06climate.com/forum.p ... d=11070&fromuid=106）。首先读取格点数据数组，然后读取states.shp文件生成us图层，对us图层所有的图元（shape，每个图元就是一个州）遍历，对于每个州，利用该州的图元maksout格点数据，然后求平均值、最大值和最小值，并打印出来（当然也可以输出到文件中，具体在网上找python输出文本文件的例子）。后面的绘图语句只是为了看看数据的分布情况。

脚本程序：

1. #Add a surfer grid data
   
2. f = addfile_surfer('D:/Temp/ascii/usgrid.dat')
   
3. tdata = f['var'][:,:]
   
4. #Read US shape file
   
5. us = shaperead('D:/Temp/map/states.shp')
   
6. #Average temporature for each state
   
7. i = 0
   
8. for rpoly in us.shapes():
   
9.     name = us.cellvalue('STATE_NAME', i)
   
10.     mdata = tdata.maskout(rpoly)
    
11.     tave = mdata.ave()
    
12.     tmin = mdata.min()
    
13.     tmax = mdata.max()
    
14.     print name + ', Ave: %.2f, Min: %.2f, Max: %.2f' %(tave, tmin, tmax)
    
15.     i += 1
    
16. 
    
17. #Plot
    
18. axesm()
    
19. geoshow('country')
    
20. geoshow(us, edgecolor=[0,0,255])
    
21. layer = contourfm(tdata,20)
    
22. title('Temporature distribution map')
    
23. colorbar(layer)

#Add a surfer grid data <br /> f = addfile_surfer('D:/Temp/ascii/usgrid.dat') <br /> tdata = f['var'][:,:] <br /> #Read US shape file <br /> us = shaperead('D:/Temp/map/states.shp') <br /> #Average temporature for each state <br /> i = 0 <br /> for rpoly in us.shapes(): <br />     name = us.cellvalue('STATE_NAME', i) <br />     mdata = tdata.maskout(rpoly) <br />     tave = mdata.ave() <br />     tmin = mdata.min() <br />     tmax = mdata.max() <br />     print name + ', Ave: %.2f, Min: %.2f, Max: %.2f' %(tave, tmin, tmax) <br />     i += 1 <br /> <br /> #Plot <br /> axesm() <br /> geoshow('country') <br /> geoshow(us, edgecolor=[0,0,255]) <br /> layer = contourfm(tdata,20) <br /> title('Temporature distribution map') <br /> colorbar(layer)

随缘

真是及时雨啊，前两天刚看到有人提问呢

killerlrving

顶~~~~~~~~~~~~~~~~~~~~~~~~~~~~

zyfws

从grads读取的数据场怎么做指定时次内的平均，比如一个72小时内的逐小时数据场，需要做24小时内的平均该怎么弄，我怎么都提示维数不对或没有ave这个方法，比如：
getV=f['ps']
ave=getV[0:23,:,:].ave()这样就不行，请问王老师该怎么弄？

MeteoInfo

zyfws 发表于 2016-2-17 11:52
从grads读取的数据场怎么做指定时次内的平均，比如一个72小时内的逐小时数据场，需要做24小时内的平均该怎 ...

参考此贴：MeteoInfoLab脚本示例：沿特定维平均
http://bbs.06climate.com/forum.p ... 039&fromuid=106
(出处: 气象家园)

磊磊

太牛掰了  正好需要

aas314

大赞 学习了

liuxiunuli

真是及时雨，大赞 学习了

hxbbing

顶  学习学习

vortexroc

王老师，这个示例是针对格点数据的，有没有针对站点数据的统计方法，比如把全国地面气象站的降雨量进行分省的统计，得到逐小时的各省降雨量合计？