MeteoInfoLab脚本示例：获取气团轨迹每个节点的气象数据

良辰

MeteoInfo 发表于 2016-8-5 16:45
参考此贴：MeteoInfoLab脚本示例：计算涡度并插值
http://bbs.06climate.com/forum.php?mod=viewthread& ...

1. # Set working directory
   
2. trajDir = 'D:/Temp/HYSPLIT'
   
3. meteoDir = r'D:\Temp\arl'
   
4. # Open trjactory data file
   
5. print 'Open trajectory data file ...'
   
6. trajfn = os.path.join(trajDir, 'traj_20150331')
   
7. print 'Trajectory file: ' + trajfn
   
8. trajf = addfile_hytraj(trajfn)
   
9. # Create trajectory layer
   
10. trajLayer = trajf.trajlayer()
    
11. # Open meteorological data file
    
12. print 'Open meteorological data file...'
    
13. meteofn = os.path.join(meteoDir, 'gdas1.mar15.w5')
    
14. print 'Meteorological file: ' + meteofn
    
15. meteof = addfile(meteofn)
    
16. # Get meteorological data along trajectory
    
17. print 'Get meteorological data along trajectory...'
    
18. outfn = os.path.join(trajDir, 'traj.txt')
    
19. outf = open(outfn, 'w')
    
20. outf.write('Lon,Lat,Time,Heigh,pres,TEMP,RELH,VORT,DIVG,es\n')
    
21. #uwnd = meteof['UWND'][:,:,:,:]
    
22. #vwnd = meteof['VWND'][:,:,:,:]
    
23. tvar = 'TEMP'
    
24. rhvar = 'RELH'
    
25. uwnd = meteof['UWND'][:,:,[10, 80],[50,140]]
    
26. vwnd = meteof['VWND'][:,:,[10, 80],[50,140]]
    
27. vort = hcurl(uwnd, vwnd)
    
28. divg = hdivg(uwnd, vwnd)
    
29. 
    
30. #######假相当位温的计算#####
    
31. es = 6.1078*exp(17.2693882*(tvar-273.16)/(tvar-35.86))
    
32. qq = rvar*(0.62197*es/(prs-0.378*es))/100.
    
33. e = prs*qq/(0.62197+qq)+1e-10
    
34. tlcl = 55.0+2840.0/(3.5*log(t0)-log(e)-4.805)
    
35. theta = t0*pow((1000/prs),(0.2854*(1.0-0.28*qq)))
    
36. eqt = theta*exp(((3376./tlcl)-2.54)*qq*(1.0+0.81*qq))
    
37. thse = eqt-273.15
    
38. print thse
    
39. ###########湿位涡###########grads程序，没有做相应修改###
    
40. f = 2*7.292*sin(lat*3.14159/180.0)*0.00001
    
41. g = 9.8
    
42. dp = 100*(lev(z-1)-lev(z+1))
    
43. deqt = eqt(z-1)-eqt(z+1)
    
44. du = uwnd(z-1)-uwnd(z+1)
    
45. dv = vwnd(z-1)-vwnd(z+1)
    
46. dx = 2.0*6370949.0*cos(lat*3.14159/180.0)*3.14159/180.0
    
47. dy = 2.0*6370949.0*3.14159/180.0
    
48. dtx = cdiff(eqt,x)
    
49. dty = cdiff(eqt,y)
    
50. pv1 = -g*(vort+f)*deqt/dp  
    
51. pv2 = g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))
    
52. pv = pv1+pv2
    
53. 
    
