Source code for pycif.plugins.models.lagrangian.io.native2inputs_adj

import numpy as np
import datetime
import os
import pandas as pd
import xarray as xr
import glob

from logging import info, debug

from ..utils.flexpart_header import read_header, Flexpartheader
from .inputs.inicond_ad import inicond_contribution_ad
from .inputs.fluxes_ad import flux_contribution_ad




def native2inputs_adj(
        self, data2dump, input_type,
        di, df, runsubdir, mode='fwd',
        check_transforms=False, **kwargs):
    """Reads outputs to pycif objects.

       Does nothing for now as we instead read FLEXPART output 
       inside loop over observations in obsoper.py 

    """

    ddi = min(di, df)

    for spec in self.chemistry.acspecies.attributes:
        trid = (input_type, spec)
        if trid not in data2dump:
            continue

        dataobs = self.dataobs[ddi][spec]
        nobs = len(dataobs)
        subdir = ddi.strftime(self.footprint_dir_format)

        # Initialize header
        if self.footprint_type == "STILT":
            fp_header_nest = Flexpartheader()
            fp_header_glob = Flexpartheader()
            fp_header_nest.outheight = [1]

        else:
            # Ref station ID for header
            ref_header = getattr(
                self, "ref_header_ID",
                dataobs.head(1)["metadata"]['station'].values[0].upper())
            
            fp_header_glob = None
            header_nest = "header"
            if self.domain.nested:
                fp_header_glob = read_header(
                    self,
                    os.path.join(
                        self.run_dir_glob,
                        ref_header,
                        subdir, 'header'))
                header_nest = "header_nest"

            if self.force_read_nest and not self.domain.nested:
                header_nest = "header_nest"

            fp_header_nest = read_header(
                self,
                os.path.join(
                    self.run_dir_nest,
                    dataobs["metadata"].head(1)['station']
                    .values[0].upper(),
                    subdir, header_nest))

            fp_header_init = None
            if self.read_background:
                fp_header_init = read_header(
                    self,
                    os.path.join(
                        self.run_dir_bg,
                        dataobs["metadata"].head(1)['station']
                        .values[0].upper(),
                        subdir, 'header'))

            # Nest domain definition
            ix1 = self.domain.ix1
            ix2 = self.domain.ix2
            iy1 = self.domain.iy1
            iy2 = self.domain.iy2

        # Save to datastore for debugging purposes
        obs_ghg = np.nan * np.empty(nobs)
        obs_bkg = np.nan * np.empty(nobs)
        obs_sim = np.nan * np.empty(nobs)
        obs_model = np.nan * np.empty(nobs)
        obs_check = np.nan * np.empty(nobs)
        obs_bkgerr = np.nan * np.empty(nobs)
        obs_err = np.nan * np.empty(nobs)

        info(f"di, df: {di}, {df}, {datetime.datetime.now()}")

        # Apply sensitivity to background if required
        if self.read_background and input_type == "inicond":
            ini_dates = pd.DatetimeIndex(
                self.datainicond[("inicond", spec)][ddi]["spec"].time.values
            ).to_pydatetime()
            inicond_sensit = inicond_contribution_ad(
                self, mode, dataobs, fp_header_init, spec, ddi)
            data2dump[trid]["data"][ddi]["adj_out"] = xr.DataArray(
                inicond_sensit,
                coords={"time": ini_dates},
                dims=("time", "lev", "lat", "lon"))

        # Apply fluxes contribution
        if self.read_surface_sensitivity and input_type == "flux":
            flux = {"spec": self.dataflx[("flux", spec)][ddi]["spec"]}
            dataflx = flux["spec"]
            flx_dates = pd.DatetimeIndex(
                dataflx.time.values).to_pydatetime()
            flx_sensit = flux_contribution_ad(
                self, mode, dataobs, fp_header_nest, fp_header_glob, spec, ddi)

            data2dump[trid]["data"][ddi]["adj_out"] = xr.DataArray(
                flx_sensit,
                coords={"time": np.append(
                    self.input_dates[ddi][0]
                    - pd.to_timedelta(self.backward_trajdays),
                    self.input_dates[ddi])},
                dims=("time", "lev", "lat", "lon"))

        continue

    return data2dump