MeteoInfoLab脚本示例：计算水平螺旋度

MeteoInfo

参考此贴（http://bbs.06climate.com/forum.php?mod=viewthread&tid=37089）尝试了用MeteoInfoLab编写计算水平螺旋度的脚本，由于我不是学气象出身，只是照猫画虎，希望有高手验证一下是否正确。

对平均风速的计算根据river的建议试了两种方法，1、先计算U/V分量的平均值，再利用U/V分量平均值计算平均风速；2、先计算全风速（所有U/V数组元素都计算风速），再平均。结果对于Pattern没有影响，但对具体的值有影响。

脚本程序：

1. print 'Open data files...'
   
2. f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc')
   
3. f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc')
   
4. 
   
5. print 'Calculate average wind field from 850 to 600 hpa...'
   
6. tidx = 173    # Jun 23, 2011
   
7. t = f_uwnd.gettime(tidx)
   
8. level = '850:600'    # 850 - 600 hPa
   
9. lat = '15:55'
   
10. lon = '70:135'
    
11. uvar = f_uwnd['uwnd']
    
12. vvar = f_vwnd['vwnd']
    
13. uwnd = uvar[tidx,level,lat,lon]
    
14. vwnd = vvar[tidx,level,lat,lon]
    
15. uc = uwnd.ave()
    
16. vc = vwnd.ave()sp = magnitude(uwnd, vwnd)
    
17. speed = sp.ave()
    
18. #speed = magnitude(uc, vc)
    
19. direc = atan2(vc, uc)
    
20. cdirec = direc-40./180.*3.14159
    
21. if cdirec < -3.14159:
    
22.     cdirec=cdirec+3.14159*2
    
23. cuc = speed*cos(cdirec)
    
24. cvc = speed*sin(cdirec)
    
25. 
    
26. print 'Calculate horizontal helicity...'
    
27. Hrs = 0.
    
28. zHrs = 0.
    
29. hh = 6
    
30. while hh <= 10:
    
31.     zHrs = (uvar[tidx,hh+1,lat,lon]-cuc)*(vvar[tidx,hh,lat,lon]-cvc)- \
    
32.         (uvar[tidx,hh,lat,lon]-cuc)*(vvar[tidx,hh+1,lat,lon]-cvc)
    
33.     Hrs = Hrs+zHrs
    
34.     hh += 1
    
35. Hrs = Hrs[::-1,:]
    
36. 
    
37. print 'Plot...'
    
38. axesm()
    
39. mlayer = shaperead('D:/Temp/map/country1.shp')
    
40. geoshow(mlayer, edgecolor='black')
    
41. layer = contourfm(Hrs, 20)
    
42. title('Horizontal helicity (' + t.strftime('%Y-%m-%d') + ')')
    
43. colorbar(layer)

print 'Open data files...'<br /> f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc')<br /> f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc')<br /> <br /> print 'Calculate average wind field from 850 to 600 hpa...'<br /> tidx = 173    # Jun 23, 2011<br /> t = f_uwnd.gettime(tidx)<br /> level = '850:600'    # 850 - 600 hPa<br /> lat = '15:55'<br /> lon = '70:135'<br /> uvar = f_uwnd['uwnd']<br /> vvar = f_vwnd['vwnd']<br /> uwnd = uvar[tidx,level,lat,lon]<br /> vwnd = vvar[tidx,level,lat,lon]<br /> uc = uwnd.ave()<br /> vc = vwnd.ave()sp = magnitude(uwnd, vwnd)<br /> speed = sp.ave()<br /> #speed = magnitude(uc, vc)<br /> direc = atan2(vc, uc)<br /> cdirec = direc-40./180.*3.14159<br /> if cdirec < -3.14159:<br />     cdirec=cdirec+3.14159*2<br /> cuc = speed*cos(cdirec)<br /> cvc = speed*sin(cdirec)<br /> <br /> print 'Calculate horizontal helicity...'<br /> Hrs = 0.<br /> zHrs = 0.<br /> hh = 6<br /> while hh <= 10:<br />     zHrs = (uvar[tidx,hh+1,lat,lon]-cuc)*(vvar[tidx,hh,lat,lon]-cvc)- \<br />         (uvar[tidx,hh,lat,lon]-cuc)*(vvar[tidx,hh+1,lat,lon]-cvc)<br />     Hrs = Hrs+zHrs<br />     hh += 1<br /> Hrs = Hrs[::-1,:]<br /> <br /> print 'Plot...'<br /> axesm()<br /> mlayer = shaperead('D:/Temp/map/country1.shp')<br /> geoshow(mlayer, edgecolor='black')<br /> layer = contourfm(Hrs, 20)<br /> title('Horizontal helicity (' + t.strftime('%Y-%m-%d') + ')')<br /> colorbar(layer)

脚本运行结果图：

第一种计算平均风速的结果：


第二种计算平均风速的结果：

兰溪之水

好像每种语言都有这个函数，atan2(y,x)是坐标原点为起点，指向(x,y)的射线在坐标平面上与x轴正方向之间的角的角度, 范围介于 -pi 到 pi 之间（不包括 -pi），这个函数在求角度方面很有用，可以省很多判断。

river

我在想求平均风那一块，平均风矢量只能是先求各分量的平均，然后合成；那平均风速就可以有两种方法，第一种就像王老师这样，第二种就是先求出全风速，然后再平均。哪一个更合理呢？

MeteoInfo

兰溪之水 发表于 2015-10-8 07:02
好像每种语言都有这个函数，atan2(y,x)是坐标原点为起点，指向(x,y)的射线在坐标平面上与x轴正方向之间的角 ...

感谢你的提示，的确是这样。

MeteoInfo

river 发表于 2015-10-8 07:36
我在想求平均分那一块，平均风矢量只能是先求各分量的平均，然后合成；那平均风速就可以有两种方法，第一种 ...

感觉你说的第二种更合理一些，我再修改试试。

river

MeteoInfo 发表于 2015-10-8 09:47
感觉你说的第二种更合理一些，我再修改试试。

随便用了几个数字，看了看两者间的差异，发现还是比较小的。不知道数据多了的话，这个差异会不会扩大···

catharine0607

还有一个问题，风暴的传播速度通常取为1.5~7.0km气层间的平均风速的75%且向右偏转40°，因此，speed还应乘以0.75.

geden

请教一下老师

Hrs = Hrs[::-1,:]

请问这一行的作用是什么呢？为什么要翻过来