MeteoInfoLab脚本示例：wrfout转arl

MeteoInfo

HYSPLIT模式需要特定格式（ARL）的气象数据来驱动，这里给出了将wrf模式输出的netCDf格式数据转为ARL格式的示例脚本。ARL格式中变量的名称和单位都是固定的，很多变量需要经过单位转换再写入ARL数据文件中。脚本里有详细的注释。
脚本程序：

1. #--------------------------------------------------------        # Author: Yaqiang Wang                                          
   
2. # Date: 2015-12-2                                            
   
3. # Purpose: Convert WRF out netCDF data to ARL data  
   
4. # Note: Sample                                                   
   
5. #-----------------------------------------------------------
   
6. #---- Set data folder
   
7. datadir = 'E:/Temp/Yaqiang'
   
8. #---- Set output data file
   
9. outfn = os.path.join(datadir, 'test_wrfout1.arl')
   
10. if os.path.exists(outfn):
    
11.     os.remove(outfn)
    
12. #---- Read a netCDF data file
    
13. infn = os.path.join(datadir, 'wrfout_d01_2009-06-08_00_00_00')
    
14. print infn
    
15. inf = addfile(infn)
    
16. print 'NetCDF file has been opened...'
    
17. #---- Set output ARL data file
    
18. arlf = addfile(outfn, 'c', dtype='arl')
    
19. #---- Set variable and level list
    
20. wvar2d = ['HGT','PSFC','PBLH','UST','SWDOWN','HFX','LH','T2','U10','V10','RAINNC']
    
21. wvar3d = ['P','T','U','V','W','QVAPOR']
    
22. avar2d = ['SHGT','PRSS','PBLH','USTR','DSWF','SHTF','LHTF','T02M','U10M','V10M','TPPA']
    
23. avar3d = ['PRES','TEMP','UWND','VWND','WWND','SPHU']
    
24. wv = inf['P']
    
25. nx = wv.dimlen(wv.ndim - 1)
    
26. ny = wv.dimlen(wv.ndim - 2)
    
27. levels = wv.dimvalue(wv.ndim - 3)
    
28. nz = len(levels)
    
29. arlf.setlevels(levels)
    
30. arlf.set2dvar(avar2d)
    
31. for l in levels:
    
32.     arlf.set3dvar(avar3d)
    
33. #---- Constant for poisson's equation
    
34. cp = 1004.0         # J/kg/K; specific heat
    
35. grav = 9.81         # m/s**2; gravity
    
36. rdry = 287.04       # J/kg/K; gas constant
    
37. rovcp = rdry / cp   # constant for poisson's equation
    
38. #---- Write ARL data file
    
39. arlf.setx(wv.dimvalue(wv.ndim - 1))
    
40. arlf.sety(wv.dimvalue(wv.ndim - 2))
    
41. tNum = inf.timenum()
    
42. fhour = 0
    
43. for t in range(0, tNum):
    
44.     print 'Time index: ' + str(t)
    
45.     atime = inf.gettime(t)   
    
46.     print atime.strftime('%Y-%m-%d %H:00')
    
47.     dhead = arlf.getdatahead(inf.proj, 'AWRF', 1, fhour)  
    
48.     #Pre-write index record without checksum - will be over-write latter
    
49.     arlf.writeindexrec(atime, dhead)
    
50.     #Checksum list
    
51.     ksumlist = []
    
52.     # Write 2d variables
    
53.     ksums = []
    
54.     for avname,wvname in zip(avar2d, wvar2d):        
    
55.         #print avname + ' ' + wvname
    
56.         gdata = inf[wvname][t,:,:]
    
57.         if avname == 'PRSS':
    
58.             gdata = gdata * 0.01
    
59.         elif avname == 'TPPA':
    
60.             gdata = gdata * 0.001
    
61.         ksum = arlf.writedatarec(atime, 0, avname, fhour, 99, gdata)
    
62.         ksums.append(ksum)
    
63.     ksumlist.append(ksums)
    
64.     # Write 3d variables
    
65.     for lidx in range(0, nz):
    
66.         ksums = []
    
67.         #print lidx
    
68.         pp = inf['P'][t,lidx,:,:]
    
69.         pb = inf['PB'][t,lidx,:,:]
    
70.         pres = pp + pb        
    
71.         uwnd = inf['U'][t,lidx,:,:]               
    
72.         vwnd = inf['V'][t,lidx,:,:]        
    
73.         temp = inf['T'][t,lidx,:,:]
    
74.         #potential to ambient temperature
    
75.         temp = (temp + 300.) * (pres / 100000.) ** rovcp        
    
76.         sphu = inf['QVAPOR'][t,lidx,:,:]        
    
