Source code for pycif.plugins.obsparsers.tropomi_official.parse

from logging import info, debug
import xarray as xr
import pandas as pd
import numpy as np
import datetime
import netCDF4 as nc
from ....utils.datastores.empty import init_empty




def do_parse(
    self,
    obs_file,
    **kwargs
):
    """Parse function for a file from template observations

    Args:
        obs_file (str) :
            Path to input file

    Returns:
        pandas.DataFrame :
            Dataframe from input file df[parameter][station]

    """

    info("Parsing official TROPOMI data "
         "from observation file: {}".format(obs_file))
    
    ddi, ddf = self.datei, self.datef

    product_data = xr.open_dataset(obs_file, group="PRODUCT")
    obs_dates = product_data["time_utc"].values.astype(
        np.datetime64).astype(datetime.datetime)
    self.data_date_max = obs_dates.max()

    if ddi > obs_dates.max() or ddf < obs_dates.min():
        debug(f"File {obs_file} with date {obs_dates.min()} "
              f"is outside simulation time period {ddi} / {ddf}")
        return init_empty()

    # Filter out data outside domain
    longitude = product_data.variables['longitude'].values
    latitude = product_data.variables['latitude'].values
    mask_domain = np.ones_like(longitude, dtype=bool)
    if hasattr(self, "crop_domain"):
        crop_domain = self.crop_domain
        mask_domain = ((longitude < crop_domain.xmax) &
                       (longitude > crop_domain.xmin) &
                       (latitude < crop_domain.ymax) &
                       (latitude > crop_domain.ymin))
        if mask_domain.sum() == 0:
            return init_empty()
    
    # Check number of levels
    nlevread = product_data.dims["layer"]
    if nlevread != self.nlayers:
        raise Exception(
            f"The file {obs_file} does not have the expected number of "
            f"averaging kernel layers: {nlevread} instead of "
            f"{self.nlayers}. \n"
            f"Please check your files or change the parameter 'nlayers' "
            f"in the Yaml (expert users only).")
    
    # Valid data
    scanline = product_data.dims["scanline"]
    ground_pixel = product_data.dims["ground_pixel"]
    time = product_data.dims["time"]
    obs = product_data.variables['methane_mixing_ratio'].values
    mask_valid = product_data.variables['qa_value'].values >= self.qa_value
    mask_valid = mask_valid & ~np.isnan(obs)

    if (mask_domain & mask_valid).sum() == 0:
        return init_empty()
    
    # Filter out data not in simulation window
    mask_time = (obs_dates >= ddi) & (obs_dates <= ddf)
    mask_all = mask_time[..., np.newaxis] & mask_valid & mask_domain
    if mask_all.sum() == 0:
        return init_empty()

    # Now process actual data
    debug(f"Data found in file {obs_file}")
    list_basic_cols = [
        'date', 'duration', 'station', 'network', 'parameter',
        'lon', 'lat', 'obs', 'obserror']

    subdata = pd.DataFrame(columns=list_basic_cols)
    subdata['date'] = np.repeat(
        obs_dates[..., np.newaxis], ground_pixel, axis=2)[mask_all]
    subdata = subdata.assign(duration=self.obs_duration)
    subdata['lon'] = longitude[mask_all]
    subdata['lat'] = latitude[mask_all]
    subdata = subdata.assign(station=self.stationID)
    subdata = subdata.assign(network=self.networkID)
    subdata = subdata.assign(parameter=self.parameter)
    
    subdata["obs"] = obs[mask_all].astype(float)
    subdata["obserror"] = product_data.variables[
        'methane_mixing_ratio_precision'].values[mask_all].astype(float)
    
    # Load aks and obs
    detailed_results_grp = xr.open_dataset(
        obs_file, group='PRODUCT/SUPPORT_DATA/DETAILED_RESULTS')
    ak = np.float64(
        detailed_results_grp.variables['column_averaging_kernel']
        .values[mask_all, :]
    )

    # Pressure levels
    # particular case of TROPOM CH4: create pressure grid from psurf and delta
    input_data_grp = xr.open_dataset(
        obs_file, group='PRODUCT/SUPPORT_DATA/INPUT_DATA')
    psurf = input_data_grp.variables['surface_pressure'].values[mask_all]
    deltap = input_data_grp.variables['pressure_interval'].values[mask_all]
    pavg0 = np.array([psurf - deltap * lev for lev in range(self.nlayers + 1)])
    pavg = np.float64(np.transpose(pavg0))
    
    # prior profile in molecules_percm2
    qa0 = np.float64(
        input_data_grp.variables['methane_profile_apriori'].values[mask_all, :]
    )
    conv = input_data_grp.variables[
        'methane_profile_apriori'
    ].attrs['multiplication_factor_to_convert_to_molecules_percm2']
    qa0 *= conv
    
    # information specific to TROPOMI
    # dry air subcolumns - required if not provided by the CTM (to update if O
    dvair = np.float64(
        input_data_grp.variables['dry_air_subcolumns'].values[mask_all, :]
    )
    conv = input_data_grp.variables[
        'dry_air_subcolumns'
    ].attrs['multiplication_factor_to_convert_to_molecules_percm2']
    dvair *= conv

    ds = xr.Dataset({'qa0': (['index', 'level'], qa0),
                     'ak': (['index', 'level'], ak),
                     'pavg0': (['index', 'level_pressure'], pavg),
                     'dryair': (['index', 'level'], dvair)},
                    coords={'index': subdata.index,
                            'level': range(self.nlayers),
                            'level_pressure': range(self.nlayers + 1)})

    subdata = subdata.to_xarray()
    subdata = xr.merge([ds, subdata])

    debug(f"Observations after filter: {len(subdata.index)}")

    return subdata