可以计算一整年Trajectory的脚本

MeteoInfo

利用HYSPLIT GUI或者TrajStat可以一次计算一个月的气团轨迹，但对于要大批量计算长时间序列的轨迹来说还是不够方便。这里给出一个能够一次计算一年轨迹的脚本，会极大地减少工作量。需要最新的MeteoInfo Java版软件和TrajStat Java版plugin。

1. #-----------------------------------------------------
   
2. # Author: Yaqiang Wang
   
3. # Date: 2014-7-11
   
4. # Purpose: Run TrajStat - calculate trajectories
   
5. # Note: Sample
   
6. #-----------------------------------------------------
   
7. import sys
   
8. sys.path.append('D:/MyProgram/java/MeteoInfoDev/MeteoInfo/plugins/TrajStat/TrajStat.jar')
   
9. import os
   
10. import time
    
11. from datetime import datetime
    
12. from datetime import timedelta
    
13. from calendar import monthrange
    
14. from trajstat.trajectory import TrajUtil
    
15. from trajstat.trajectory import TrajConfig
    
16. 
    
17. # Set working directory
    
18. metDir = 'P:/MetData/MetData'
    
19. outDir = 'O:/trajectory/bishkek/trajdata'
    
20. 
    
21. # Get month abstract string
    
22. def getmonthstr(m):  
    
23.   mmm = 'jan'
    
24.   if m == 1:
    
25.     mmm = 'jan'
    
26.   elif m == 2:
    
27.     mmm = 'feb'
    
28.   elif m == 3:
    
29.     mmm = 'mar'
    
30.   elif m == 4:
    
31.     mmm = 'apr'
    
32.   elif m == 5:
    
33.     mmm = 'may'
    
34.   elif m == 6:
    
35.     mmm = 'jun'
    
36.   elif m == 7:
    
37.     mmm = 'jul'
    
38.   elif m == 8:
    
39.     mmm = 'aug'
    
40.   elif m == 9:
    
41.     mmm = 'sep'
    
42.   elif m == 10:
    
43.     mmm = 'oct'
    
44.   elif m == 11:
    
45.     mmm = 'nov'
    
46.   elif m == 12:
    
47.     mmm = 'dec'
    
48. 
    
49.   return mmm
    
50. 
    
51. # Get GDAS1 meteorological data files by time
    
52. def getmeteofiles(y, m):
    
53.   ystr = str(y)[2:]
    
54.   mmm = getmonthstr(m)
    
55.   fns = []
    
56.   # The meteo files of this month
    
57.   cdir = os.path.join(metDir, str(y))
    
58.   for i in range(1,6):
    
59.     fn = os.path.join(cdir, 'gdas1.' + mmm + ystr + '.w' + str(i))
    
60.     if os.path.exists(fn):
    
61.       fns.append(fn)
    
62. 
    
63.   # The last two meteo files of last month  
    
64.   m = m - 1
    
65.   if m == 0:
    
66.     m = 12
    
67.     y = y - 1
    
68.     ystr = str(y)[2:]
    
69.     cdir = os.path.join(metDir, str(y))
    
70.   mmm = getmonthstr(m)
    
71.   fn = os.path.join(cdir, 'gdas1.' + mmm + ystr + '.w4')
    
72.   fns.append(fn)
    
73.   fn = os.path.join(cdir, 'gdas1.' + mmm + ystr + '.w5')
    
74.   if not os.path.exists(fn):
    
75.     fn = os.path.join(cdir, 'gdas1.' + mmm + ystr + '.w3')
    
76.   fns.append(fn)
    
77. 
    
78.   return fns
    
79. 
    
80. # Set parameters
    
81. print 'Set paramters...'
    
82. year = 2011
    
83. tc = TrajConfig()
    
84. tc.setStartDay(1)
    
85. tc.setStartHours('00 06 12 18')
    
86. tc.setLocation('42.68 74.69 1000.00')
    
87. tc.setRunHours(-96)
    
88. tc.setTopOfModel(10000.0)
    
89. tc.setVertical(0)
    
90. tc.setOutPath(outDir)
    
91. 
    
92. # Loop - one year calculation
    
93. print 'Trajectory calculation...'
    
94. for month in range(1, 13):
    
95.   print 'Month: ' + str(month)  
    
96.   eday = monthrange(year, month)[1]
    
97.   tc.setEndDay(eday)
    
98.   tc.initStartTime(year, month)
    
99.   fns = getmeteofiles(year, month)
    
100.   tc.setMeteoFiles(fns)
     
101.   #TrajUtil.trajCal(tc)
     
102.   print 'Convert trajectory endpoint files to tgs files...'
     
103.   TrajUtil.trajToTGS(tc)
     
104.   print 'Join tgs files...'
     
105.   tgsfn = TrajUtil.joinTGSFiles(tc)
     
106.   print 'Convert tgs file to shape file...'
     
107.   shpfn = tgsfn.replace('.tgs','.shp')
     
108.   TrajUtil.convertToShapeFile(tgsfn, shpfn)
     
109.   print 'Remove trajectory files...'
     
110.   TrajUtil.removeTrajFiles(tc)
     
111. 
     
112. print 'Finish !'

