grib2转arl过程中遇到问题，求指点！

良辰

MeteoInfo 发表于 2016-8-1 15:52
那只能在结果文件中写21层了。

1. #--------------------------------------------------------
   
2. # Author: Yaqiang Wang
   
3. # Date: 2016-2-4
   
4. # Purpose: Convert GRIB data to ARL data
   
5. # Note: Sample
   
6. #-----------------------------------------------------------
   
7. #---- Set data folder
   
8. datadir = 'D:/Temp'
   
9. #---- Set output data file
   
10. outfn = os.path.join(datadir, 'arl/fnl_20150218_06_00.grib2.arl')
    
11. #if os.path.exists(outfn):
    
12. #    os.remove(outfn)
    
13. #---- Read a GRIB data file
    
14. infn = os.path.join(datadir, 'grib/fnl_20150218_06_00.grib2')
    
15. print infn
    
16. inf = addfile(infn)
    
17. print 'GRIB data file has been opened...'
    
18. #---- Set output ARL data file
    
19. arlf = addfile(outfn, 'c', dtype='arl')
    
20. #---- Set variable and level list
    
21. gvar2d = ['Pressure_surface','Temperature_surface','u-component_of_wind_height_above_ground',\
    
22.     'v-component_of_wind_height_above_ground']
    
23. gvar3d = ['Geopotential_height_isobaric','Temperature_isobaric',\
    
24.     'u-component_of_wind_isobaric','v-component_of_wind_isobaric',\
    
25.     'Vertical_velocity_pressure_isobaric','Relative_humidity_isobaric']
    
26. avar2d = ['PRSS','T02M','U10M','V10M']
    
27. avar3d = ['HGTS','TEMP','UWND','VWND','WWND','RELH']
    
28. gv = inf['Geopotential_height_isobaric']
    
29. nx = gv.dimlen(gv.ndim - 1)
    
30. ny = gv.dimlen(gv.ndim - 2)
    
31. #levels = gv.dimvalue(gv.ndim - 3)[::-1]
    
32. levels = [1000,975,950,925,900,850,800,750,700,\
    
33.   650,600,550,500,450,400,350,300,250,200,150,100]
    
34. nz = len(levels)
    
35. arlf.setlevels(levels)
    
36. arlf.set2dvar(avar2d)
    
37. for l in levels:
    
38.     arlf.set3dvar(avar3d)
    
39. 
    
40. #---- Write ARL data file
    
41. arlf.setx(gv.dimvalue(gv.ndim - 1))
    
42. arlf.sety(gv.dimvalue(gv.ndim - 2))
    
43. tNum = inf.timenum()
    
44. fhour = 0
    
45. for t in range(0, tNum):
    
46.     print 'Time index: ' + str(t)
    
47.     atime = inf.gettime(t)
    
48.     print atime.strftime('%Y-%m-%d %H:00')
    
49.     dhead = arlf.getdatahead(inf.proj, 'RSMC', 2, fhour)
    
50.     #Pre-write index record without checksum - will be over-write latter
    
51.     arlf.writeindexrec(atime, dhead)
    
52.     #Checksum list
    
53.     ksumlist = []
    
54.     # Write 2d variables
    
55.     ksums = []
    
56.     for avname,gvname in zip(avar2d, gvar2d):
    
57.         print avname + ' ' + gvname
    
58.         if avname == 'U10M' or avname == 'V10M':
    
59.             gdata = inf[gvname][t,0,:,:]
    
60.         else:
    
61.             gdata = inf[gvname][t,:,:]
    
62.         if avname == 'PRSS':
    
63.             gdata = gdata * 0.01
    
64.         ksum = arlf.writedatarec(atime, 0, avname, fhour, 99, gdata)
    
65.         ksums.append(ksum)
    
66.     ksumlist.append(ksums)
    
67.     # Write 3d variables
    
68.     for lidx in range(0, nz):
    
69.         ksums = []
    
70.         llidx = nz - lidx - 1
    
71.         print lidx
    
72.         print llidx
    
73.         for avname,gvname in zip(avar3d, gvar3d):
    
74.             print avname + ' ' + gvname
    
75.             gdata = inf[gvname][t,llidx,:,:]   
    
76.             if avname == 'WWND':
    
77.                gdata = gdata * 0.01      
    
78.             ksum = arlf.writedatarec(atime, lidx + 1, avname, fhour, 99, gdata)           
    
79.             ksums.append(ksum)
    
80.         ksumlist.append(ksums)
    
81.     #Re-write index record with checksum
    
82.     arlf.writeindexrec(atime, dhead, ksumlist)
    
83.     fhour += 1
    
84. arlf.close()
    
85. print 'Finished!'

复制代码

王老师，我又被虐了！计算出来后三维变量都是一样的！哪个循环错了！呜呜。写入的都是V10M

良辰

良辰 发表于 2016-8-1 20:30
王老师，我又被虐了！计算出来后三维变量都是一样的！哪个循环错了！呜呜。写入的都是V10M

王老师，您能帮我调试一下吗？fnl资料grib2格式的我放在云盘了。两个时次的，希望一起转ARL。我这调试和过关一样，您指点我过一关，总是麻烦您！跪求帮忙调试！
链接: http://pan.baidu.com/s/1hrL8yS8 密码: 9way

良辰

良辰

MeteoInfo 发表于 2016-8-1 15:52
那只能在结果文件中写21层了。

王老师，您能帮我调试一下吗？fnl资料grib2格式的我放在云盘了。两个时次的，希望一起转ARL。我这调试和过关一样，您指点我过一关，总是麻烦您！跪求帮忙调试！
链接: http://pan.baidu.com/s/1hrL8yS8 密码: 9way