MeteoInfoLab脚本示例：计算涡度并插值

MeteoInfo

从4维（时间、高度、纬度、经度）气象格点数据中读取风场U/V分量4维数组，然后用vort函数计算涡度，读取气团轨迹节点的4维坐标，将涡度数据插值到轨迹节点上（用到interpn函数）。

1. # Set working directory
   
2. trajDir = 'D:/Temp/HYSPLIT'
   
3. meteoDir = r'U:\data\ARL\2015'
   
4. 
   
5. # Open trjactory data file
   
6. print 'Open trajectory data file ...'
   
7. trajfn = os.path.join(trajDir, 'traj_20150331')
   
8. print 'Trajectory file: ' + trajfn
   
9. trajf = addfile_hytraj(trajfn)
   
10. 
    
11. # Create trajectory layer
    
12. trajLayer = trajf.trajlayer()
    
13. 
    
14. # Open meteorological data file
    
15. print 'Open meteorological data file...'
    
16. meteofn = os.path.join(meteoDir, 'gdas1.mar15.w5')
    
17. print 'Meteorological file: ' + meteofn
    
18. meteof = addfile(meteofn)
    
19. 
    
20. # Get meteorological data along trajectory
    
21. print 'Get meteorological data along trajectory...'
    
22. outfn = os.path.join(trajDir, 'traj.txt')
    
23. outf = open(outfn, 'w')
    
24. outf.write('Lon,Lat,Time,Heigh,pres,VORT\n')
    
25. #uwnd = meteof['UWND'][:,:,:,:]
    
26. #vwnd = meteof['VWND'][:,:,:,:]
    
27. uwnd = meteof['UWND'][:,:,[10, 80],[50,140]]
    
28. vwnd = meteof['VWND'][:,:,[10, 80],[50,140]]
    
29. vort = hcurl(uwnd, vwnd)
    
30. idx = 0
    
31. for tline in trajLayer.shapes():   
    
32.     t = trajLayer.cellvalue('Date', idx)
    
33.     h = trajLayer.cellvalue('Hour', idx)
    
34.     t.replace(hour=h)
    
35.     for ps in tline.getPoints():  
    
36.         lon = ps.X
    
37.         lat = ps.Y
    
38.         z = ps.M
    
39.         pres = ps.Z
    
40.         v = vort.interpn([t,pres,lat,lon])
    
41.         print 'lon: %.3f; lat: %.3f; time: %s; height: %.1f; pres: %.1f; VORT: %f' % (lon,lat,t.strftime('%Y%m%d_%H:%M'),z,pres,v)
    
42.         line = '%.3f,%.3f,%s,%.1f,%.1f,%f' % (lon,lat,t.strftime('%Y%m%d_%H:%M'),z,pres,v)
    
43.         outf.write(line + '\n')
    
44.         t = t + datetime.timedelta(hours=-1)
    
45.     idx += 1
    
46. outf.close()
    
47. print 'Finish...'

# Set working directory<br /> trajDir = 'D:/Temp/HYSPLIT'<br /> meteoDir = r'U:\data\ARL\2015'<br /> <br /> # Open trjactory data file<br /> print 'Open trajectory data file ...'<br /> trajfn = os.path.join(trajDir, 'traj_20150331')<br /> print 'Trajectory file: ' + trajfn<br /> trajf = addfile_hytraj(trajfn)<br /> <br /> # Create trajectory layer<br /> trajLayer = trajf.trajlayer()<br /> <br /> # Open meteorological data file<br /> print 'Open meteorological data file...'<br /> meteofn = os.path.join(meteoDir, 'gdas1.mar15.w5')<br /> print 'Meteorological file: ' + meteofn<br /> meteof = addfile(meteofn)<br /> <br /> # Get meteorological data along trajectory<br /> print 'Get meteorological data along trajectory...'<br /> outfn = os.path.join(trajDir, 'traj.txt')<br /> outf = open(outfn, 'w')<br /> outf.write('Lon,Lat,Time,Heigh,pres,VORT\n')<br /> #uwnd = meteof['UWND'][:,:,:,:]<br /> #vwnd = meteof['VWND'][:,:,:,:]<br /> uwnd = meteof['UWND'][:,:,[10, 80],[50,140]]<br /> vwnd = meteof['VWND'][:,:,[10, 80],[50,140]]<br /> vort = hcurl(uwnd, vwnd)<br /> idx = 0<br /> for tline in trajLayer.shapes():    <br />     t = trajLayer.cellvalue('Date', idx)<br />     h = trajLayer.cellvalue('Hour', idx)<br />     t.replace(hour=h)<br />     for ps in tline.getPoints():  <br />         lon = ps.X<br />         lat = ps.Y <br />         z = ps.M<br />         pres = ps.Z<br />         v = vort.interpn([t,pres,lat,lon])<br />         print 'lon: %.3f; lat: %.3f; time: %s; height: %.1f; pres: %.1f; VORT: %f' % (lon,lat,t.strftime('%Y%m%d_%H:%M'),z,pres,v)<br />         line = '%.3f,%.3f,%s,%.1f,%.1f,%f' % (lon,lat,t.strftime('%Y%m%d_%H:%M'),z,pres,v)<br />         outf.write(line + '\n')<br />         t = t + datetime.timedelta(hours=-1)<br />     idx += 1<br /> outf.close()<br /> print 'Finish...'

良辰

万分感谢王老师！真想做您的学生！

MeteoInfo

良辰 发表于 2016-8-6 12:05
万分感谢王老师！真想做您的学生！

interpn是新加的函数（相当于多维线性插值），挺费劲的。

欢迎报考我的研究生，或者做访问学者。

良辰

MeteoInfo 发表于 2016-8-6 12:22
interpn是新加的函数（相当于多维线性插值），挺费劲的。

欢迎报考我的研究生，或者做访问学者。

好的，一定积极争取！有机会欢迎王老师来长春市气象局来做报告。