MeteoInfo脚本示例：GRIB to ARL

MeteoInfo

chenxl 发表于 2015-6-30 14:57
老师，我想请教3个关于wrf-out数据错列的问题。1 wrf输出数据出现错列是由于什么导致的呢?可不可以通过wrf ...

1、错列应该是WRF本身设计成这样，是否能够通过设置修改输出我也不知道（没转过WRF）。是否错列，你自己用MeteoInfo打开数据可以看变量的维是否有错列维（有staggered字符）。

2、应该和错列关系不大，你可以看看两个数据的风场差别大不大。

chenxl

MeteoInfo 发表于 2015-6-30 15:13
1、错列应该是WRF本身设计成这样，是否能够通过设置修改输出我也不知道（没转过WRF）。是否错列，你自己 ...

谢谢老师。我自己再摸索摸索。

红色的兔子

王老师，您好，我近期工作也是用格点资料转化成arl资料
我尝试了转化了一下，发现可以转化生成arl数据，但是我对比了一下，arl资料和原始的grads资料有错误
地表资料没有问题，但是上层资料数据是错误的，与grads资料不对应

http://pan.baidu.com/s/1i39xW7z，提取码：25mb
这是我的资料，griball.grd是原数据，ctl是说明文件，test.arl是我用资料生成的文件

下面是我的脚本

#--------------------------------------------------------        
# Author: Yaqiang Wang                                          
# Date: 2014-10-24                                               
# Purpose: Convert Grads data to ARL data  
# Note: Sample                                                   
#-----------------------------------------------------------
from org.meteoinfo.data.meteodata import MeteoDataInfo
from org.meteoinfo.data.meteodata.arl import ARLDataInfo
from org.meteoinfo.data.meteodata.arl import DataLabel
import os

#---- Set directories
dataDir = "F:/shuiqi/ncep"

#---- Set output data file
outFile = os.path.join(dataDir, 'test.arl')

#---- Read a GRADS data file
mydata = MeteoDataInfo()
infile = os.path.join(dataDir, 'griball.ctl')
print infile
mydata.openGrADSData(infile)
print 'GRADS file has been opened...'


#---- Set output ARL data info
arlDI = ARLDataInfo()


#---- Set variable and level list
gvars = ['PRSS','T02M','U10M','V10M','HGTS','TEMP','UWND','VWND','WWND','RELH']
avars = ['PRSS','T02M','U10M','V10M','HGTS','TEMP','UWND','VWND','WWND','RELH']
levels = [0,1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10]
for l in levels:
  arlDI.levels.add(l)
  if l == 0:
    arlDI.LevelVarList.add(['PRSS','T02M','U10M','V10M'])
#  elif l == [1000,925,850,700,600,500,400,300]:
#    arlDI.LevelVarList.add(['HGTS','TEMP','UWND','VWND','WWND','RELH'])
#  elif l == [250,200,150,100]:
#    arlDI.LevelVarList.add(['HGTS','TEMP','UWND','VWND','WWND'])
  else:
    arlDI.LevelVarList.add(['HGTS','TEMP','UWND','VWND','WWND','RELH'])   
#---- Write ARL data file

dataInfo = mydata.getDataInfo()
arlDI.createDataFile(outFile)
arlDI.X = dataInfo.getXDimension().getValues()
arlDI.Y = dataInfo.getYDimension().getValues()
variables = dataInfo.getVariables()
tNum = dataInfo.getTimeNum()
for t in range(0, tNum):
  mydata.setTimeIndex(t)
  atime = dataInfo.getTimes().get(t)
  aDH = arlDI.getDataHead(mydata.getProjectionInfo(), 'NCEP', 2)
  arlDI.writeIndexRecord(atime, aDH)
  lidx = 0
  for l in arlDI.levels:
    print l
    for v in arlDI.LevelVarList[lidx]:
      vName = gvars[avars.index(v)]
      print vName
      if lidx == 0:
        mydata.setLevelIndex(lidx)
      else:
        variable = dataInfo.getVariable(vName)
        nlidx = variable.getZDimension().getDimValue().indexOf(l*100.0)
        mydata.setLevelIndex(nlidx)
      gData = mydata.getGridData(vName)
      if v == 'PRSS' or v == 'WWND':
        gData = gData.div(100)      
      aDL = DataLabel(atime)
      aDL.setLevel(lidx)
      aDL.setVarName(v)
      aDL.setGrid(99)
      aDL.setForecast(0)
      arlDI.writeGridData(aDL, gData)
    lidx += 1
arlDI.closeDataFile()

print 'Finished!'

不知道是什么地方出现问题了，不知道您能不能帮助修改一下。
谢谢老师了

红色的兔子

希望能尽快得到回复，谢谢啦

MeteoInfo

红色的兔子 发表于 2015-7-1 14:42
希望能尽快得到回复，谢谢啦

首先你的ctl文件中需要加上options yrev，因为你的数据是Y轴反向的，从PRSS变量的图形很容易看出来。

修改后的脚本程序如下：

1. #--------------------------------------------------------        
   
2. # Author: Yaqiang Wang                                          
   
3. # Date: 2014-10-24                                               
   
4. # Purpose: Convert Grads 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/test"
    
14. 
    
15. #---- Set output data file
    
16. outFile = os.path.join(dataDir, 'test.arl')
    
17. 
    
18. #---- Read a GRADS data file
    
