将MERRA2数据转为ARL数据报错

LLLLLYH

代码结合MeteoInfo转wrfout为arl和hysplit中mer2arl.f，能转为arl数据，但是用addfile打开的时候出现了报错，在HYSPLIT中也不能正常使用

import os

datadir = '/Users/luoyh23/hysplit/data2arl/mer2arl'

prefix = 'MERRA2_100.tavg'
suffix = '.nc4'

yyyy = 1980
mon = 6
dd = 22
for day in range(dd,dd+1):

    mm='%02d' % mon
    daystr='%02d' % day
    cdate = str(yyyy)+mm+daystr
    outfn = os.path.join(datadir, 'mer2arl'+cdate+'.arl')

    infn1 = os.path.join(datadir, prefix+'3_3d_asm_Nv.'+cdate+suffix)
    infn2 = os.path.join(datadir, prefix+'1_2d_flx_Nx.'+cdate+suffix)
    infn3 = os.path.join(datadir, prefix+'1_2d_slv_Nx.'+cdate+suffix)
    infn4 = os.path.join(datadir, prefix+'1_2d_int_Nx.'+cdate+suffix)

    inf1 = addfile(infn1)
    inf2 = addfile(infn2)
    inf3 = addfile(infn3)
    inf4 = addfile(infn4)

    arlf = addfile(outfn, 'c', dtype='arl')

    wvar2d = ['PHIS','PS','PRECTOTCORR','PBLH','USTAR','Z0M','SWNETSRF','HFLUX','EFLUX','T2M','U10M','V10M']
    avar2d = ['SHGT','PRSS','PRT3','PBLH','USTR','RGHS','DSWF','SHTF','LHTF','T02M','U10M','V10M']

    wvar3d = ['PL','U','V','OMEGA','T','QV']
    avar3d = ['PRES','UWND','VWND','WWND','TEMP','SPHU']

    wv = inf1['PL']
    nx = wv.dimlen(wv.ndim - 1)
    ny = wv.dimlen(wv.ndim - 2)

    levels = wv.dimvalue(wv.ndim - 3)
    nz = len(levels)
    arlf.setlevels(levels)
    arlf.set2dvar(avar2d)
    for l in levels:
        arlf.set3dvar(avar3d)

    arlf.setx(wv.dimvalue(wv.ndim - 1))
    arlf.sety(wv.dimvalue(wv.ndim - 2))

    tNum1 = inf1.timenum()
    tNum2 = inf2.timenum()
    tNum3 = inf3.timenum()
    tNum4 = inf4.timenum()

    fhour = 3

    #loop through each time step
    for t in range(0, tNum1):
        #print 'Time index: ' + str(t)
        ksumlist = []
        ksums = []
        atime1 = inf1.gettime(t)
        dhead = arlf.getdatahead(inf1.proj,'MER2', 2, fhour)
        arlf.writeindexrec(atime1, dhead)

        #hourly data for flux vatiables
        for k in range(0,tNum2):
            atime2 = inf2.gettime(k)
            if atime1==atime2 or atime2 > atime1:
                iflx = k
                break

        #hourly data for 2d variables
        for k in range(0,tNum3):
            atime3 = inf3.gettime(k)
            if atime1==atime3 or atime3 > atime1:
                islv=k
                iint=k
                break

        #extract 2d variables
        for avname,wvname in zip(avar2d, wvar2d):
            if wvname in ['PS','T2M','U10M','V10M']:
                gdata = inf3[wvname][islv,:,:]
            elif wvname in ['SWNETSRF','HFLUX']:
                gdata = inf4[wvname][iint,:,:]
            elif wvname in ['PRECTOTCORR','PBLH','USTAR','Z0M','EFLUX']:
                gdata = inf2[wvname][iflx,:,:]
            elif wvname == 'PHIS':
                gdata = inf1[wvname][t,:,:]

            ksum = arlf.writedatarec(atime1, 0, avname, fhour, 99, gdata)
            ksums.append(ksum)

        ksumlist.append(ksums)

        #write 3d variables
        for lidx in range(0, nz):
            ksums = []

            for avname,wvname in zip(avar3d, wvar3d):
                gdata = inf1[wvname][t,lidx,:,:]
                ksum = arlf.writedatarec(atime1, lidx + 1, avname, fhour, 99, gdata)
                ksums.append(ksum)

            ksumlist.append(ksums)

        arlf.writeindexrec(atime1, dhead, ksumlist)
        fhour += 0

    arlf.close()
    print 'Finished!'

