MeteoInfo脚本示例：GRIB to ARL

MeteoInfo

之前发过帖子提供了将GRIB数据转为ARL格式数据（HYSPLIT模式所需气象数据格式）的IronPython脚本程序（MeteoInfo C#版），这里再给出一个Jython示例脚本（MeteoInfo Java版）：

1. #--------------------------------------------------------        
   
2. # Author: Yaqiang Wang                                          
   
3. # Date: 2014-10-24                                               
   
4. # Purpose: Convert GRIB data to ARL data  
   
5. # Note: Sample                                                   
   
6. #-----------------------------------------------------------
   
7. from org.meteoinfo.data.meteodata import MeteoDataInfo
   
8. from org.meteoinfo.data.meteodata.arl import ARLDataInfo
   
9. from org.meteoinfo.data.meteodata.arl import DataLabel
   
10. import os
    
11. 
    
12. #---- Set directories
    
13. dataDir = "D:/Temp"
    
14. 
    
15. #---- Set output data file
    
16. outFile = os.path.join(dataDir, 'arl/test2.arl')
    
17. 
    
18. #---- Read a GRIB data file
    
19. mydata = MeteoDataInfo()
    
20. infile = os.path.join(dataDir, 'grib/201001011800.pgbh06.gdas.20100101-20100105.grb2')
    
21. print infile
    
22. mydata.openNetCDFData(infile)
    
23. print 'GRIB file has been opened...'
    
24. 
    
25. #---- Set output ARL data info
    
26. arlDI = ARLDataInfo()
    
27. 
    
28. #---- Set variable and level list
    
29. gvars = ['Pressure_surface','Temperature_height_above_ground',\
    
30.   'u-component_of_wind_height_above_ground','v-component_of_wind_height_above_ground',\
    
31.   'Geopotential_height_isobaric','Temperature_isobaric',\
    
32.   'u-component_of_wind_isobaric','v-component_of_wind_isobaric','Vertical_velocity_pressure_isobaric',\
    
33.   'Relative_humidity_isobaric']
    
34. avars = ['PRSS','T02M','U10M','V10M','HGTS','TEMP','UWND','VWND','WWND','RELH']
    
35. levels = [0,1000,975,950,925,900,875,850,825,800,775,750,700,\
    
36.   650,600,550,500,450,400,350,300,250,225,200,175,150,\
    
37.   125,100,70,50,30,20,10,7,5,3,2,1]
    
38. for l in levels:
    
39.   arlDI.levels.add(l)
    
40.   if l == 0:
    
41.     arlDI.LevelVarList.add(['PRSS','T02M','U10M','V10M'])
    
42.   else:
    
43.     arlDI.LevelVarList.add(['HGTS','TEMP','UWND','VWND','WWND','RELH'])
    
44. 
    
45. #---- Write ARL data file
    
46. dataInfo = mydata.getDataInfo()
    
47. arlDI.createDataFile(outFile)
    
48. arlDI.X = dataInfo.getXDimension().getValues()
    
49. arlDI.Y = dataInfo.getYDimension().getValues()
    
50. variables = dataInfo.getVariables()
    
51. tNum = dataInfo.getTimeNum()
    
52. for t in range(0, tNum):
    
53.   mydata.setTimeIndex(t)
    
54.   atime = dataInfo.getTimes().get(t)
    
55.   aDH = arlDI.getDataHead(mydata.getProjectionInfo(), 'FNL1', 2)
    
56.   arlDI.writeIndexRecord(atime, aDH)
    
57.   lidx = 0
    
58.   for l in arlDI.levels:
    
59.     print l
    
60.     for v in arlDI.LevelVarList[lidx]:
    
61.       vName = gvars[avars.index(v)]
    
62.       print vName
    
63.       if lidx == 0:
    
64.         mydata.setLevelIndex(lidx)
    
65.       else:
    
66.         variable = dataInfo.getVariable(vName)
    