#-----------------------------------------------------<br /> # Author: Yaqiang Wang<br /> # Date: 2014-7-11<br /> # Purpose: Run TrajStat - calculate trajectories<br /> # Note: Sample<br /> #-----------------------------------------------------<br /> import sys<br /> sys.path.append('D:/MyProgram/java/MeteoInfoDev/MeteoInfo/plugins/TrajStat/TrajStat.jar')<br /> import os<br /> import time<br /> from datetime import datetime<br /> from datetime import timedelta<br /> from calendar import monthrange<br /> from trajstat.trajectory import TrajUtil<br /> from trajstat.trajectory import TrajConfig<br /> <br /> # Set working directory<br /> metDir = 'P:/MetData/MetData'<br /> outDir = 'O:/trajectory/bishkek/trajdata'<br /> <br /> # Get month abstract string<br /> def getmonthstr(m):  <br />   mmm = 'jan'<br />   if m == 1:<br />     mmm = 'jan'<br />   elif m == 2:<br />     mmm = 'feb'<br />   elif m == 3:<br />     mmm = 'mar'<br />   elif m == 4:<br />     mmm = 'apr'<br />   elif m == 5:<br />     mmm = 'may'<br />   elif m == 6:<br />     mmm = 'jun'<br />   elif m == 7:<br />     mmm = 'jul'<br />   elif m == 8:<br />     mmm = 'aug'<br />   elif m == 9:<br />     mmm = 'sep'<br />   elif m == 10:<br />     mmm = 'oct'<br />   elif m == 11:<br />     mmm = 'nov'<br />   elif m == 12:<br />     mmm = 'dec'<br /> <br />   return mmm<br /> <br /> # Get GDAS1 meteorological data files by time<br /> def getmeteofiles(y, m):<br />   ystr = str(y)[2:]<br />   mmm = getmonthstr(m)<br />   fns = []<br />   # The meteo files of this month<br />   cdir = os.path.join(metDir, str(y))<br />   for i in range(1,6):<br />     fn = os.path.join(cdir, 'gdas1.' + mmm + ystr + '.w' + str(i))<br />     if os.path.exists(fn):<br />       fns.append(fn)<br /> <br />   # The last two meteo files of last month  <br />   m = m - 1<br />   if m == 0:<br />     m = 12<br />     y = y - 1<br />     ystr = str(y)[2:]<br />     cdir = os.path.join(metDir, str(y))<br />   mmm = getmonthstr(m)<br />   fn = os.path.join(cdir, 'gdas1.' + mmm + ystr + '.w4')<br />   fns.append(fn)<br />   fn = os.path.join(cdir, 'gdas1.' + mmm + ystr + '.w5')<br />   if not os.path.exists(fn):<br />     fn = os.path.join(cdir, 'gdas1.' + mmm + ystr + '.w3')<br />   fns.append(fn)<br /> <br />   return fns<br /> <br /> # Set parameters<br /> print 'Set paramters...'<br /> year = 2011<br /> tc = TrajConfig()<br /> tc.setStartDay(1)<br /> tc.setStartHours('00 06 12 18')<br /> tc.setLocation('42.68 74.69 1000.00')<br /> tc.setRunHours(-96)<br /> tc.setTopOfModel(10000.0)<br /> tc.setVertical(0)<br /> tc.setOutPath(outDir)<br /> <br /> # Loop - one year calculation<br /> print 'Trajectory calculation...'<br /> for month in range(1, 13):<br />   print 'Month: ' + str(month)  <br />   eday = monthrange(year, month)[1]<br />   tc.setEndDay(eday)<br />   tc.initStartTime(year, month)<br />   fns = getmeteofiles(year, month)<br />   tc.setMeteoFiles(fns)<br />   #TrajUtil.trajCal(tc)<br />   print 'Convert trajectory endpoint files to tgs files...'<br />   TrajUtil.trajToTGS(tc)<br />   print 'Join tgs files...'<br />   tgsfn = TrajUtil.joinTGSFiles(tc)<br />   print 'Convert tgs file to shape file...'<br />   shpfn = tgsfn.replace('.tgs','.shp')<br />   TrajUtil.convertToShapeFile(tgsfn, shpfn)<br />   print 'Remove trajectory files...'<br />   TrajUtil.removeTrajFiles(tc)<br /> <br /> print 'Finish !'

rceclx

沙发！王老师新作！

cbsdly792

多谢分享。。。

piwa

1、要算一阵年的反向轨迹，然后进行聚类，必须在metroinfo里面用上述程序做吗？目前在trajstat中不能直接算吗？
2、进行反向轨迹计算，并在此基础上进行PSCF分析的数据是在Gridded Meteorological Data Archives中下载WRF ARW (Weather Research & Forecasting Model for CAPTEX Tracer Experiment)中的数据么？
菜鸟一只，望不要嫌弃！

piwa

数据下载位置是这个吗？

洛小舞

王老师，这个计算一年的后向轨迹的脚本，怎么用到TrajStat里面啊？只运行脚本的话，不能得到什么结果啊/(ㄒoㄒ)/~~

KIMO23

{:eb502:} 謝謝您的分享

dtlyr

修改了文件的路径，并且101.  #TrajUtil.trajCal(tc)的#去掉，结果出现
Percent complete:  99.0
Percent complete: 100.0
Complete Hysplit
HYSPLIT49 (Oct 2008) - Initialization
  NOTICE: using namelist file
不停的循环，是什么原因，如何停止脚本，并正确的计算一年轨迹呢？

SETUP.CFG              

MeteoInfo

dtlyr 发表于 2016-7-17 12:58
修改了文件的路径，并且101.  #TrajUtil.trajCal(tc)的#去掉，结果出现
Percent complete:  99.0
Perce ...

就是要循环计算呀，不然怎么能一次计算一年的轨迹呢。

dtlyr

MeteoInfo 发表于 2016-7-17 19:53
就是要循环计算呀，不然怎么能一次计算一年的轨迹呢。

如果中途要终止程序应该如何操作？
计算一年大概需要多久？运行了20多分钟都没有结果呢？
感觉code对于不熟悉jython 的新手操作起来还是有难度，如果有详细的如何运行出结果的例子就好了？