54. '''
    
55. idx = 0
    
56. for tline in trajLayer.shapes():   
    
57.     t = trajLayer.cellvalue('Date', idx)
    
58.     h = trajLayer.cellvalue('Hour', idx)
    
59.     t.replace(hour=h)
    
60.     for ps in tline.getPoints():  
    
61.         lon = ps.X
    
62.         lat = ps.Y
    
63.         z = ps.M
    
64.         pres = ps.Z
    
65.         v = vort.interpn([t,pres,lat,lon])
    
66.         d = divg.interpn([t,pres,lat,lon])
    
67.         temp = meteof.tostation(tvar, lon, lat, pres, t)
    
68.         rh = meteof.tostation(rhvar, lon, lat, pres, t)
    
69.         es = meteof.tostation(esvar, lon, lat, pres, t)
    
70.         print 'lon: %.3f; lat: %.3f; time: %s; height: %.1f; pres: %.1f; TEMP: %.1f;  RELH: %.1f;  VORT: %f; DIVG: %f; es:%f' % (lon,lat,t.strftime('%Y%m%d_%H:%M'),z,pres,temp,rh,v,d,es)
    
71.         line = '%.3f,%.3f,%s,%.1f,%.1f,%.1f,%.1f,%f,%f,%f' % (lon,lat,t.strftime('%Y%m%d_%H:%M'),z,pres,temp,rh,v,d,es)
    
72.         outf.write(line + '\n')
    
73.         t = t + datetime.timedelta(hours=-1)
    
74.     idx += 1
    
75. outf.close()
    
76. print 'Finish...'

复制代码

王老师，我想尝试进行湿位涡的计算和追踪，刚开始调试就出错，饱和水汽压都进行不下去!错误如下：如果王老师有时间帮忙指点！
Traceback (most recent call last):
  File "<iostream>", line 31, in <module>
TypeError: unsupported operand type(s) for -: 'str' and 'float'

MeteoInfo

良辰 发表于 2016-8-8 16:15
王老师，我想尝试进行湿位涡的计算和追踪，刚开始调试就出错，饱和水汽压都进行不下去!错误如下：如果 ...

tvar是字符串，需要提取温度数据数组参与计算。

良辰

好的，谢谢老师！

良辰

MeteoInfo 发表于 2016-8-8 16:16
tvar是字符串，需要提取温度数据数组参与计算。

1. # Set working directory
   
2. trajDir = 'D:/Temp/HYSPLIT'
   
3. meteoDir = r'D:\Temp\arl'
   
4. # Open trjactory data file
   
5. print 'Open trajectory data file ...'
   
6. trajfn = os.path.join(trajDir, 'traj_20150331')
   
7. print 'Trajectory file: ' + trajfn
   
8. trajf = addfile_hytraj(trajfn)
   
9. # Create trajectory layer
   
10. trajLayer = trajf.trajlayer()
    
11. # Open meteorological data file
    
12. print 'Open meteorological data file...'
    
13. meteofn = os.path.join(meteoDir, 'gdas1.mar15.w5')
    
14. print 'Meteorological file: ' + meteofn
    
15. meteof = addfile(meteofn)
    
16. # Get meteorological data along trajectory
    
17. print 'Get meteorological data along trajectory...'
    
18. outfn = os.path.join(trajDir, 'traj.txt')
    
19. outf = open(outfn, 'w')
    
20. outf.write('Lon,Lat,Time,Heigh,pres,TEMP,RELH,VORT,DIVG,es\n')
    
21. #uwnd = meteof['UWND'][:,:,:,:]
    
22. #vwnd = meteof['VWND'][:,:,:,:]
    
23. tvar = meteof['TEMP'][:,:,:,:]
    
24. rhvar = meteof['RELH'][:,:,:,:]
    
25. uwnd = meteof['UWND'][:,:,[10, 80],[50,160]]
    
26. vwnd = meteof['VWND'][:,:,[10, 80],[50,160]]
    
27. gv = meteof['HGTS']
    
28. nx = gv.dimlen(gv.ndim - 1)
    
29. ny = gv.dimlen(gv.ndim - 2)
    
30. levels = gv.dimvalue(gv.ndim - 3)[::1]
    
31. prs = levels
    
32. #vort = hcurl(uwnd, vwnd)
    
33. #divg = hdivg(uwnd, vwnd)
    
34. 
    
35. #######假相当位温的计算#####
    
36. es = 6.1078*exp(17.2693882*(tvar-273.16)/(tvar-35.86))
    
37. print es
    
38. '''
    
39. qq = rvar*(0.62197*es/(prs-0.378*es))/100.
    
40. e = prs*qq/(0.62197+qq)+1e-10
    
41. tlcl = 55.0+2840.0/(3.5*log(tvar)-log(e)-4.805)
    
42. theta = tvar*pow((1000/prs),(0.2854*(1.0-0.28*qq)))
    
43. eqt = theta*exp(((3376./tlcl)-2.54)*qq*(1.0+0.81*qq))
    
44. thse = eqt-273.15
    
45. print thse
    
46. ###########湿位涡###########grads程序，没有做相应修改###
    
47. f = 2*7.292*sin(lat*3.14159/180.0)*0.00001
    
48. g = 9.8
    
