Convert GRIB data to ARL dataΒΆ

ARL meteorological data format is specified using in HYSPLIT model. This is an example script for converting GRIB data to ARL data.

# Convert GRIB data to ARL data

#---- Set data folder
datadir = 'D:/Temp/grib'
#---- Set output data file
outfn = os.path.join(datadir, 'test_grib.arl')
#if os.path.exists(outfn):
#    os.remove(outfn)
#---- Read a GRIB data file
infn = os.path.join(datadir, 'INGVC_2802_00_024')
print infn
inf = addfile(infn)
print 'GRIB data file has been opened...'
#---- Set output ARL data file
arlf = addfile(outfn, 'c', dtype='arl')
#---- Set variable and level list
gvar2d = ['Pressure_surface','Temperature_surface','u-component_of_wind_height_above_ground',\
    'v-component_of_wind_height_above_ground']
gvar3d = ['Geopotential_isobaric','Temperature_isobaric','Pressure_Vertical_velocity_isobaric',\
    'u-component_of_wind_isobaric','v-component_of_wind_isobaric','Specific_humidity_isobaric']
avar2d = ['PRSS','T02M','U10M','V10M']
avar3d = ['HGTS','TEMP','WWND','UWND','VWND','SPHU']
gv = inf['Geopotential_isobaric']
nx = gv.dimlen(gv.ndim - 1)
ny = gv.dimlen(gv.ndim - 2)
levels = gv.dimvalue(gv.ndim - 3)[::-1]
nz = len(levels)
arlf.setlevels(levels)
arlf.set2dvar(avar2d)
for l in levels:
    arlf.set3dvar(avar3d)

#---- Write ARL data file
arlf.setx(gv.dimvalue(gv.ndim - 1))
arlf.sety(gv.dimvalue(gv.ndim - 2))
tNum = inf.timenum()
fhour = 0
for t in range(0, tNum):
    print 'Time index: ' + str(t)
    atime = inf.gettime(t)
    print atime.strftime('%Y-%m-%d %H:00')
    dhead = arlf.getdatahead(inf.proj, 'RSMC', 2, fhour)
    #Pre-write index record without checksum - will be over-write latter
    arlf.writeindexrec(atime, dhead)
    #Checksum list
    ksumlist = []
    # Write 2d variables
    ksums = []
    for avname,gvname in zip(avar2d, gvar2d):
        print avname + ' ' + gvname
        if avname == 'U10M' or avname == 'V10M':
            gdata = inf[gvname][t,0,:,:]
        else:
            gdata = inf[gvname][t,:,:]
        if avname == 'PRSS':
            gdata = gdata * 0.01
        ksum = arlf.writedatarec(atime, 0, avname, fhour, 99, gdata)
        ksums.append(ksum)
    ksumlist.append(ksums)
    # Write 3d variables
    for lidx in range(0, nz):
        ksums = []
        llidx = nz - lidx - 1
        print lidx
        print llidx
        for avname,gvname in zip(avar3d, gvar3d):
            print avname + ' ' + gvname
            gdata = inf[gvname][t,llidx,:,:]
            if avname == 'WWND':
                gdata = gdata * 0.01
            elif avname == 'SPHU':
                gdata = gdata * 1000.
            ksum = arlf.writedatarec(atime, lidx + 1, avname, fhour, 99, gdata)
            ksums.append(ksum)
        ksumlist.append(ksums)
    #Re-write index record with checksum
    arlf.writeindexrec(atime, dhead, ksumlist)
    fhour += 1
arlf.close()
print 'Finished!'

