MeteoInfoLab脚本示例:计算水平螺旋度
尝试了用MeteoInfoLab编写计算水平螺旋度的脚本,结果未经验证。
脚本程序:
print 'Open data files...' f_uwnd = addfile('D:/Temp/nc/uwnd.2011.nc') f_vwnd = addfile('D:/Temp/nc/vwnd.2011.nc') print 'Calculate average wind field from 850 to 600 hpa...' tidx = 173 # Jun 23, 2011 t = f_uwnd.gettime(tidx) level = [850,600] # 850 - 600 hPa lat = [15,55] lon = [70,135] uvar = f_uwnd['uwnd'] vvar = f_vwnd['vwnd'] uwnd = uvar[tidx,level,lat,lon] vwnd = vvar[tidx,level,lat,lon] uc = uwnd.ave() vc = vwnd.ave() sp = magnitude(uwnd, vwnd) speed = sp.ave() #speed = magnitude(uc, vc) direc = atan2(vc, uc) cdirec = direc-40./180.*3.14159 if cdirec < -3.14159: cdirec=cdirec+3.14159*2 cuc = speed*cos(cdirec) cvc = speed*sin(cdirec) print 'Calculate horizontal helicity...' Hrs = 0. zHrs = 0. hh = 6 while hh <= 10: zHrs = (uvar[tidx,hh+1,lat,lon]-cuc)*(vvar[tidx,hh,lat,lon]-cvc)- \ (uvar[tidx,hh,lat,lon]-cuc)*(vvar[tidx,hh+1,lat,lon]-cvc) Hrs = Hrs+zHrs hh += 1 Hrs = Hrs[::-1,:] print 'Plot...' axesm() mlayer = shaperead('D:/Temp/map/country1.shp') geoshow(mlayer, edgecolor='black') layer = contourfm(Hrs, 20) title('Horizontal helicity (' + t.strftime('%Y-%m-%d') + ')') colorbar(layer)