MeteoInfoLab脚本示例：风玫瑰图中绘制污染物浓度

MeteoInfo

利用观测点风速、风向、污染物浓度的时间序列数据可以在风玫瑰图中绘制出污染物浓度等值线。脚本有点复杂，其实是用了两个不同的坐标系，等值线是用直角坐标系画的，风频曲线是用极坐标系画的。风玫瑰图的坐标系和传统的极坐标系也有区别，不同坐标系有一些转换的计算。

用到的示例数据文件： pm10.txt (618 Bytes, 下载次数: 800)

2017-6-30 16:11 上传

点击文件名下载附件

1. def windrose2polar(a):
   
2.     r = 360 - a + 90
   
3.     r[r>360] = r - 360
   
4.     return r
   
5. 
   
6. #Read data (wind speed, weed direction, pm10)
   
7. fn = os.path.join(migl.get_sample_folder(), 'ASCII', 'pm10.txt')
   
8. df = DataFrame.read_table(fn, format='%3f')
   
9. ws = df['WS'].values
   
10. wd = df['WD'].values
    
11. pm10 = df['PM10'].values
    
12. N = len(ws)
    
13. 
    
14. #Convert from windrose coordinate to polar coordinate
    
15. rwd = windrose2polar(wd)
    
16. 
    
17. #Degree to radians
    
18. rwd = radians(rwd)
    
19. 
    
20. #Calculate frequency of each wind direction bin
    
21. wdbins = linspace(0.0, pi * 2, 9)
    
22. dwdbins = degrees(wdbins)
    
23. dwdbins = windrose2polar(dwdbins)
    
24. rwdbins = radians(dwdbins)
    
25. wdN = len(wdbins) - 1
    
26. theta = ones(wdN + 1)
    
27. for j in range(wdN):
    
28.     theta[j] = rwdbins[j]
    
29. theta[wdN] = theta[0]
    
30. wd = wd + 360./wdN/2
    
31. wd[wd>360] = wd - 360
    
32. wdhist = histogram(radians(wd), wdbins)[0].astype('float')     
    
33. wdhist = wdhist / N
    
34. nwdhist = wdhist.aslist()
    
35. nwdhist.append(nwdhist[0])
    
36. nwdhist = array(nwdhist)
    
37. 
    
38. #Polar coordinate to Cartesian coordinate
    
39. rwdc, wsc = pol2cart(rwd, ws)
    
40. 
    
41. #Get convexhull (minimum outer polygon of the wind points)
    
42. poly = topo.convexhull(rwdc, wsc)
    
43. 
    
44. #Get grid data
    
45. dd = 0.5
    
46. x = linspace(rwdc.min() - dd, rwdc.max() + dd, 50)
    
47. y = linspace(wsc.min() - dd, wsc.max() + dd, 50)
    
48. data = griddata((rwdc, wsc), pm10, xi=(x, y), method='idw', pointnum=5, convexhull=False)[0]
    
49. 
    
50. #---------------------------------------
    
51. #Plot figure
    
52. pos = [0.13, 0.1, 0.775, 0.775]
    
53. 
    
54. #Cartesian axes
    
55. ax = axes(position=pos, aspect='equal')
    
56. xaxis(visible=False)
    
57. yaxis(location='right', visible=False)
    
58. yaxis(location='left', shift = 50)
    
59. ylabel('Wind speed (m/s)')
    
60. levs = arange(100, 2000, 100)
    
61. cg = contourf(x, y, data, levs, cmap='BlAqGrYeOrRe', visible=False)
    
62. cg = cg.clip([poly])
    
63. ax.add_graphic(cg)
    
64. #ss = scatter(rwdc, wsc, s=6, fill=False, edgecolor='b', edgesize=1)
    
65. #patch(poly)
    
66. colorbar(cg, shrink=0.6, xshift=30, label=r'$\mu g/m^3$', labelloc='bottom')
    
67. maxv = 10
    
68. xlim(-maxv, maxv)
    
69. ylim(-maxv, maxv)
    
70. ticks = ax.get_yticks()
    
71. ax.set_yticklabels([abs(yy) for yy in ticks])
    
72. 
    
73. #Polar axes
    
74. axp = axes(position=pos, polar=True)
    
75. plot(theta, nwdhist, color='k', linewidth=2)
    
76. axp.set_rmax(1)
    
77. axp.set_rlabel_position(25.)
    
78. axp.set_rtick_locations([0.2,0.4,0.6,0.8,1])
    
79. #axp.set_rticks(['2','4','6','8','10 m/s'])
    
80. axp.set_rticks(['20%','40%','60%','80%','100%'])
    
81. axp.set_xtick_font(size=14)
    
82. axp.set_xticks(['E','NE','N','NW','W','SW','S','SE'])
    
83. title(r'$Windrose \ with \ PM_{10} \ concentrations$', fontsize=18)

