Source code for pycif.plugins.models.lagrangian.utils.flexpart_header

""" Class for reading FLEXPART binary headers

"""

import datetime
import struct
import numpy as np
from scipy.io import FortranFile




def read_header(self, filename):
    """Reading FLEXPART header file"""

    flexpart_header = Flexpartheader()
    
    with FortranFile(filename, "r") as f:
        # Deal with difference between FLEXPART 10 and previous versions
        if getattr(self, "flexversion", 10) >= 10:
            len_ver = 256
        else:
            len_ver = 13
        
        tmp = f.read_record(f'<i4,<i4,S{len_ver}')[0]
        flexpart_header.ibdate = tmp[0]
        flexpart_header.ibtime = tmp[1]
        
        loutstep, loutaver, loutsample = f.read_ints('i4')

        tmp = f.read_record('<f4', '<f4', '<i4', '<i4', '<f4', '<f4')
        flexpart_header.bndx = tmp[0][0]
        flexpart_header.bndy = tmp[1][0]
        flexpart_header.numx = tmp[2][0]
        flexpart_header.numy = tmp[3][0]
        flexpart_header.delx = tmp[4][0]
        flexpart_header.dely = tmp[5][0]

        tmp = f.read_reals('f4')
        nzgrid = struct.unpack('i', struct.pack('f', tmp[0]))[0]
        flexpart_header.outheight = np.zeros(50)
        flexpart_header.outheight[:nzgrid] = tmp[1: nzgrid + 1]

        jjjjmmdd, ihmmss = f.read_ints('i4')
        nspec, nzgrid = f.read_ints('i4')
        nspec /= 3
        nspec = int(nspec)
        for k in np.arange(nspec):
            _ = f.read_ints('<i2')
            _ = f.read_ints('<i2')
            _ = f.read_ints('<i2')

        flexpart_header.numpoint = f.read_ints()[0]
        ireleasestart = np.zeros(flexpart_header.numpoint)
        ireleaseend = np.zeros(flexpart_header.numpoint)
        xpoint = np.zeros(flexpart_header.numpoint)
        ypoint = np.zeros(flexpart_header.numpoint)
        zpoint1 = np.zeros(flexpart_header.numpoint)
        zpoint2 = np.zeros(flexpart_header.numpoint)
        npart = np.zeros(flexpart_header.numpoint)
        nkind = np.zeros(flexpart_header.numpoint)
        xmass = np.zeros((flexpart_header.numpoint, nspec))
        compoint = np.empty(flexpart_header.numpoint, dtype="|S45")

        for i in np.arange(flexpart_header.numpoint):
            line = f.read_record(np.dtype("a10"))
            ireleasestart[i] = struct.unpack('i', line[0][:4])[0]
            xpoint[i], ypoint[i], xp2, yp2, zpoint1[i], zpoint2[i] = f.read_reals("f4")
            npart[i], nkind[i] = f.read_ints("i4")
            compoint[i] = f.read_record(np.dtype("a45"))[0]
            for j in np.arange(nspec):
                _ = f.read_ints()
                _ = f.read_ints()
                xmass[i, j] = f.read_reals("f4")

        imethod = f.read_ints()

        tmp = f.read_ints('i4')
        nageclass = tmp[0]
        lage = tmp[1:1 + nageclass]

        flexpart_header.releases = np.zeros(1000)
        flexpart_header.releases[:flexpart_header.numpoint] = ireleasestart

        flexpart_header.ageclass = lage[0]
        flexpart_header.trajdays = flexpart_header.ageclass / 3600. / 24.
        flexpart_header.ndgrid = np.abs(24 * 3600 / loutstep)

        if flexpart_header.numx * flexpart_header.delx == 360:
            flexpart_header.xshift = \
                int((flexpart_header.bndx + 180) / flexpart_header.delx)
            flexpart_header.bndx = -180

        else:
            flexpart_header.xshift = 0

    return flexpart_header





class Flexpartheader(object):
    """ Class for reading FLEXPART binary headers

    """

    def __init__(self):
        """Initialize a Flexpartheader object

        """

        self.numpoint = np.nan
        self.ibdate = np.nan
        self.ibtime = np.nan
        self.releases = np.nan
        self.numx = np.nan
        self.numy = np.nan
        self.bndx = np.nan
        self.bndy = np.nan
        self.delx = np.nan
        self.dely = np.nan
        self.xshift = np.nan
        self.ageclass = np.nan
        self.trajdays = np.nan
        self.ndgrid = np.nan
        self.outheight = np.nan
        self.maxngrid = np.nan