67.         nlidx = variable.getZDimension().getDimValue().indexOf(l*100.0)
    
68.         mydata.setLevelIndex(nlidx)
    
69.       gData = mydata.getGridData(vName)
    
70.       if v == 'PRSS' or v == 'WWND':
    
71.         gData = gData.div(100)      
    
72.       aDL = DataLabel(atime)
    
73.       aDL.setLevel(lidx)
    
74.       aDL.setVarName(v)
    
75.       aDL.setGrid(99)
    
76.       aDL.setForecast(0)
    
77.       arlDI.writeGridData(aDL, gData)
    
78.     lidx += 1
    
79. 
    
80. arlDI.closeDataFile()
    
81. 
    
82. print 'Finished!'

#--------------------------------------------------------        <br /> # Author: Yaqiang Wang                                           <br /> # Date: 2014-10-24                                               <br /> # Purpose: Convert GRIB data to ARL data  <br /> # Note: Sample                                                   <br /> #-----------------------------------------------------------<br /> from org.meteoinfo.data.meteodata import MeteoDataInfo<br /> from org.meteoinfo.data.meteodata.arl import ARLDataInfo<br /> from org.meteoinfo.data.meteodata.arl import DataLabel<br /> import os<br /> <br /> #---- Set directories<br /> dataDir = "D:/Temp"<br /> <br /> #---- Set output data file<br /> outFile = os.path.join(dataDir, 'arl/test2.arl')<br /> <br /> #---- Read a GRIB data file<br /> mydata = MeteoDataInfo()<br /> infile = os.path.join(dataDir, 'grib/201001011800.pgbh06.gdas.20100101-20100105.grb2')<br /> print infile<br /> mydata.openNetCDFData(infile)<br /> print 'GRIB file has been opened...'<br /> <br /> #---- Set output ARL data info<br /> arlDI = ARLDataInfo()<br /> <br /> #---- Set variable and level list<br /> gvars = ['Pressure_surface','Temperature_height_above_ground',\<br />   'u-component_of_wind_height_above_ground','v-component_of_wind_height_above_ground',\<br />   'Geopotential_height_isobaric','Temperature_isobaric',\<br />   'u-component_of_wind_isobaric','v-component_of_wind_isobaric','Vertical_velocity_pressure_isobaric',\<br />   'Relative_humidity_isobaric']<br /> avars = ['PRSS','T02M','U10M','V10M','HGTS','TEMP','UWND','VWND','WWND','RELH']<br /> levels = [0,1000,975,950,925,900,875,850,825,800,775,750,700,\<br />   650,600,550,500,450,400,350,300,250,225,200,175,150,\<br />   125,100,70,50,30,20,10,7,5,3,2,1]<br /> for l in levels:<br />   arlDI.levels.add(l)<br />   if l == 0:<br />     arlDI.LevelVarList.add(['PRSS','T02M','U10M','V10M'])<br />   else:<br />     arlDI.LevelVarList.add(['HGTS','TEMP','UWND','VWND','WWND','RELH'])<br /> <br /> #---- Write ARL data file<br /> dataInfo = mydata.getDataInfo()<br /> arlDI.createDataFile(outFile)<br /> arlDI.X = dataInfo.getXDimension().getValues()<br /> arlDI.Y = dataInfo.getYDimension().getValues()<br /> variables = dataInfo.getVariables()<br /> tNum = dataInfo.getTimeNum()<br /> for t in range(0, tNum):<br />   mydata.setTimeIndex(t)<br />   atime = dataInfo.getTimes().get(t)<br />   aDH = arlDI.getDataHead(mydata.getProjectionInfo(), 'FNL1', 2)<br />   arlDI.writeIndexRecord(atime, aDH)<br />   lidx = 0<br />   for l in arlDI.levels:<br />     print l<br />     for v in arlDI.LevelVarList[lidx]:<br />       vName = gvars[avars.index(v)]<br />       print vName<br />       if lidx == 0:<br />         mydata.setLevelIndex(lidx)<br />       else:<br />         variable = dataInfo.getVariable(vName)<br />         nlidx = variable.getZDimension().getDimValue().indexOf(l*100.0)<br />         mydata.setLevelIndex(nlidx)<br />       gData = mydata.getGridData(vName)<br />       if v == 'PRSS' or v == 'WWND':<br />         gData = gData.div(100)      <br />       aDL = DataLabel(atime)<br />       aDL.setLevel(lidx)<br />       aDL.setVarName(v)<br />       aDL.setGrid(99)<br />       aDL.setForecast(0)<br />       arlDI.writeGridData(aDL, gData)<br />     lidx += 1<br /> <br /> arlDI.closeDataFile()<br /> <br /> print 'Finished!'

