MeteoInfo脚本示例：多站点数据平均

MeteoInfo

日常用到的站点数据（比如micaps 1类地面全要素观测）每个时次包含的站点都不一样，要同时计算多个站点数据的平均值就比较困难。MeteoInfoLab提供了Series类可以对一维数组设置index，方便此类操作。这里给出一个计算micaps 1类数据多个站点日平均温度的例子供参考。需要首先确定计算哪些站点（在文件cma_stations_base.csv中），然后在提取micaps文件中站点和温度数组，组合成为Series数据类型，再从中提取出设定站点数据（如果某个设定站点在micaps数据中不存在则为nan），放入一个数据列表中，最后对数据列表进行平均即可。

1. #Get stations
   
2. stfn = os.path.join('D:/Working/Stations/cma_stations_base.csv')
   
3. table = readtable(stfn, delimiter=',', format='%4s%3f')
   
4. stids = table[u'区站号']
   
5. lons = table[u'经度']
   
6. lats = table[u'纬度']
   
7. 
   
8. #Read data
   
9. datadir = r'U:\data\micaps\2016\plot'
   
10. year = 2016
    
11. month = 12
    
12. day = 1
    
13. t = datetime.datetime(year, month, day)
    
14. temps = []
    
15. for hour in range(2, 24, 3):
    
16.     t = t.replace(hour=hour)
    
17.     fn = os.path.join(datadir, t.strftime('%y%m%d%H') + '.000')
    
18.     print fn
    
19.     f = addfile_micaps(fn)
    
20.     st = f1['Stid'][:]
    
21.     temp = f1['Temperature'][:]   
    
22.     #Make series data
    
23.     s = series(temp, index=st)
    
24.     #Extract data by stations
    
25.     s = s[stids]
    
26.     temps.append(s.data)
    
27. 
    
28. #Average
    
29. ave = mean(temps)
    
30. 
    
31. #Plot
    
32. axesm(bgcolor=(204,255,255))
    
33. lworld = shaperead('D:/Temp/map/country1.shp')
    
34. lchina = shaperead('D:/Temp/map/bou2_4p.shp')
    
35. geoshow(lworld, edgecolor='k', facecolor=(255,251,195))
    
36. geoshow(lchina, edgecolor='k')
    
37. levs = arange(0, 35, 2)
    
38. layer = scatterm(lons, lats, ave, levs)
    
39. colorbar(layer)
    
40. yticks([20,30,40,50])
    
41. title('Daily averged temperature (' + t.strftime('%Y-%m-%d') + ')')
    
42. xlim(72, 136)
    
43. ylim(16, 55)

#Get stations<br /> stfn = os.path.join('D:/Working/Stations/cma_stations_base.csv')<br /> table = readtable(stfn, delimiter=',', format='%4s%3f')<br /> stids = table[u'区站号']<br /> lons = table[u'经度']<br /> lats = table[u'纬度']<br /> <br /> #Read data<br /> datadir = r'U:\data\micaps\2016\plot'<br /> year = 2016<br /> month = 12<br /> day = 1<br /> t = datetime.datetime(year, month, day)<br /> temps = []<br /> for hour in range(2, 24, 3):<br />     t = t.replace(hour=hour)<br />     fn = os.path.join(datadir, t.strftime('%y%m%d%H') + '.000')<br />     print fn<br />     f = addfile_micaps(fn)<br />     st = f1['Stid'][:]<br />     temp = f1['Temperature'][:]    <br />     #Make series data<br />     s = series(temp, index=st)<br />     #Extract data by stations<br />     s = s[stids]<br />     temps.append(s.data)<br /> <br /> #Average<br /> ave = mean(temps)<br /> <br /> #Plot<br /> axesm(bgcolor=(204,255,255))<br /> lworld = shaperead('D:/Temp/map/country1.shp')<br /> lchina = shaperead('D:/Temp/map/bou2_4p.shp')<br /> geoshow(lworld, edgecolor='k', facecolor=(255,251,195))<br /> geoshow(lchina, edgecolor='k')<br /> levs = arange(0, 35, 2)<br /> layer = scatterm(lons, lats, ave, levs)<br /> colorbar(layer)<br /> yticks([20,30,40,50])<br /> title('Daily averged temperature (' + t.strftime('%Y-%m-%d') + ')')<br /> xlim(72, 136)<br /> ylim(16, 55)

zihan6688

MeteoInfo脚本很强大

ddd212

学习学习。。。。

liuxiunuli

谢谢老师分享，学习中

guoaxing

谢谢老师分享，学习中

lynpatty

王老师请问一下有没有降水阈值相关的脚本呢？

MeteoInfo

lynpatty 发表于 2018-6-9 12:40
王老师请问一下有没有降水阈值相关的脚本呢？

获取降水数组，用percentile函数得到降水阈值。比如降水量数组是 a ，95%的阈值是 v = percentile(a, 95)

mikado

#Get stations
stfn = os.path.join('D:/Working/Stations/cma_stations_base.csv')
table = readtable(stfn, delimiter=',', format='%4s%3f')
stids = table[u'区站号']
lons = table[u'经度']
lats = table[u'纬度']

王老师可以展示一下csv文件的内容吗？
看起来有三列数据，但是format='%4s%3f'看起来是定义了两类格式，
而且区站号一般是5位字符，这里的%4s不知道怎么理解

最近在试用readtable读取包含“区站号 经度 纬度 站名”的文本文件，看了家园的很多样例和官网的介绍，还是不太明白这些参数的作用，主要是format和usecols

MeteoInfo

mikado 发表于 2020-7-24 02:24
#Get stations
stfn = os.path.join('D:/Working/Stations/cma_stations_base.csv')
table = readtable(s ...

format='%4s%3f'值的数据前4列是字符型，后续3列是浮点型。