49. dp = 100*(lev(z-1)-lev(z+1))
    
50. deqt = eqt(z-1)-eqt(z+1)
    
51. du = uwnd(z-1)-uwnd(z+1)
    
52. dv = vwnd(z-1)-vwnd(z+1)
    
53. dx = 2.0*6370949.0*cos(lat*3.14159/180.0)*3.14159/180.0
    
54. dy = 2.0*6370949.0*3.14159/180.0
    
55. dtx = cdiff(eqt,x)
    
56. dty = cdiff(eqt,y)
    
57. pv1 = -g*(vort+f)*deqt/dp  
    
58. pv2 = g*((dv/dp)*(dtx/dx)-(du/dp)*(dty/dy))
    
59. pv = pv1+pv2
    
60. '''
    
61. idx = 0
    
62. for tline in trajLayer.shapes():   
    
63.     t = trajLayer.cellvalue('Date', idx)
    
64.     h = trajLayer.cellvalue('Hour', idx)
    
65.     t.replace(hour=h)
    
66.     for ps in tline.getPoints():  
    
67.         lon = ps.X
    
68.         lat = ps.Y
    
69.         z = ps.M
    
70.         pres = ps.Z
    
71.         v = vort.interpn([t,pres,lat,lon])
    
72.         d = divg.interpn([t,pres,lat,lon])
    
73. #        temp = meteof.tostation(tvar, lon, lat, pres, t)
    
74. #        rh = meteof.tostation(rhvar, lon, lat, pres, t)
    
75. #        es = meteof.tostation(esvar, lon, lat, pres, t)
    
76.         print 'lon: %.3f; lat: %.3f; time: %s; height: %.1f; pres: %.1f; TEMP: %.1f;  RELH: %.1f;  VORT: %f; DIVG: %f' % (lon,lat,t.strftime('%Y%m%d_%H:%M'),z,pres,temp,rh,v,d)
    
77.         line = '%.3f,%.3f,%s,%.1f,%.1f,%.1f,%.1f,%f,%f,%f' % (lon,lat,t.strftime('%Y%m%d_%H:%M'),z,pres,temp,rh,v,d)
    
78.         outf.write(line + '\n')
    
79.         t = t + datetime.timedelta(hours=-1)
    
80.     idx += 1
    
81. outf.close()
    
82. print 'Finish...'

复制代码

王老师，我想计算位涡。提取温度数据数组参与计算了，但是还是出错，错误如下：
File "<iostream>", line 36, in <module>
  File "D:\Temp\MeteoInfo\pylib\mipylib\dimarray.py", line 210, in __rmul__
    return DimArray.__mul__(self, other)
  File "D:\Temp\MeteoInfo\pylib\mipylib\dimarray.py", line 206, in __mul__
    r = DimArray(self.array.__mul__(other), self.dims, self.fill_value, self.proj)
  File "D:\Temp\MeteoInfo\pylib\mipylib\miarray.py", line 166, in __mul__
    r = MIArray(ArrayMath.mul(self.array, other))
王老师，有时间帮我看一下呗？

MeteoInfo

良辰 发表于 2016-8-16 09:14
王老师，我想计算位涡。提取温度数据数组参与计算了，但是还是出错，错误如下：
File "", line 36, in ...

参考此贴：MeteoInfoLab脚本示例：计算湿位涡
http://bbs.06climate.com/forum.p ... 931&fromuid=106
(出处: 气象家园)

需要下载最新版本的MeteoInfo。

良辰

MeteoInfo 发表于 2016-8-16 23:12
参考此贴：MeteoInfoLab脚本示例：计算湿位涡
http://bbs.06climate.com/forum.php?mod=viewthread&tid= ...

非常感谢王老师！

MeteoInfo

良辰 发表于 2016-8-17 16:16
非常感谢王老师！

你那里有条件的话最好验证一下湿位温的计算结果是否正确。

良辰

MeteoInfo 发表于 2016-8-17 16:33
你那里有条件的话最好验证一下湿位温的计算结果是否正确。

好的，我有再分析资料画一下，进行对比！我还要计算一下位涡。这几天弄好了，一起发上来！

wanGKun

哎哟，我羞涩了:$
不过，这里真是一个气象人待下去的家

zhaojunri1020

@meteoinfo 王老师 做pscf分析的时候我应该怎么求或者怎么在meteoinfo里边看所设网格内所以轨迹点数？