def windrose2polar(a):<br />     r = 360 - a + 90<br />     r[r>360] = r - 360<br />     return r<br /> <br /> #Read data (wind speed, weed direction, pm10)<br /> fn = os.path.join(migl.get_sample_folder(), 'ASCII', 'pm10.txt')<br /> df = DataFrame.read_table(fn, format='%3f')<br /> ws = df['WS'].values<br /> wd = df['WD'].values<br /> pm10 = df['PM10'].values<br /> N = len(ws)<br /> <br /> #Convert from windrose coordinate to polar coordinate<br /> rwd = windrose2polar(wd)<br /> <br /> #Degree to radians<br /> rwd = radians(rwd)<br /> <br /> #Calculate frequency of each wind direction bin<br /> wdbins = linspace(0.0, pi * 2, 9)<br /> dwdbins = degrees(wdbins)<br /> dwdbins = windrose2polar(dwdbins)<br /> rwdbins = radians(dwdbins)<br /> wdN = len(wdbins) - 1<br /> theta = ones(wdN + 1)<br /> for j in range(wdN):<br />     theta[j] = rwdbins[j]<br /> theta[wdN] = theta[0]<br /> wd = wd + 360./wdN/2<br /> wd[wd>360] = wd - 360<br /> wdhist = histogram(radians(wd), wdbins)[0].astype('float')     <br /> wdhist = wdhist / N<br /> nwdhist = wdhist.aslist()<br /> nwdhist.append(nwdhist[0])<br /> nwdhist = array(nwdhist)<br /> <br /> #Polar coordinate to Cartesian coordinate<br /> rwdc, wsc = pol2cart(rwd, ws)<br /> <br /> #Get convexhull (minimum outer polygon of the wind points)<br /> poly = topo.convexhull(rwdc, wsc)<br /> <br /> #Get grid data<br /> dd = 0.5<br /> x = linspace(rwdc.min() - dd, rwdc.max() + dd, 50)<br /> y = linspace(wsc.min() - dd, wsc.max() + dd, 50)<br /> data = griddata((rwdc, wsc), pm10, xi=(x, y), method='idw', pointnum=5, convexhull=False)[0]<br /> <br /> #---------------------------------------<br /> #Plot figure<br /> pos = [0.13, 0.1, 0.775, 0.775]<br /> <br /> #Cartesian axes<br /> ax = axes(position=pos, aspect='equal')<br /> xaxis(visible=False)<br /> yaxis(location='right', visible=False)<br /> yaxis(location='left', shift = 50)<br /> ylabel('Wind speed (m/s)')<br /> levs = arange(100, 2000, 100)<br /> cg = contourf(x, y, data, levs, cmap='BlAqGrYeOrRe', visible=False)<br /> cg = cg.clip([poly])<br /> ax.add_graphic(cg)<br /> #ss = scatter(rwdc, wsc, s=6, fill=False, edgecolor='b', edgesize=1)<br /> #patch(poly)<br /> colorbar(cg, shrink=0.6, xshift=30, label=r'$\mu g/m^3$', labelloc='bottom')<br /> maxv = 10<br /> xlim(-maxv, maxv)<br /> ylim(-maxv, maxv)<br /> ticks = ax.get_yticks()<br /> ax.set_yticklabels([abs(yy) for yy in ticks])<br /> <br /> #Polar axes<br /> axp = axes(position=pos, polar=True)<br /> plot(theta, nwdhist, color='k', linewidth=2)<br /> axp.set_rmax(1)<br /> axp.set_rlabel_position(25.)<br /> axp.set_rtick_locations([0.2,0.4,0.6,0.8,1])<br /> #axp.set_rticks(['2','4','6','8','10 m/s'])<br /> axp.set_rticks(['20%','40%','60%','80%','100%'])<br /> axp.set_xtick_font(size=14)<br /> axp.set_xticks(['E','NE','N','NW','W','SW','S','SE'])<br /> title(r'$Windrose \ with \ PM_{10} \ concentrations$', fontsize=18)

estgree

MeteoInfo 发表于 2018-1-10 18:46
下载最新版本的MeteoInfo

老师您好，按照这个脚本，一直显示的 run script...  但是没有出图，是什么原因呀

SYNJU

王老师，您好。请问我如果要把PM10、PM2.5、SO2、NO2、CO、O3空气质量六参的污染玫瑰做成一张图，有什么实现方法，尝试了一下subplo函数，发现实现不了。

mayfin

还有个问题想请教下，如果在某个风向风速下有很多PM10的数值，请问填色的是指该风向风速下的PM10浓度平均值还是最大值？

lyysms

请问老师，风速超过10, 内部图形超过了圆圈，这个应如何处理呢？@MeteoInfo

billy12323

好东东，顶一个

irregular

这个有点意思。。。

river

漂亮！好东西啊，就是零基础学起来太难了

MeteoInfo

river 发表于 2017-6-30 22:34
漂亮！好东西啊，就是零基础学起来太难了

这个在数据处理、绘图语言里已经是比较容易学的了。

戏子丶

哇 我之前看到有人在问这个问题怎么画。想不到王老师这么快就搞定了。

戏子丶

另外就是，王老师，如果我学了py后，弄这个是不是就很轻松啊？

rceclx

看上去很漂亮的！向王老师学习！学习！再学习！

MeteoInfo

戏子丶 发表于 2017-7-1 11:13
另外就是，王老师，如果我学了py后，弄这个是不是就很轻松啊？

是的

xubaisen00000

大家好，我想请教下，你们写的窗体 不断滚动鼠标滑轮放大mapview 不会提示“System.OverflowException”类型的未经处理的异常在 System.Drawing.dll 中发生 这个错误吗 我看了好几个例子都有这个错误

xubaisen00000

还有请教下王老师，怎么控制MeteInfo中mapview的放大倍数，谢谢！

MeteoInfo

xubaisen00000 发表于 2017-7-3 10:04
还有请教下王老师，怎么控制MeteInfo中mapview的放大倍数，谢谢！

请不要在帖子里问不相干的问题，可以重开一个帖子提问。