MeteoInfoLab脚本示例：计算气团轨迹

MeteoInfo

这里假设计算机上已经安装了HYSPLIT（或者TrajStat），具备了气团轨迹的计算环境。其实就是写CONTROL文件，再调用hyts_std.exe计算气团轨迹。MeteoInfoLab基于Jython，下面的脚本中并未用到MeteoInfLab中封装的函数，因此是纯粹的Python代码。

计算一天的轨迹：

1. #-----------------------------------------------------# Author: Yaqiang Wang
   
2. # Date: 2015-9-30
   
3. # Purpose: Run HYSPLIT trajectory model
   
4. # Note: Sample
   
5. #-----------------------------------------------------
   
6. 
   
7. # Set working directory
   
8. metDir = 'D:/Temp/arl'
   
9. outDir = 'D:/Temp/HYSPLIT'
   
10. workingDir = 'C:/hysplit4/working'
    
11. os.chdir(workingDir)
    
12. print 'Current directory: ' + os.getcwd()
    
13. 
    
14. # Set parameters
    
15. lon = '115.2'
    
16. lat = '40.1'
    
17. shour = '06'
    
18. heights = ['100.0','500.0','1000.0']
    
19. hnum = len(heights)
    
20. hours = '-48'
    
21. vertical = '0'
    
22. top = '10000.0'
    
23. 
    
24. # Set meteorological data files
    
25. fns = []
    
26. fn = 'gdas1.jul09.w5'
    
27. fns.append(fn)
    
28. 
    
29. # Set start/end time
    
30. stime = datetime.datetime(2009,7,31)
    
31. 
    
32. # Write CONTROL file
    
33. ctFile = './CONTROL'
    
34. print stime.strftime('%Y-%m-%d ') + shour + ':00'
    
35. ctf = open(ctFile, 'w')
    
36. ctf.write(stime.strftime('%y %m %d ') + shour + "\n")
    
37. ctf.write(str(hnum) + '\n')
    
38. for j in range(0,hnum):
    
39.     ctf.write(lat + ' ' + lon + ' ' + heights[j] + '\n')
    
40. ctf.write(hours + '\n')
    
41. ctf.write(vertical + '\n')
    
42. ctf.write(top + '\n')
    
43. fnnum = len(fns)
    
44. ctf.write(str(fnnum) + '\n')
    
45. for j in range(0,fnnum):
    
46.     ctf.write(metDir + '/' + '\n')
    
47.     ctf.write(fns[j] + '\n')
    
48. ctf.write(outDir + '/' + '\n')
    
49. outfn = stime.strftime('traj_%Y%m%d')
    
50. ctf.write(outfn)
    
51. ctf.close()
    
52. 
    
53. # Calculate trajectories
    
54. os.system('c:/hysplit4/exec/hyts_std.exe')
    
55. 
    
56. print 'Finish...'

#-----------------------------------------------------# Author: Yaqiang Wang<br /> # Date: 2015-9-30<br /> # Purpose: Run HYSPLIT trajectory model<br /> # Note: Sample<br /> #-----------------------------------------------------<br /> <br /> # Set working directory<br /> metDir = 'D:/Temp/arl'<br /> outDir = 'D:/Temp/HYSPLIT'<br /> workingDir = 'C:/hysplit4/working'<br /> os.chdir(workingDir)<br /> print 'Current directory: ' + os.getcwd()<br /> <br /> # Set parameters<br /> lon = '115.2'<br /> lat = '40.1'<br /> shour = '06'<br /> heights = ['100.0','500.0','1000.0']<br /> hnum = len(heights)<br /> hours = '-48'<br /> vertical = '0'<br /> top = '10000.0'<br /> <br /> # Set meteorological data files<br /> fns = []<br /> fn = 'gdas1.jul09.w5'<br /> fns.append(fn)<br /> <br /> # Set start/end time<br /> stime = datetime.datetime(2009,7,31)<br /> <br /> # Write CONTROL file<br /> ctFile = './CONTROL'<br /> print stime.strftime('%Y-%m-%d ') + shour + ':00'<br /> ctf = open(ctFile, 'w')<br /> ctf.write(stime.strftime('%y %m %d ') + shour + "\n")<br /> ctf.write(str(hnum) + '\n')<br /> for j in range(0,hnum):<br />     ctf.write(lat + ' ' + lon + ' ' + heights[j] + '\n')<br /> ctf.write(hours + '\n')<br /> ctf.write(vertical + '\n')<br /> ctf.write(top + '\n')<br /> fnnum = len(fns)<br /> ctf.write(str(fnnum) + '\n')<br /> for j in range(0,fnnum):<br />     ctf.write(metDir + '/' + '\n')<br />     ctf.write(fns[j] + '\n')<br /> ctf.write(outDir + '/' + '\n')<br /> outfn = stime.strftime('traj_%Y%m%d')<br /> ctf.write(outfn)<br /> ctf.close()<br /> <br /> # Calculate trajectories<br /> os.system('c:/hysplit4/exec/hyts_std.exe')<br /> <br /> print 'Finish...'