Convert ERA5 grib data to ARL data. The optional field DIFF is recognized as the difference between the original data and the packed data: DIFF=ORIGINAL-PACKED. The effect of this is to increase the precision of variables that have this additional field. When the DIFF field is read by HYSPLIT, the values are added to the data field resulting in values closer to the original data. Currently only DIFW and DIFR (vertical velocity and precipitation) are recognized as valid DIFF fields (https://ready.arl.noaa.gov/hysplitusersguide/S141.htm ).

# Convert ERA5 GRIB data to ARL data

#---- Set data folder
datadir = r'D:/Temp/grib'
#---- Set output data file
outfn = os.path.join(datadir, 'test_era5_grib_diff.arl')
#---- Read a GRIB data file
infn3d = addfile('{}/ERA5_2017.Aug22.3dpl.grib'.format(datadir))
infn2d = addfile('{}/ERA5_2017.Aug22.2dpl.all.grib  '.format(datadir))

print 'GRIB data file has been opened...'
#---- Set output ARL data file
arlf = addfile(outfn, 'c', dtype='arl')

#---- Set variable and level list
#---- Variable names in ERA5 data file
gvar2d = ['Surface_pressure_surface','2_metre_temperature_surface','10_metre_U_wind_component_surface',\
    '10_metre_V_wind_component_surface','Boundary_layer_height_surface','Convective_available_potential_energy_surface',\
    'Instantaneous_eastward_turbulent_surface_stress_surface','Instantaneous_northward_turbulent_surface_stress_surface']

gvar3d = ['Geopotential_isobaric','Temperature_isobaric','Vertical_velocity_isobaric',\
    'U_component_of_wind_isobaric','V_component_of_wind_isobaric','Relative_humidity_isobaric']

#---- Corresponding variable names in ARL data file
avar2d = ['PRSS','T02M','U10M','V10M','PBLH','CAPE','UMOF','VMOF']
avar3d = ['HGTS','TEMP','WWND','UWND','VWND','RELH']
#--- Add DIFF fields - difference between the original data and the packed data
avar3d_diff = list(avar3d)
avar3d_diff.append('DIFW')

#---- Set parameters of ARL data file
gv = infn3d['Geopotential_isobaric']
nx = gv.dimlen(gv.ndim - 1)
ny = gv.dimlen(gv.ndim - 2)
levels = gv.dimvalue(gv.ndim - 3)[::-1]
nz = len(levels)
arlf.setlevels(levels)
arlf.set2dvar(avar2d)
for l in levels:
    arlf.set3dvar(avar3d_diff)
arlf.setx(gv.dimvalue(gv.ndim - 1))
arlf.sety(gv.dimvalue(gv.ndim - 2))

#---- Write ARL data file
tNum = infn3d.timenum()
fhour = 0
for t in range(0, tNum):
    print 'Time index: ' + str(t)
    atime = infn3d.gettime(t)
    print atime.strftime('%Y-%m-%d %H:00')
    dhead = arlf.getdatahead(infn3d.proj, 'RSMC', 2, fhour)
    #Pre-write index record without checksum - will be over-write latter
    arlf.writeindexrec(atime, dhead)
    #Checksum list
    ksumlist = []
    # Write 2d variables
    ksums = []
    for avname,gvname in zip(avar2d, gvar2d):
        gdata = infn2d[gvname][t,:,:]
        if avname == 'PRSS':
            gdata = gdata * 0.01
        ksum = arlf.writedatarec(atime, 0, avname, fhour, 99, gdata)
        ksums.append(ksum)
    ksumlist.append(ksums)
    # Write 3d variables
    for lidx in range(0, nz):
        ksums = []
        llidx = nz - lidx - 1
        print(lidx,llidx)
        for avname,gvname in zip(avar3d, gvar3d):
            gdata = infn3d[gvname][t,llidx,:,:]
            if avname == 'WWND':
                gdata = gdata * 0.01
                difw = arlf.diff_origin_pack(gdata)
            elif avname == 'SPHU':
                gdata = gdata * 1000.
            elif avname == 'HGTS':
                gdata = gdata / 9.80665
            ksum = arlf.writedatarec(atime, lidx + 1, avname, fhour, 99, gdata)
            ksums.append(ksum)
        ksum = arlf.writedatarec(atime, lidx + 1, 'DIFW', fhour, 99, difw)
        ksums.append(ksum)
        ksumlist.append(ksums)
    #Re-write index record with checksum
    arlf.writeindexrec(atime, dhead, ksumlist)
    fhour += 1
arlf.close()
print 'Finished!'

Comparing ERA5 data with converted ARL data. The two data array are not exactly consistant due to the lossy compression algorithm of ARL data format.

ddir = 'D:/Temp/grib'
f_era5_3d = addfile(os.path.join(ddir, 'ERA5_2017.Aug22.3dpl.grib'))
w1 = f_era5_3d['Vertical_velocity_isobaric'][0,-1]    #Pa s**-1
f = addfile(os.path.join(ddir, 'test_era5_grib_diff.arl'))
vname = 'WWND'
w = f[vname][0,0]    #hPa s**-1
difw = f['DIFW'][0,0]
w2 = w + difw
w = w * 100
w2 = w2 * 100

subplot(2,2,1,axestype='map')
geoshow('country')
levs = arange(-1, 1, 0.02)
layer1 = imshow(w1, levs)
colorbar(layer1)
title('ERA5 ({})'.format(vname))

subplot(2,2,2,axestype='map')
geoshow('country')
layer = imshow(w2, levs)
colorbar(layer)
title('ARL + DIFF ({})'.format(vname))

subplot(2,2,3,axestype='map')
geoshow('country')
layer2 = imshow(w1 - w, 20)
colorbar(layer2)
title('ERA5 - ARL ({})'.format(vname))

subplot(2,2,4,axestype='map')
geoshow('country')
layer2 = imshow(w1 - w2, 20)
colorbar(layer2)
title('ERA5 - ARL + DIFF ({})'.format(vname))

savefig('D:/Temp/test/era2arl_diff.png', 650, 370)
../../../_images/era5_arl.png