77.         wwnd = inf['W'][t,lidx+1,:,:]
    
78.         #convert vertical velocity from m/s to hPa/s using omega = -W g rho
    
79.         wwnd = -wwnd * grav * pres * 0.01 / (rdry * temp * (1.0 + 0.6077 * sphu))
    
80. 
    
81.         pres = pres * 0.01
    
82.         ksum = arlf.writedatarec(atime, lidx + 1, 'PRES', fhour, 99, pres)
    
83.         ksums.append(ksum)
    
84.         ksum = arlf.writedatarec(atime, lidx + 1, 'TEMP', fhour, 99, temp)
    
85.         ksums.append(ksum)
    
86.         ksum = arlf.writedatarec(atime, lidx + 1, 'UWND', fhour, 99, uwnd[:,:nx])
    
87.         ksums.append(ksum)
    
88.         ksum = arlf.writedatarec(atime, lidx + 1, 'VWND', fhour, 99, vwnd[:ny,:])
    
89.         ksums.append(ksum)
    
90.         ksum = arlf.writedatarec(atime, lidx + 1, 'WWND', fhour, 99, wwnd)
    
91.         ksums.append(ksum)      
    
92.         ksum = arlf.writedatarec(atime, lidx + 1, 'SPHU', fhour, 99, sphu)
    
93.         ksums.append(ksum)
    
94.         ksumlist.append(ksums)
    
95.     #Re-write index record with checksum
    
96.     arlf.writeindexrec(atime, dhead, ksumlist)
    
97.     fhour += 6
    
98. arlf.close()
    
99. print 'Finished!'

#--------------------------------------------------------        # Author: Yaqiang Wang                                           <br /> # Date: 2015-12-2                                            <br /> # Purpose: Convert WRF out netCDF data to ARL data  <br /> # Note: Sample                                                   <br /> #-----------------------------------------------------------<br /> #---- Set data folder<br /> datadir = 'E:/Temp/Yaqiang'<br /> #---- Set output data file<br /> outfn = os.path.join(datadir, 'test_wrfout1.arl')<br /> if os.path.exists(outfn):<br />     os.remove(outfn)<br /> #---- Read a netCDF data file<br /> infn = os.path.join(datadir, 'wrfout_d01_2009-06-08_00_00_00')<br /> print infn<br /> inf = addfile(infn)<br /> print 'NetCDF file has been opened...'<br /> #---- Set output ARL data file<br /> arlf = addfile(outfn, 'c', dtype='arl')<br /> #---- Set variable and level list<br /> wvar2d = ['HGT','PSFC','PBLH','UST','SWDOWN','HFX','LH','T2','U10','V10','RAINNC']<br /> wvar3d = ['P','T','U','V','W','QVAPOR']<br /> avar2d = ['SHGT','PRSS','PBLH','USTR','DSWF','SHTF','LHTF','T02M','U10M','V10M','TPPA']<br /> avar3d = ['PRES','TEMP','UWND','VWND','WWND','SPHU']<br /> wv = inf['P']<br /> nx = wv.dimlen(wv.ndim - 1)<br /> ny = wv.dimlen(wv.ndim - 2)<br /> levels = wv.dimvalue(wv.ndim - 3)<br /> nz = len(levels)<br /> arlf.setlevels(levels)<br /> arlf.set2dvar(avar2d)<br /> for l in levels:<br />     arlf.set3dvar(avar3d)<br /> #---- Constant for poisson's equation<br /> cp = 1004.0         # J/kg/K; specific heat<br /> grav = 9.81         # m/s**2; gravity<br /> rdry = 287.04       # J/kg/K; gas constant<br /> rovcp = rdry / cp   # constant for poisson's equation<br /> #---- Write ARL data file<br /> arlf.setx(wv.dimvalue(wv.ndim - 1))<br /> arlf.sety(wv.dimvalue(wv.ndim - 2))<br /> tNum = inf.timenum()<br /> fhour = 0<br /> for t in range(0, tNum):<br />     print 'Time index: ' + str(t)<br />     atime = inf.gettime(t)   <br />     print atime.strftime('%Y-%m-%d %H:00') <br />     dhead = arlf.getdatahead(inf.proj, 'AWRF', 1, fhour)  <br />     #Pre-write index record without checksum - will be over-write latter<br />     arlf.writeindexrec(atime, dhead)<br />     #Checksum list<br />     ksumlist = []<br />     # Write 2d variables<br />     ksums = []<br />     for avname,wvname in zip(avar2d, wvar2d):        <br />         #print avname + ' ' + wvname<br />         gdata = inf[wvname][t,:,:]<br />         if avname == 'PRSS':<br />             gdata = gdata * 0.01<br />         elif avname == 'TPPA':<br />             gdata = gdata * 0.001<br />         ksum = arlf.writedatarec(atime, 0, avname, fhour, 99, gdata)<br />         ksums.append(ksum)<br />     ksumlist.append(ksums)<br />     # Write 3d variables<br />     for lidx in range(0, nz):<br />         ksums = []<br />         #print lidx<br />         pp = inf['P'][t,lidx,:,:]<br />         pb = inf['PB'][t,lidx,:,:]<br />         pres = pp + pb        <br />         uwnd = inf['U'][t,lidx,:,:]               <br />         vwnd = inf['V'][t,lidx,:,:]        <br />         temp = inf['T'][t,lidx,:,:]<br />         #potential to ambient temperature<br />         temp = (temp + 300.) * (pres / 100000.) ** rovcp        <br />         sphu = inf['QVAPOR'][t,lidx,:,:]        <br />         wwnd = inf['W'][t,lidx+1,:,:]<br />         #convert vertical velocity from m/s to hPa/s using omega = -W g rho<br />         wwnd = -wwnd * grav * pres * 0.01 / (rdry * temp * (1.0 + 0.6077 * sphu))<br /> <br />         pres = pres * 0.01<br />         ksum = arlf.writedatarec(atime, lidx + 1, 'PRES', fhour, 99, pres)<br />         ksums.append(ksum)<br />         ksum = arlf.writedatarec(atime, lidx + 1, 'TEMP', fhour, 99, temp)<br />         ksums.append(ksum)<br />         ksum = arlf.writedatarec(atime, lidx + 1, 'UWND', fhour, 99, uwnd[:,:nx])<br />         ksums.append(ksum)<br />         ksum = arlf.writedatarec(atime, lidx + 1, 'VWND', fhour, 99, vwnd[:ny,:])<br />         ksums.append(ksum)<br />         ksum = arlf.writedatarec(atime, lidx + 1, 'WWND', fhour, 99, wwnd)<br />         ksums.append(ksum)       <br />         ksum = arlf.writedatarec(atime, lidx + 1, 'SPHU', fhour, 99, sphu)<br />         ksums.append(ksum)<br />         ksumlist.append(ksums)<br />     #Re-write index record with checksum<br />     arlf.writeindexrec(atime, dhead, ksumlist)<br />     fhour += 6<br /> arlf.close()<br /> print 'Finished!'