计算多日的气团轨迹：

1. import calendar
   
2. 
   
3. # Set working directory
   
4. metDir = 'U:/data/ARL'
   
5. outDir = 'D:/Temp/HYSPLIT'
   
6. workingDir = 'C:/hysplit4/working'
   
7. os.chdir(workingDir)
   
8. print 'Current directory: ' + os.getcwd()
   
9. 
   
10. # Set parameters
    
11. lon = '115.2'
    
12. lat = '40.1'
    
13. shour = '06'
    
14. heights = ['100.0','500.0','1000.0']
    
15. hnum = len(heights)
    
16. hours = '-48'
    
17. vertical = '0'
    
18. top = '10000.0'
    
19. 
    
20. # Get GDAS1 meteorological data files by time
    
21. def getmeteofiles(t):
    
22.     ystr = t.strftime('%y')
    
23.     mdir = metDir + '/%s' % t.strftime('%Y')
    
24.     mmm = miutil.dateformat(t, 'MMM', 'us_en').lower()
    
25.     mdirs = []
    
26.     fns = []
    
27.     # The meteo files of this month
    
28.     for i in range(1,6):
    
29.         fn = 'gdas1.' + mmm + ystr + '.w' + str(i)
    
30.         if os.path.exists(os.path.join(mdir, fn)):
    
31.             mdirs.append(mdir)
    
32.             fns.append(fn)
    
33. 
    
34.     # The last two meteo files of last month
    
35.     days = calendar.monthrange(t.year, t.month)[1]
    
36.     t = t - datetime.timedelta(days=days)
    
37.     ystr = t.strftime('%y')
    
38.     mdir = metDir + '/%s' % t.strftime('%Y')
    
39.     mmm = miutil.dateformat(t, 'MMM', 'us_en').lower()
    
40.     fn = 'gdas1.' + mmm + ystr + '.w4'
    
41.     mdirs.append(mdir)
    
42.     fns.append(fn)
    
43.     fn = 'gdas1.' + mmm + ystr + '.w5'
    
44.     if os.path.exists(os.path.join(mdir, fn)):
    
45.         mdirs.append(mdir)
    
46.         fns.append(fn)
    
47.     else:
    
48.         mdirs.append(mdir)
    
49.         fns.append('gdas1.' + mmm + ystr + '.w3')
    
50. 
    
51.     return mdirs, fns
    
52. 
    
53. # Set start/end time
    
54. stime = datetime.datetime(2014,9,1)
    
55. etime = datetime.datetime(2014,10,1)
    
56. 
    
57. # Loop
    
58. ctFile = './CONTROL'
    
59. while stime < etime:
    
60.     print stime.strftime('%Y-%m-%d ') + shour + ':00'
    
61.     ctf = open(ctFile, 'w')
    
62.     ctf.write(stime.strftime('%y %m %d ') + shour + "\n")
    
63.     ctf.write(str(hnum) + '\n')
    
64.     for j in range(0,hnum):
    
65.         ctf.write(lat + ' ' + lon + ' ' + heights[j] + '\n')
    
66.     ctf.write(hours + '\n')
    
67.     ctf.write(vertical + '\n')
    
68.     ctf.write(top + '\n')
    
69.     mdirs, fns = getmeteofiles(stime)
    
70.     fnnum = len(fns)
    
71.     ctf.write(str(fnnum) + '\n')
    
72.     for mdir,fn in zip(mdirs,fns):
    
73.         ctf.write(mdir + '/' + '\n')
    
74.         ctf.write(fn + '\n')
    
75.     ctf.write(outDir + '/' + '\n')
    
76.     outfn = stime.strftime('traj_%Y%m%d')
    
77.     ctf.write(outfn)
    
78.     ctf.close()
    
79.     os.system('c:/hysplit4/exec/hyts_std.exe')
    
80. 
    
81.     stime = stime + datetime.timedelta(days=1)
    
82. 
    
83. print 'Finish...'