f = addfile(outfn)

最后打开该arl文件的时候出现了以下的报错

LLLLLYH

已解决，原因是垂直层颠倒了，手动颠倒后就能用了，顺便调整了各个变量的单位转换，能运行的代码如下，希望能给以后转换MERRA2的朋友参考

import os

datadir = '/Users/luoyh23/hysplit/data2arl/mer2arl'

prefix = 'MERRA2_100.tavg'
suffix = '.nc4'

yyyy = 1980
mon = 6
dd = 22

#转换事件发生前10天的数据
for day in range(dd-10,dd+1):

    mm='%02d' % mon
    daystr='%02d' % day
    cdate = str(yyyy)+mm+daystr
    outfn = os.path.join(datadir, 'mer2arl'+cdate+'.arl')

    infn1 = os.path.join(datadir, prefix+'3_3d_asm_Nv.'+cdate+suffix)
    infn2 = os.path.join(datadir, prefix+'1_2d_flx_Nx.'+cdate+suffix)
    infn3 = os.path.join(datadir, prefix+'1_2d_slv_Nx.'+cdate+suffix)
    infn4 = os.path.join(datadir, prefix+'1_2d_int_Nx.'+cdate+suffix)

    inf1 = addfile(infn1)
    inf2 = addfile(infn2)
    inf3 = addfile(infn3)
    inf4 = addfile(infn4)

    arlf = addfile(outfn, 'c', dtype='arl')

    conv2d = [0.102, 0.01, 0.06, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
    wvar2d = ['PHIS','PS','PRECTOTCORR','PBLH','USTAR','Z0M','SWNETSRF','HFLUX','EFLUX','T2M','U10M','V10M']
    avar2d = ['SHGT','PRSS','PRT3','PBLH','USTR','RGHS','DSWF','SHTF','LHTF','T02M','U10M','V10M']

    conv3d = [0.01, 1.0, 1.0, 1.0, 0.01, 1.0]
    wvar3d = ['PL','U','V','OMEGA','T','QV']
    avar3d = ['PRES','UWND','VWND','WWND','TEMP','SPHU']
   
    wv = inf1['PL']
    nx = wv.dimlen(wv.ndim - 1)
    ny = wv.dimlen(wv.ndim - 2)

    levels = wv.dimvalue(wv.ndim - 3)[::-1]
    nz = len(levels)
    arlf.setlevels(levels)
    arlf.set2dvar(avar2d)
    for l in levels:
        arlf.set3dvar(avar3d)

    arlf.setx(wv.dimvalue(wv.ndim - 1))
    arlf.sety(wv.dimvalue(wv.ndim - 2))

    tNum1 = inf1.timenum()
    tNum2 = inf2.timenum()
    tNum3 = inf3.timenum()
    tNum4 = inf4.timenum()

    fhour = 3

    #loop through each time step
    for t in range(0, tNum1):
        #print 'Time index: ' + str(t)
        datamin = 90 + t*180

        atime1 = inf1.gettime(t)
        dhead = arlf.getdatahead(inf1.proj,'MER2', 2, 0)
        arlf.writeindexrec(atime1, dhead)
   
        #hourly data for flux vatiables
        for k in range(0,tNum2):
            atime2 = inf2.gettime(k)
            if atime1==atime2 or atime2 > atime1:
                iflx = k
                break

        #hourly data for 2d variables
        for k in range(0,tNum3):
            atime3 = inf3.gettime(k)
            if atime1==atime3 or atime3 > atime1:
                islv=k
                iint=k
                break

        ksumlist = []
        ksums = []
        #extract 2d variables
        for coeff,avname,wvname in zip(conv2d,avar2d, wvar2d):
            if wvname in ['PS','T2M','U10M','V10M']:
                gdata = coeff*inf3[wvname][islv,:,:]
            elif wvname in ['SWNETSRF','HFLUX']:
                gdata = coeff*inf4[wvname][iint,:,:]
            elif wvname in ['PBLH','USTAR','Z0M','EFLUX']:
                gdata = coeff*inf2[wvname][iflx,:,:]
            elif wvname == 'PHIS':
                gdata = coeff*inf1[wvname][t,:,:]

            #special for PRT3
            if wvname == 'PRECTOTCORR':
                gdata = coeff*(inf2[wvname][iflx,:,:]+inf2[wvname][iflx-1,:,:]+inf2[wvname][iflx+1,:,:])/3.
                 
            ksum = arlf.writedatarec(atime1, 0, avname, fhour, 99, gdata)
            ksums.append(ksum)

        ksumlist.append(ksums)

        #write 3d variables
        for lidx0 in range(0, nz):
            ksums = []
            lidx = nz - lidx0 -1
   
            for coeff,avname,wvname in zip(conv3d,avar3d, wvar3d):
                gdata = coeff*inf1[wvname][t,lidx,:,:]
                ksum = arlf.writedatarec(atime1, lidx + 1, avname, fhour, 99, gdata)
                ksums.append(ksum)

                ksumlist.append(ksums)

        arlf.writeindexrec(atime1, dhead, ksumlist)
        fhour += 1
    arlf.close()
    print 'Finished!'
   

LLLLLYH

链接: https://pan.baidu.com/s/1aN46Q0NJFS8_P_wk7rAvJA 提取码: vkda
该链接中为以上代码所用到的19800622MERRA2数据

MeteoInfo

LLLLLYH 发表于 2022-2-26 16:17
链接: https://pan.baidu.com/s/1aN46Q0NJFS8_P_wk7rAvJA 提取码: vkda
该链接中为以上代码所用到的198006 ...

把已经生成的转换文件删除再试试

LLLLLYH

MeteoInfo 发表于 2022-2-26 17:54
把已经生成的转换文件删除再试试

感谢王老师的回复，每次都是删除后再运行，还是有这个问题

LLLLLYH

MeteoInfo 发表于 2022-2-26 17:54
把已经生成的转换文件删除再试试

王老师，这里还想请教一下您，我去查阅了源代码，里面显示的主要是getdatahead的问题，并且我看dimdatafile.py里的这个函数有五个参数，我就在最后加上了一个参数，结果没有改变，然后我输出了每次的arlf.getdatahead，发现都是空的（<bound method DimDataFile.getdatahead of None>），我有点不太明白为什么可以一开始就可以arlf.getdatahead而不需要先写入，但是我不太确定是不是这里的问题。

MeteoInfo

LLLLLYH 发表于 2022-2-27 20:59
已解决，原因是垂直层颠倒了，手动颠倒后就能用了，顺便调整了各个变量的单位转换，能运行的代码如下，希望 ...

赞一个

小徐同学

    楼主您好，我按照附件里的程序对2006年1月1日的数据进行转换，发现转换后的数据跟原场不匹配，具体表现为二维数据（U10m、V10m等）不同格点数值最后几位为相同值，三维数据-9999值的点数值发生变化并且其他格点为同一数值。
    请问楼主知道怎么解决吗？{:eb303:}{:eb303:}（程序数据和转换后的数据放在下面的网盘链接里了）
链接：https://pan.baidu.com/s/1KDcFDH5UAtzyJlS1eGkxvA?
提取码：1234

LLLLLYH

小徐同学 发表于 2022-5-12 21:15
楼主您好，我按照附件里的程序对2006年1月1日的数据进行转换，发现转换后的数据跟原场不匹配，具体表现 ...

不好意思，我也遇到了这个问题，当时没发现，后来追踪的时候发现了，但是因为时间比较紧急就换成ERA5数据来做了，目前主要觉得可能是经度需要从-180-180转到0-360，你可以试试这样能不能有效，还有关于地形的部分，MERRA2里面涉及地形的格点数据写的NAN，但是ARL识别不了，所以我才改成9999，我也不太确定ARL是怎么识别缺测的。如果你有结果了还希望可以在这里分享一下