19. mydata = MeteoDataInfo()
    
20. infile = os.path.join(dataDir, 'griball.ctl')
    
21. print infile
    
22. mydata.openGrADSData(infile)
    
23. print 'GRADS file has been opened...'
    
24. 
    
25. 
    
26. #---- Set output ARL data info
    
27. arlDI = ARLDataInfo()
    
28. 
    
29. 
    
30. #---- Set variable and level list
    
31. gvars = ['PRSS','T02M','U10M','V10M','HGTS','TEMP','UWND','VWND','WWND','RELH']
    
32. avars = ['PRSS','T02M','U10M','V10M','HGTS','TEMP','UWND','VWND','WWND','RELH']
    
33. levels = [0,1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10]
    
34. for l in levels:
    
35.     arlDI.levels.add(l)
    
36.     if l == 0:
    
37.         arlDI.LevelVarList.add(['PRSS','T02M','U10M','V10M'])
    
38.     elif l >= 300:
    
39.         arlDI.LevelVarList.add(['HGTS','TEMP','UWND','VWND','WWND','RELH'])
    
40.     else:
    
41.         arlDI.LevelVarList.add(['HGTS','TEMP','UWND','VWND','WWND'])
    
42.     #else:
    
43.     #    arlDI.LevelVarList.add(['HGTS','TEMP','UWND','VWND','WWND','RELH'])   
    
44. #---- Write ARL data file
    
45. 
    
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(), 'NCEP', 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*1.0)
    
68.                 if nlidx < 0:
    
69.                     print 'No level: ' + str(l)
    
70.                     continue
    
71.                 mydata.setLevelIndex(nlidx)
    
72.             gData = mydata.getGridData(vName)
    
73.             if v == 'PRSS' or v == 'WWND':
    
74.                 gData = gData.div(100)      
    
75.             aDL = DataLabel(atime)
    
76.             aDL.setLevel(lidx)
    
77.             aDL.setVarName(v)
    
78.             aDL.setGrid(99)
    
79.             aDL.setForecast(0)
    
80.             arlDI.writeGridData(aDL, gData)
    
81.         lidx += 1
    
82. arlDI.closeDataFile()
    
83. 
    
84. print 'Finished!'

#--------------------------------------------------------        <br /> # Author: Yaqiang Wang                                           <br /> # Date: 2014-10-24                                               <br /> # Purpose: Convert Grads 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/test"<br /> <br /> #---- Set output data file<br /> outFile = os.path.join(dataDir, 'test.arl')<br /> <br /> #---- Read a GRADS data file<br /> mydata = MeteoDataInfo()<br /> infile = os.path.join(dataDir, 'griball.ctl')<br /> print infile<br /> mydata.openGrADSData(infile)<br /> print 'GRADS file has been opened...'<br /> <br /> <br /> #---- Set output ARL data info<br /> arlDI = ARLDataInfo()<br /> <br /> <br /> #---- Set variable and level list<br /> gvars = ['PRSS','T02M','U10M','V10M','HGTS','TEMP','UWND','VWND','WWND','RELH']<br /> avars = ['PRSS','T02M','U10M','V10M','HGTS','TEMP','UWND','VWND','WWND','RELH']<br /> levels = [0,1000,925,850,700,600,500,400,300,250,200,150,100,70,50,30,20,10]<br /> for l in levels:<br />     arlDI.levels.add(l)<br />     if l == 0:<br />         arlDI.LevelVarList.add(['PRSS','T02M','U10M','V10M'])<br />     elif l >= 300:<br />         arlDI.LevelVarList.add(['HGTS','TEMP','UWND','VWND','WWND','RELH'])<br />     else:<br />         arlDI.LevelVarList.add(['HGTS','TEMP','UWND','VWND','WWND'])<br />     #else:<br />     #    arlDI.LevelVarList.add(['HGTS','TEMP','UWND','VWND','WWND','RELH'])    <br /> #---- Write ARL data file<br /> <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(), 'NCEP', 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*1.0)<br />                 if nlidx < 0:<br />                     print 'No level: ' + str(l)<br />                     continue<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 /> arlDI.closeDataFile()<br /> <br /> print 'Finished!'

红色的兔子

MeteoInfo 发表于 2015-7-1 22:22
首先你的ctl文件中需要加上options yrev，因为你的数据是Y轴反向的，从PRSS变量的图形很容易看出来。
...

谢谢王老师，万分感谢

lovechang1314

请问有没有输出到micaps4类格式的没，我看文档里面没有micaps4INFO类里面没有输出的？

MeteoInfo

lovechang1314 发表于 2015-7-24 17:17
请问有没有输出到micaps4类格式的没，我看文档里面没有micaps4INFO类里面没有输出的？

目前没有这样的功能。

自由行

脚本一直都是很难研究的问题，受益匪浅

红色的兔子

王老师，我最近使用您修改的脚本转化数据。之前转化成功了，但是最近我换了一组数据之后，结果又出现了错误了。我运行脚本之后，能生成arl文件，但是生成的数据用meteoinfo和hysplit软件都不能打开。

链接: http://pan.baidu.com/s/1gdJGkJD ，提取密码: f7fg
这是我的资料，griball199105.grd是原数据，griball199105.ctl是说明文件，FGOALS199105.arl是我用资料生成的文件,1.py是修改过的脚本。
不知道错误出现在哪里
希望您抽空看一下。万分感谢