转换后的数据和HYSPLIT官方提供的arw2arl程序转换的结果是一致的。

梁的丰

{:5_213:}{:5_213:}{:5_213:}{:5_213:}{:5_213:}{:5_213:}{:5_213:}

river

赞一个，太棒了！

○畏，奶牛

王老师您好，我重新安装了官网上的最新版本（MeteoInfo 1.2.8R15 - Java (for all systems)），然后按照你的脚本示例进行了修改，但是报错（见图片），感觉很奇怪啊

hbm

谢谢王老师，又学到了新的知识

MeteoInfo

○畏，奶牛 发表于 2015-12-2 16:32
王老师您好，我重新安装了官网上的最新版本（MeteoInfo 1.2.8R15 - Java (for all systems)），然后按照你 ...

在MeteoInfoLab中运行。参考此贴：http://bbs.06climate.com/forum.p ... &extra=page%3D1

MeteoInfo

○畏，奶牛 发表于 2015-12-2 16:32
王老师您好，我重新安装了官网上的最新版本（MeteoInfo 1.2.8R15 - Java (for all systems)），然后按照你 ...

在wrf转arl功能中有一个bug，转多时次数据会有问题，需要重新下载MeteoInfo最新版本，另外1楼脚本也有更新，包含了更多的变量。

aas314

好贴好东西

○畏，奶牛

MeteoInfo 发表于 2015-12-3 23:58
在wrf转arl功能中有一个bug，转多时次数据会有问题，需要重新下载MeteoInfo最新版本，另外1楼脚本也有更 ...

哦哦，谢谢王老师的提醒！昨天听您的换到在MeteoInfoLab中运行就不报错了，但是确实有时次的问题。今天又下载了官网的最新版本：MeteoInfo 1.2.8R16 - Java (for all systems)，对照您1楼更新的脚本进行了修改，并在MeteoInfoLab中运行，没有报错（见图1）。
但是当我运行hysplit的时候在Run Model这一步卡住了（见图2），我在Setup Run设置的是和以前一样的。然后我用MeteoInfo想打开转换的数据看一下，发现也不能打开，有报错（见图3）。看图1好像转换成功了，但是还是有错吗？
文件链接: http://pan.baidu.com/s/1A3xp0 密码: w86j

○畏，奶牛

图片的顺序还是错了，顺序变成了：图2，图3，图1