llg2002

1、如果输入格式是二进制，并且数据是单要素单层次的多个数据写成一个arl文件，是不是需要读二进制的函数，并且进行多次循环？
2、二进制数据的时间要自己赋给变量？不像grib数据通过函数就可以获取？
3、二进制数据的经纬度是不是也要赋给变量？
4、aDL.setGrid(99)是不是格数为99，这参数是不是可以修改？
5、aDL.setForecast(0)指的是预报时效？可修改？
6、aDH = arlDI.getDataHead(mydata.getProjectionInfo(), 'FNL1', 2)的2是什么意思？可修改吗？

llg2002

3、二进制数据的经纬度间隔是不是也要赋给变量？

llg2002

Sample Surface Level Parameters中的Total precipitation (6-h) m TPP6要素，如果资料的时间分辨率为3小时，是不是降水要素为3小时累积量，要素名改为TPP3？

MeteoInfo

llg2002 发表于 2015-1-13 10:43
1、如果输入格式是二进制，并且数据是单要素单层次的多个数据写成一个arl文件，是不是需要读二进制的函数， ...

你能提供一个（或几个）二进制数据的示例文件吗？当然还要有数据格式的说明。数据本身是否包含经纬度、高度、时间等信息？

MeteoInfo

llg2002 发表于 2015-1-13 11:01
Sample Surface Level Parameters中的Total precipitation (6-h) m TPP6要素，如果资料的时间分辨率为3小时 ...

1、2、3：是的
4、这个Grid应该没有什么作用，不是格点数，而是格点的一个所谓的符号。
5、是预报时次，不过也不重要，重要的是数据的时间。
6、这个是指的垂直分层的类型：
   1 - sigma
    2 - pressure
    3 - terrain
    4 - hybrid
7、需要知道数据的经纬度（或者投影的详细信息）
8、没见过TPP3这样的变量，不知道模式认不认，可以累加成TPP6

llg2002

8、没见过TPP3这样的变量，不知道模式认不认，可以累加成TPP6

如果数据的时间分辨率是小于6小时，降水要素的分辨率仍为6小时，那就会出现有些时次没有降水要素。

MeteoInfo

llg2002 发表于 2015-1-13 14:57
如果数据的时间分辨率是小于6小时，降水要素的分辨率仍为6小时，那就会出现有些时次没有降水要素。

这个问题你可以去HYSPLIT官方论坛上去问一下：https://hysplitbbs.arl.noaa.gov/index.php

MeteoInfo

llg2002 发表于 2015-1-13 14:57
如果数据的时间分辨率是小于6小时，降水要素的分辨率仍为6小时，那就会出现有些时次没有降水要素。

我刚上HYSPLIT官网上看了一下，TPP3是可以的，EDAS40数据里就有这个变量：http://ready.arl.noaa.gov/edas40.php

llg2002

你能提供一个（或几个）二进制数据的示例文件吗？当然还要有数据格式的说明。数据本身是否包含经纬度、高度、时间等信息？

配上ctl，用aDataInfo = GrADSDataInfo()，aDataInfo.ReadDataInfo(infile)应该能读取数据