GRIB格式转为ARL格式

MeteoInfo · 发表于 2013-1-2 13:35:23

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

MeteoInfo支持多种气象数据格式，比如netCDF, GRIB 1&2, ARL等，但是缺乏写ARL数据的能力，因此最近做了一些这方面的工作，有了一些初步结果，还需要进一步完善。ARL格式用HYSPLIT模式的人应该比较清楚。需要MeteoInfo最新文件（见置顶帖子）。

转换需要写脚本来完成。测试了一个FNL的GRIB2数据，转换后的ARL和原GRIB2数据是吻合的，当然还需要更多的测试。

示例脚本代码如下：

#--------------------------------------------------------
# Author: Yaqiang Wang
# Date: 2012-12-31
# Purpose: Convert GRIB data to ARL data
# Note: Sample
#-----------------------------------------------------------
import clr
from System import *
from System.Collections.Generic import *
clr.AddReference("MeteoInfoC.dll")
from MeteoInfoC import *
from MeteoInfoC.Data import *
from MeteoInfoC.Data.MeteoData import *
#---- Set directories
dataDir = "D:\\Temp\\"
#---- Set output data file
outFile = dataDir + 'arl\\test1.arl'
#---- Read a GRIB data file
mydata = MeteoDataInfo()
infile = dataDir + 'grib\\fnl_20110416_00_00'
mydata.OpenGRIBData(infile)
#---- Set output ARL data info
arlDI = ARLDataInfo()
#---- Set variable and level list
gvars = ['Pressure@surface','Temperature@height_above_ground',\
'U-component_of_wind@height_above_ground','V-component_of_wind@height_above_ground',\
'Total_precipitation@surface','Geopotential_height@pressure','Temperature@pressure',\
'U-component_of_wind@pressure','V-component_of_wind@pressure','Vertical_velocity@pressure',\
'Relative_humidity@pressure']
avars = ['PRSS','T02M','U10M','V10M','TPP6','HGTS','TEMP','UWND','VWND','WWND','RELH']
levels = [0,10,20,30,50,70,100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,\
850,900,925,950,975,1000]
for l in levels:
arlDI.levels.Add(l)
if l == 0:
arlDI.LevelVarList.Add(List[str](['PRSS','T02M','U10M','V10M','TPP6']))
elif l < 100:
arlDI.LevelVarList.Add(List[str](['HGTS','TEMP','UWND','VWND']))
else:
arlDI.LevelVarList.Add(List[str](['HGTS','TEMP','UWND','VWND','WWND','RELH']))
#---- Write ARL data file
arlDI.CreateDataFile(outFile)
arlDI.X = mydata.GetX()
arlDI.Y = mydata.GetY()
variables = mydata.GetVariables()
tNum = mydata.GetTimeNumber()
for t in range(0, tNum):
mydata.TimeIndex = t
atime = mydata.GetTime(t)
print atime.ToString("yyyy-MM-dd HH")
aDH = arlDI.GetDataHead(mydata.ProjInfo, 'FNL1', 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.LevelIndex = lidx
else:
variable = mydata.GetVariable(vName)
nlidx = variable.Levels.IndexOf(l)
mydata.LevelIndex = nlidx
gData = mydata.GetGridData(vName)
aDL = DataLabel(atime)
aDL.Level = lidx
aDL.Variable = v
aDL.Grid = 99
aDL.Forecast = 0
arlDI.WriteGridData(aDL, gData)
lidx += 1
arlDI.CloseDataFile()
print 'Finished!'

#--------------------------------------------------------         
# Author: Yaqiang Wang                                           
# Date: 2012-12-31                                                
# Purpose: Convert GRIB data to ARL data   
# Note: Sample                                                    
#----------------------------------------------------------- 
import clr 
from System import * 
from System.Collections.Generic import * 
clr.AddReference("MeteoInfoC.dll") 
from MeteoInfoC import * 
from MeteoInfoC.Data import * 
from MeteoInfoC.Data.MeteoData import * 
 
#---- Set directories 
dataDir = "D:\\Temp\\" 
 
#---- Set output data file 
outFile = dataDir + 'arl\\test1.arl' 
 
#---- Read a GRIB data file 
mydata = MeteoDataInfo() 
infile = dataDir + 'grib\\fnl_20110416_00_00' 
mydata.OpenGRIBData(infile) 
 
#---- Set output ARL data info 
arlDI = ARLDataInfo() 
 
#---- Set variable and level list 
gvars = ['Pressure@surface','Temperature@height_above_ground',\ 
    'U-component_of_wind@height_above_ground','V-component_of_wind@height_above_ground',\ 
    'Total_precipitation@surface','Geopotential_height@pressure','Temperature@pressure',\ 
    'U-component_of_wind@pressure','V-component_of_wind@pressure','Vertical_velocity@pressure',\ 
    'Relative_humidity@pressure'] 
avars = ['PRSS','T02M','U10M','V10M','TPP6','HGTS','TEMP','UWND','VWND','WWND','RELH'] 
levels = [0,10,20,30,50,70,100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,\ 
    850,900,925,950,975,1000] 
for l in levels: 
    arlDI.levels.Add(l) 
    if l == 0: 
        arlDI.LevelVarList.Add(List[str](['PRSS','T02M','U10M','V10M','TPP6'])) 
    elif l < 100: 
        arlDI.LevelVarList.Add(List[str](['HGTS','TEMP','UWND','VWND'])) 
    else: 
        arlDI.LevelVarList.Add(List[str](['HGTS','TEMP','UWND','VWND','WWND','RELH'])) 
 
#---- Write ARL data file 
arlDI.CreateDataFile(outFile) 
arlDI.X = mydata.GetX() 
arlDI.Y = mydata.GetY() 
variables = mydata.GetVariables() 
tNum = mydata.GetTimeNumber() 
for t in range(0, tNum): 
    mydata.TimeIndex = t 
    atime = mydata.GetTime(t) 
    print atime.ToString("yyyy-MM-dd HH") 
    aDH = arlDI.GetDataHead(mydata.ProjInfo, 'FNL1', 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.LevelIndex = lidx 
            else: 
                variable = mydata.GetVariable(vName) 
                nlidx = variable.Levels.IndexOf(l) 
                mydata.LevelIndex = nlidx 
            gData = mydata.GetGridData(vName) 
            aDL = DataLabel(atime) 
            aDL.Level = lidx 
            aDL.Variable = v 
            aDL.Grid = 99 
            aDL.Forecast = 0 
            arlDI.WriteGridData(aDL, gData) 
        lidx += 1 
 
arlDI.CloseDataFile() 
 
print 'Finished!'

topmad · 发表于 2013-1-2 14:03:20

王老师过节还不休息啊...

棉花糖 · 发表于 2013-1-2 14:53:07

好好向王老师学习

freekiller · 发表于 2013-1-2 19:38:01

赞，很实用的

喜欢伽菲 · 发表于 2013-1-2 20:23:33

谢谢楼主啦

MeteoInfo · 发表于 2013-1-5 13:40:40

具体在格式转化时可以先在MeteoInfo中打开GRIB数据，并查看其变量信息，选择某变量可以查看其分层信息，这些信息可以用于script程序中。ARL数据的变量信息可以在ARL官方网站上查找。

具体转换时还需要考虑GRIB和ARL数据中相应变量可能单位不同，比如此例中GRIB气压数据的单位是Pa，而ARL中是hPa，在写格点数据前需要先进行单位变换的计算，修改后的代码和图形如下：

#--------------------------------------------------------
# Author: Yaqiang Wang
# Date: 2012-12-31
# Purpose: Convert GRIB data to ARL data
# Note: Sample
#-----------------------------------------------------------
import clr
from System import *
from System.Collections.Generic import *
clr.AddReference("MeteoInfoC.dll")
from MeteoInfoC import *
from MeteoInfoC.Data import *
from MeteoInfoC.Data.MeteoData import *
#---- Set directories
dataDir = "D:\\Temp\\"
#---- Set output data file
outFile = dataDir + 'arl\\test1.arl'
#---- Read a GRIB data file
mydata = MeteoDataInfo()
infile = dataDir + 'grib\\fnl_20110416_00_00'
mydata.OpenGRIBData(infile)
#---- Set output ARL data info
arlDI = ARLDataInfo()
#---- Set variable and level list
gvars = ['Pressure@surface','Temperature@height_above_ground',\
'U-component_of_wind@height_above_ground','V-component_of_wind@height_above_ground',\
'Total_precipitation@surface','Geopotential_height@pressure','Temperature@pressure',\
'U-component_of_wind@pressure','V-component_of_wind@pressure','Vertical_velocity@pressure',\
'Relative_humidity@pressure']
avars = ['PRSS','T02M','U10M','V10M','TPP6','HGTS','TEMP','UWND','VWND','WWND','RELH']
levels = [0,10,20,30,50,70,100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,\
850,900,925,950,975,1000]
for l in levels:
arlDI.levels.Add(l)
if l == 0:
arlDI.LevelVarList.Add(List[str](['PRSS','T02M','U10M','V10M','TPP6']))
elif l < 100:
arlDI.LevelVarList.Add(List[str](['HGTS','TEMP','UWND','VWND']))
else:
arlDI.LevelVarList.Add(List[str](['HGTS','TEMP','UWND','VWND','WWND','RELH']))
#---- Write ARL data file
arlDI.CreateDataFile(outFile)
arlDI.X = mydata.GetX()
arlDI.Y = mydata.GetY()
variables = mydata.GetVariables()
tNum = mydata.GetTimeNumber()
for t in range(0, tNum):
mydata.TimeIndex = t
atime = mydata.GetTime(t)
print atime.ToString("yyyy-MM-dd HH")
aDH = arlDI.GetDataHead(mydata.ProjInfo, 'FNL1', 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.LevelIndex = lidx
else:
variable = mydata.GetVariable(vName)
nlidx = variable.Levels.IndexOf(l)
mydata.LevelIndex = nlidx
gData = mydata.GetGridData(vName)
if v == 'PRSS' or v == 'WWND':
gData = gData / 100
elif v == 'TPP6':
gData = gData / 1000
aDL = DataLabel(atime)
aDL.Level = lidx
aDL.Variable = v
aDL.Grid = 99
aDL.Forecast = 0
arlDI.WriteGridData(aDL, gData)
lidx += 1
arlDI.CloseDataFile()
print 'Finished!'

#--------------------------------------------------------         
# Author: Yaqiang Wang                                           
# Date: 2012-12-31                                                
# Purpose: Convert GRIB data to ARL data   
# Note: Sample                                                    
#----------------------------------------------------------- 
import clr 
from System import * 
from System.Collections.Generic import * 
clr.AddReference("MeteoInfoC.dll") 
from MeteoInfoC import * 
from MeteoInfoC.Data import * 
from MeteoInfoC.Data.MeteoData import * 
 
#---- Set directories 
dataDir = "D:\\Temp\\" 
 
#---- Set output data file 
outFile = dataDir + 'arl\\test1.arl' 
 
#---- Read a GRIB data file 
mydata = MeteoDataInfo() 
infile = dataDir + 'grib\\fnl_20110416_00_00' 
mydata.OpenGRIBData(infile) 
 
#---- Set output ARL data info 
arlDI = ARLDataInfo() 
 
#---- Set variable and level list 
gvars = ['Pressure@surface','Temperature@height_above_ground',\ 
    'U-component_of_wind@height_above_ground','V-component_of_wind@height_above_ground',\ 
    'Total_precipitation@surface','Geopotential_height@pressure','Temperature@pressure',\ 
    'U-component_of_wind@pressure','V-component_of_wind@pressure','Vertical_velocity@pressure',\ 
    'Relative_humidity@pressure'] 
avars = ['PRSS','T02M','U10M','V10M','TPP6','HGTS','TEMP','UWND','VWND','WWND','RELH'] 
levels = [0,10,20,30,50,70,100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,\ 
    850,900,925,950,975,1000] 
for l in levels: 
    arlDI.levels.Add(l) 
    if l == 0: 
        arlDI.LevelVarList.Add(List[str](['PRSS','T02M','U10M','V10M','TPP6'])) 
    elif l < 100: 
        arlDI.LevelVarList.Add(List[str](['HGTS','TEMP','UWND','VWND'])) 
    else: 
        arlDI.LevelVarList.Add(List[str](['HGTS','TEMP','UWND','VWND','WWND','RELH'])) 
 
#---- Write ARL data file 
arlDI.CreateDataFile(outFile) 
arlDI.X = mydata.GetX() 
arlDI.Y = mydata.GetY() 
variables = mydata.GetVariables() 
tNum = mydata.GetTimeNumber() 
for t in range(0, tNum): 
    mydata.TimeIndex = t 
    atime = mydata.GetTime(t) 
    print atime.ToString("yyyy-MM-dd HH") 
    aDH = arlDI.GetDataHead(mydata.ProjInfo, 'FNL1', 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.LevelIndex = lidx 
            else: 
                variable = mydata.GetVariable(vName) 
                nlidx = variable.Levels.IndexOf(l) 
                mydata.LevelIndex = nlidx 
            gData = mydata.GetGridData(vName) 
            if v == 'PRSS' or v == 'WWND': 
                gData = gData / 100 
            elif v == 'TPP6': 
                gData = gData / 1000 
            aDL = DataLabel(atime) 
            aDL.Level = lidx 
            aDL.Variable = v 
            aDL.Grid = 99 
            aDL.Forecast = 0 
            arlDI.WriteGridData(aDL, gData) 
        lidx += 1 
 
arlDI.CloseDataFile() 
 
print 'Finished!'

dolces · 发表于 2013-1-8 18:37:05

王老师，提取ARL数据写成新的ARL数据应该也差不多吧？直接改这个script行不行？

比如说我想提取GDAP0P5数据的单数层，然后再写成一个新的ARL文件，是不是在levels=[]里把需要的层次列出来写到新的文件里就可以了？

还有就是GDAS0P5文件垂直方向是气压-sigma混合坐标，在读取的时候是不是还是跟读气压坐标的数据是一样的？

MeteoInfo · 发表于 2013-1-8 20:04:25

dolces 发表于 2013-1-8 18:37

登录/注册后可看大图

王老师，提取ARL数据写成新的ARL数据应该也差不多吧？直接改这个script行不行？

比如说我想提取GDAP0P5数 ...

应该是可以的，你需要试试。注意Level循环时的处理。

如果原数据是混合坐标，那么提取的数据也是一样的。改变垂直坐标需要垂直插值，目前MeteoInfo里还没有现成的函数，在脚本里做应该也是可以的，不过估计比较繁琐。

dolces · 发表于 2013-1-14 20:31:03

MeteoInfo 发表于 2013-1-8 20:04

登录/注册后可看大图

应该是可以的，你需要试试。注意Level循环时的处理。

如果原数据是混合坐标，那么提取的数据也是一样的 ...

王老师，我改好了script，运行也没问题，可是生成的ARL文件meteoinfo不能打开，提示‘input string was not a correct format’，不知道是什么原因？

另外，改写ARL数据还有个问题，ARL数据的地面层标题不是数字，是用Surface来标记的，这种情况怎么在level循环里把地面层的数据写进新的ARL文件里？

附件里是我的script，麻烦您帮忙看一下。

arlextract.py (3.71 KB, 下载次数: 6)

MeteoInfo · 发表于 2013-1-15 09:42:28

dolces 发表于 2013-1-14 20:31

登录/注册后可看大图

王老师，我改好了script，运行也没问题，可是生成的ARL文件meteoinfo不能打开，提示‘input string was n ...

我下载一个0.5度的数据试试吧。

		自动登录	找回密码
密码			立即注册