import calendar<br /> <br /> # Set working directory<br /> metDir = 'U:/data/ARL'<br /> outDir = 'D:/Temp/HYSPLIT'<br /> workingDir = 'C:/hysplit4/working'<br /> os.chdir(workingDir)<br /> print 'Current directory: ' + os.getcwd()<br /> <br /> # Set parameters<br /> lon = '115.2'<br /> lat = '40.1'<br /> shour = '06'<br /> heights = ['100.0','500.0','1000.0']<br /> hnum = len(heights)<br /> hours = '-48'<br /> vertical = '0'<br /> top = '10000.0'<br /> <br /> # Get GDAS1 meteorological data files by time<br /> def getmeteofiles(t):<br />     ystr = t.strftime('%y')<br />     mdir = metDir + '/%s' % t.strftime('%Y')<br />     mmm = miutil.dateformat(t, 'MMM', 'us_en').lower()<br />     mdirs = []<br />     fns = []<br />     # The meteo files of this month<br />     for i in range(1,6):<br />         fn = 'gdas1.' + mmm + ystr + '.w' + str(i)<br />         if os.path.exists(os.path.join(mdir, fn)):<br />             mdirs.append(mdir)<br />             fns.append(fn)<br /> <br />     # The last two meteo files of last month<br />     days = calendar.monthrange(t.year, t.month)[1]<br />     t = t - datetime.timedelta(days=days)<br />     ystr = t.strftime('%y')<br />     mdir = metDir + '/%s' % t.strftime('%Y')<br />     mmm = miutil.dateformat(t, 'MMM', 'us_en').lower()<br />     fn = 'gdas1.' + mmm + ystr + '.w4'<br />     mdirs.append(mdir)<br />     fns.append(fn)<br />     fn = 'gdas1.' + mmm + ystr + '.w5'<br />     if os.path.exists(os.path.join(mdir, fn)):<br />         mdirs.append(mdir)<br />         fns.append(fn)<br />     else:<br />         mdirs.append(mdir)<br />         fns.append('gdas1.' + mmm + ystr + '.w3')<br /> <br />     return mdirs, fns<br /> <br /> # Set start/end time<br /> stime = datetime.datetime(2014,9,1)<br /> etime = datetime.datetime(2014,10,1)<br /> <br /> # Loop<br /> ctFile = './CONTROL'<br /> while stime < etime:<br />     print stime.strftime('%Y-%m-%d ') + shour + ':00'<br />     ctf = open(ctFile, 'w')<br />     ctf.write(stime.strftime('%y %m %d ') + shour + "\n")<br />     ctf.write(str(hnum) + '\n')<br />     for j in range(0,hnum):<br />         ctf.write(lat + ' ' + lon + ' ' + heights[j] + '\n')<br />     ctf.write(hours + '\n')<br />     ctf.write(vertical + '\n')<br />     ctf.write(top + '\n')<br />     mdirs, fns = getmeteofiles(stime)<br />     fnnum = len(fns)<br />     ctf.write(str(fnnum) + '\n')<br />     for mdir,fn in zip(mdirs,fns):<br />         ctf.write(mdir + '/' + '\n')<br />         ctf.write(fn + '\n')<br />     ctf.write(outDir + '/' + '\n')<br />     outfn = stime.strftime('traj_%Y%m%d')<br />     ctf.write(outfn)<br />     ctf.close()<br />     os.system('c:/hysplit4/exec/hyts_std.exe')<br /> <br />     stime = stime + datetime.timedelta(days=1)<br /> <br /> print 'Finish...'

sharrow

感谢王老师分享

hbm

王老师，按照代码运行  ctf.write(lat + ' ' + lon + ' ' + heights<i> + '\n')
出现了错误提示，“TypeError: cannot concatenate 'str' and 'list' objects”，谢谢王老师了

MeteoInfo

hbm 发表于 2015-12-6 14:18
王老师，按照代码运行  ctf.write(lat + ' ' + lon + ' ' + heights + '\n')
出现了错误提示，“TypeError ...

论坛对代码的处理有点小问题，heights<i>是错误的，需要把尖括号<>改为方括号[]。也可以参考这里：http://www.meteothinker.com/exam ... ctory/traj_cal.html

hbm

谢谢，王老师，我再试试

二子搞研究

王老师，我把您一日的那个程序稍加改了一下，运行之后没有任何错误提示，但也没有生成当日轨迹图，求解

MeteoInfo

二子搞研究 发表于 2016-1-21 22:08
王老师，我把您一日的那个程序稍加改了一下，运行之后没有任何错误提示，但也没有生成当日轨迹图，求解{:eb ...

这个脚本只是计算轨迹，没有绘图的语句。

二子搞研究

MeteoInfo 发表于 2016-1-21 22:32
这个脚本只是计算轨迹，没有绘图的语句。

那老师如果想批量绘制图片的话是不是需要自己写程序了？

MeteoInfo

二子搞研究 发表于 2016-1-21 22:37
那老师如果想批量绘制图片的话是不是需要自己写程序了？

在这里找：MeteoInfoLab脚本汇总贴
http://bbs.06climate.com/forum.p ... 151&fromuid=106
(出处: 气象家园)

dtlyr

MeteoInfo 发表于 2016-1-21 22:32
这个脚本只是计算轨迹，没有绘图的语句。

脚本计算的轨迹不是像软件版的生成tgs文件吗？请问老师有没有关于meteinfolab脚本的教程？