grib2转arl,结果错误，求王老师指点！

良辰

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['Vertical_velocity_pressure_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 = [100000,97500,95000,92500,90000,85000,80000,75000,70000,\
    
33. # 65000,60000,55000,50000,45000,40000,35000,30000,25000,20000,15000,10000]
    
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 nz
    
72. print lidx
    
73. print llidx
    
74. for avname,gvname in zip(avar3d, gvar3d):
    
75. print avname + ' ' + gvname
    
76. gdata = inf[gvname][t,llidx,:,:]
    
77. print gdata
    
78. if avname == 'WWND':
    
79. gdata = gdata * 0.01
    
80. ksum = arlf.writedatarec(atime, lidx + 1, avname, fhour, 99, gdata)
    
81. ksums.append(ksum)
    
82. ksumlist.append(ksums)
    
83. #Re-write index record with checksum
    
84. arlf.writeindexrec(atime, dhead, ksumlist)
    
85. fhour += 1
    
86. arlf.close()
    
87. print 'Finished!'

复制代码

王老师，程序脚本运行正常，三维的数据输出错误！gdata貌似有问题。有时间能帮忙看一下吗？另外，ARL数据写入有什么单位、顺序等要求吗？