"""Interpolation methods in the time domain."""
from pathlib import Path
import numpy as np
import xarray as xr
from scipy.interpolate import RegularGridInterpolator
from .read_netcdf import get_evolution_type
from .utils import convert_mass_or_annual_rate
# from scipy.interpolate import interp1d
# CONSTANTS
# Translation table from evolution keys to inventory keys
KEY_TABLE = {
"fuel": "fuel",
"EI_CO2": "CO2",
"EI_H2O": "H2O",
"EI_NOx": "NOx",
"dis_per_fuel": "distance",
}
[docs]
def interpolate(
config: dict, years: np.ndarray, val_dict: dict
) -> tuple[np.ndarray, dict]:
"""Interpolate values from discrete years to time_range set in config.
Args:
config (dict): Configuration dictionary from config
years (numpy.ndarray): Numpy array with discrete years
val_dict (dict): Dictionary with time series numpy arrays, keys are species
Returns:
tuple[numpy.ndarray, dict]: Time range over which interpolation takes place,
and interpolated values, keys are the same as in val_dict
"""
# TODO extend this function with several interpolation methods defined in config
time_range, interp_dict = interp_linear(config, years, val_dict)
return time_range, interp_dict
[docs]
def interp_linear(
config: dict,
years: np.ndarray,
val_dict: dict,
bounds_error: bool = True,
fill_value: float | None = np.nan,
) -> tuple[np.ndarray, dict]:
"""Interpolate linearly values from discrete years to time_range set in config.
Args:
config (dict): Configuration dictionary from config
years (numpy.ndarray): Numpy array with discrete years
val_dict (dict): Dictionary with time series numpy arrays, keys are species
bounds_error (bool, optional): See documentation of scipy.interpolate.interp1d
fill_value (float or None, optional):
See documentation of scipy.interpolate.RegularGridInterpolator
Raises:
IndexError: if too few discrete years or inappropriate options set
Returns:
numpy.ndarray, dict: Time range over which interpolation takes place,
and interpolated values, keys are the same as in val_dict
"""
time_config = config["time"]["range"]
time_range = np.arange(time_config[0], time_config[1], time_config[2], dtype=int)
interp_dict = {}
for key, values in val_dict.items():
# Under certain circumstances, interpolation works also for 1 given inventory
if len(years) > 1 or (len(years) == 1 and isinstance(fill_value, (float, int))):
# interp_func = interp1d(
# years,
# values,
# kind="linear",
# bounds_error=bounds_error,
# fill_value=fill_value,
# )
# interp_dict[key] = interp_func(time_range)
interp_func = RegularGridInterpolator(
[years],
values,
method="linear",
bounds_error=bounds_error,
fill_value=fill_value,
)
new_grid = np.ix_(time_range)
interp_dict[key] = interp_func(new_grid)
# only one discrete year given,
# repeat value array times number of years in time_range
elif len(years) == 1 and fill_value is None:
# interp_dict[key] = np.ones(len(time_range)) * values
interp_dict[key] = np.squeeze(np.array([values] * len(time_range)))
else:
raise IndexError(
"Interpolation not possible! Too few discrete values or "
"inappropriate options set for function interp_linear"
)
return time_range, interp_dict
[docs]
def adjust_inventories(config: dict, inv_dict: dict) -> dict:
"""Determine evolution_type and adjust inventories depending on type.
Data variables to be adjusted: fuel, species emissions and distance.
Args:
config (dict): Configuration dictionary from config
inv_dict (dict): Dictionary of xarray Datasets, keys are years of
input inventories
Raises:
ValueError: if no valid evolution_type in evolution file
Returns:
dict: normalized/scaled/unmodified dictionary of xarray Datasets,
keys are years of input inventories
"""
evolution_type = get_evolution_type(config)
if evolution_type == "norm":
out_dict = norm_inventories(config, inv_dict)
elif evolution_type == "scaling":
out_dict = scale_inventories(config, inv_dict)
elif evolution_type is False:
out_dict = inv_dict
else:
raise ValueError(
"Invalid evolution_type "
"evolution_type can be either 'scaling', 'norm' or False"
)
return out_dict
[docs]
def apply_evolution(
config: dict,
val_dict: dict,
inv_dict: dict,
inventories_adjusted: bool = False,
) -> tuple[np.ndarray, dict]:
"""Determine evolution_type, apply normalization, scaling or no evolution.
Args:
config (dict): Configuration dictionary from config
val_dict (dict): Dictionary with time series numpy arrays, keys are species
inv_dict (dict): Dictionary of xarray Datasets, keys are years of
input inventories
inventories_adjusted (bool): True if inventories have been adjusted beforehand
through adjust_inventories(config, inv_dict), defaults to False
Raises:
ValueError: if no valid evolution_type in evolution file
Returns:
numpy.ndarray, dict: time_range and normalized/scaled/unmodified
dictionary, ``{spec: numpy.ndarray}``
"""
evolution_type = get_evolution_type(config)
if evolution_type == "scaling":
time_range, out_dict = apply_scaling(
config, val_dict, inv_dict, inventories_adjusted
)
elif evolution_type == "norm":
time_range, out_dict = apply_norm(config, val_dict, inv_dict)
elif evolution_type is False:
time_range, out_dict = apply_no_evolution(config, val_dict, inv_dict)
else:
raise ValueError(
"Invalid evolution_type "
"evolution_type can be either 'scaling', 'norm' or False"
)
return time_range, out_dict
[docs]
def apply_scaling(
config: dict, val_dict: dict, inv_dict: dict, inventories_adjusted: bool
) -> tuple[np.ndarray, dict]:
"""Apply scaling on dictionary of time series.
Scaling factors are from the evolution file; time series and scaling
factors are interpolated on time_range before multiplication.
TODO: implement scaling for individual species.
Args:
config (dict): Configuration dictionary from config
val_dict (dict): Dictionary with time series numpy arrays, keys are species
inv_dict (dict): Dictionary of xarray Datasets, keys are years of
input inventories
inventories_adjusted (bool): True if inventories have been adjusted beforehand
through adjust_inventories(config, inv_dict), False otherwise
Returns:
numpy.ndarray, dict: time_range and scaled dictionary
"""
# Get inventory years
inv_years = np.array(list(inv_dict.keys()))
# Interpolate scaling factors linearly on time_range
time_range, evo_interp_dict = interp_evolution(config)
# If inventories have been adjusted beforehand, normalize scaling
# factors to inv_years
if inventories_adjusted:
# Filter evo_interp_dict to inv_years
evo_filtered_dict = filter_to_inv_years(inv_years, time_range, evo_interp_dict)
# Interpolate linearly evo_filtered_dict
_time_range, evo_filtered_interp_dict = interp_linear(
config, inv_years, evo_filtered_dict
)
# Normalize evolution scaling factors to evo_filtered_interp_dict
evo_norm_scaling_arr = np.divide(
evo_interp_dict["scaling"], evo_filtered_interp_dict["scaling"]
)
evo_interp_dict = {"scaling": evo_norm_scaling_arr}
# Interpolate time series data linearly on time_range
_time_range, interp_dict = interp_linear(config, inv_years, val_dict)
# Multiply time series data by scaling factors
out_dict = {}
for spec, series_arr in interp_dict.items():
# adapted for multiplication of multi-dimensional arrays
shape_tp = np.shape(np.transpose(series_arr))
scaling = np.transpose(np.resize(evo_interp_dict["scaling"], shape_tp))
out_dict[spec] = np.multiply(scaling, series_arr)
return time_range, out_dict
[docs]
def apply_norm(config: dict, val_dict: dict, inv_dict: dict) -> tuple[np.ndarray, dict]:
"""Apply normalization on time series.
Gets data from evolution file and inventories, interpolates time series
and evolution data over time_range, then calculates normalization
factors and applies them (multiplication).
Args:
config (dict): Configuration dictionary from config
val_dict (dict): Dictionary with time series numpy array, keys are species
inv_dict (dict): Dictionary of xarray Datasets, keys are years of
input inventories
Returns:
numpy.ndarray, dict: time_range and normalized dictionary
"""
out_dict = {}
# Interpolate evolution over time_range
_time_range, evo_interp_dict = interp_evolution(config)
# Calculate fuel sums and emission indices from inventories
# {"fuel": np.ndarray, "EI_CO2": np.ndarray, ..}
# TODO This step might be redundant if normalization is applied beforehand
# on input inventories. In this case, ei_inv_dict valures are exactly
# the values in evo_interp_dict for inventory years.
_inv_years, _inv_sum_dict, ei_inv_dict = calc_inv_quantities(config, inv_dict)
# Filter emission indices dictionary to those species specified in time evolution
ei_inv_dict = filter_dict_to_evo_keys(config, ei_inv_dict)
# Interpolate emission indices from inventories over time_range
# {"fuel": np.ndarray, "EI_CO2": np.ndarray, ..}
inv_years = np.array(list(inv_dict.keys()))
# TODO Check bounds, extrapolate
time_range, ei_inv_interp_dict = interp_linear(
config,
inv_years,
ei_inv_dict,
bounds_error=False,
fill_value=None,
)
# Calculate normalization factors, over time_range
norm_dict = calc_norm(evo_interp_dict, ei_inv_interp_dict)
norm_fuel = norm_dict["fuel"]
# Interpolate time series data on time_range
_time_range, interp_dict = interp_linear(
config,
inv_years,
val_dict,
bounds_error=False,
fill_value=None,
)
# Apply normalization factors to interp_dict
for spec, series_arr in interp_dict.items():
# adapted for multiplication of multi-dimensional arrays
shape_tp = np.shape(np.transpose(series_arr))
if spec in norm_dict:
norm = np.transpose(np.resize(norm_dict[spec], shape_tp))
normalized = np.multiply(norm, series_arr)
else:
norm = np.transpose(np.resize(norm_fuel, shape_tp))
normalized = np.multiply(norm, series_arr)
out_dict[spec] = normalized
return time_range, out_dict
[docs]
def apply_no_evolution(
config: dict, val_dict: dict, inv_dict: dict
) -> tuple[np.ndarray, dict]:
"""Apply interpolation only on time series data over time_range.
Args:
config (dict): Configuration dictionary from config
val_dict (dict): Dictionary with time series numpy arrays, keys are species
inv_dict (dict): Dictionary of xarray Datasets, keys are years of
input inventories
Returns:
numpy.ndarray, dict: time_range and dictionary of time series data
"""
time_config = config["time"]["range"]
time_range = np.arange(time_config[0], time_config[1], time_config[2], dtype=int)
# Get inventory years
inv_years = np.array(list(inv_dict.keys()))
# Interpolate time series data on time_range
_time_range, interp_dict = interp_linear(
config, inv_years, val_dict, bounds_error=False, fill_value=0.0
)
return time_range, interp_dict
[docs]
def interp_evolution(config: dict) -> tuple[np.ndarray, dict]:
"""Interpolate values in evolution file over time_range.
Args:
config (dict): Configuration dictionary from config
Returns:
array, dict: time_range and Dictionary with interpolated values
"""
time_dir = config["time"]["dir"]
file_name = config["time"]["file"]
file_path = Path(time_dir) / file_name
evolution = xr.load_dataset(file_path)
evo_years = evolution.time.values
evo_dict: dict[str, float | np.ndarray] = {}
for raw_key, var in evolution.items():
key = str(raw_key)
arr = var.values
if key == "fuel":
evo_dict[key] = convert_mass_or_annual_rate(arr, var.attrs["units"], "kg")
else:
evo_dict[key] = arr
time_range, evo_interp_dict = interp_linear(config, evo_years, evo_dict)
return time_range, evo_interp_dict
[docs]
def calc_inv_quantities(config: dict, inv_dict: dict) -> tuple[np.ndarray, dict, dict]:
"""Calculate inventory quantities: fuel sums, emission sums and emission indices.
Sums and emission indices are only calculated from species included in config.
Args:
config (dict): Configuration dictionary from config
inv_dict (dict): Dictionary of xarray Datasets, keys are years of inventories
Returns:
numpy.ndarray, dict, dict: Array of inventory years,
Dictionary of arrays of summed inventory emissions, keys are species
Dictionary of arrays of fuel sums and inventory emission indices,
keys are data variable names of evolution file
"""
#
# Translation table from evolution keys to inventory keys
evo_inv_table = KEY_TABLE
# Invert translation table: inventory keys to evolution keys
inv_evo_table = {v: k for k, v in evo_inv_table.items()}
# Initialize output array and dictionaries
# Inventory sums, keys are species, values are arrays over years
inv_sum_dict: dict = {}
# Emission indices, keys are species, values are arrays over years
ei_inv_dict = {}
# Initialize lists in inv_sum_dict from inventory species array defined in config
# TODO If inventories have missing species for different years, this is an issue!
# Check for missing species in read_netcdf.py -> open_inventories()
spec_lst = ["fuel", *config["species"]["inv"]]
for spec in spec_lst:
inv_sum_dict[spec] = []
# Iterate over inventories
for inv in inv_dict.values():
for spec, data_arr in inv.items():
# Skip coordinates
if spec in ["lon", "lat", "plev", "ac"] or spec not in spec_lst:
pass
else:
spec_sum = data_arr.sum().values.item()
inv_sum_dict[spec].append(spec_sum)
# Convert lists into numpy arrays
inv_years = np.array(list(inv_dict.keys()))
for spec, sum_arr in inv_sum_dict.items():
inv_sum_dict[spec] = np.array(sum_arr)
fuel_sum_arr = inv_sum_dict["fuel"]
# Add array of inventory fuel sums to emission index dictionary
ei_inv_dict["fuel"] = fuel_sum_arr
# Calculate emission indices for each species
for spec, sum_arr in inv_sum_dict.items():
if spec != "fuel":
# Calculate emission index for spec (array over inventory years)
ei_arr = np.divide(sum_arr, fuel_sum_arr)
# Get right evolution key from translation table
evo_key = inv_evo_table[spec]
# Add array of emission indices for spec to emission index dictionary
ei_inv_dict[evo_key] = ei_arr
return inv_years, inv_sum_dict, ei_inv_dict
[docs]
def filter_dict_to_evo_keys(config: dict, inp_dict: dict) -> dict:
"""Filter input dictionary to items matching with keys from time evolution file.
This function loads the time evolution file as an xarray dataset, and
iterates over its keys. An item from the input dictionary is transferred
to the output dictionary if the corresponding input key matches with a
key from the time evolution.
Args:
config (dict): Configuration dictionary
inp_dict (dict): Input dictionary
Returns:
dict: Dictionary with filtered items
Raises:
KeyError: If no matches are found between both sets of keys.
"""
time_dir = config["time"]["dir"]
file_name = config["time"]["file"]
file_path = Path(time_dir) / file_name
# Output dictionary
out_dict = {}
# Load the time evolution file as an xarray dataset
evolution = xr.load_dataset(file_path)
# Iterate over keys in evolution
for evo_key in evolution:
# Check if evolution key is present in the input dictionary
if evo_key in inp_dict:
# If it is, add its corresponding item to the output dictionary
out_dict[evo_key] = inp_dict[evo_key]
# If no matching data variables were found, raise a KeyError
if not out_dict:
raise KeyError(
"No matches found between keys in time evolution file and input dictionary!"
)
return out_dict
[docs]
def calc_norm(evo_dict: dict, ei_inv_dict: dict) -> dict:
"""Calculate normalization factors, either by fuel use only.
Or combined with evolution of emission index of species.
Args:
evo_dict (dict): Evolution,
``{"fuel": numpy.ndarray, "EI_CO2": numpy.ndarray, ..}``
ei_inv_dict (dict): Emission indices, calculated from inventories,
``{"fuel": numpy.ndarray, "EI_CO2": numpy.ndarray, ..}``
Returns:
dict: Dictionary of normalization factors, keys are "fuel" and species,
``{"fuel": numpy.ndarray (norm_fuel = evo_fuel / inv_fuel),
"CO2": numpy.ndarray (norm_fuel * evo_EI / inv_EI), ..}``
"""
norm_dict = {}
# Translation table from evolution to inventory variable names
key_table = KEY_TABLE
evo_fuel = evo_dict["fuel"]
inv_fuel = ei_inv_dict["fuel"]
norm_fuel = np.divide(evo_fuel, inv_fuel)
for key, inv_emi_index in ei_inv_dict.items():
if key == "fuel":
norm_arr = norm_fuel
else:
evo_emi_index = evo_dict[key]
norm_spec = np.divide(evo_emi_index, inv_emi_index)
norm_arr = np.multiply(norm_fuel, norm_spec)
norm_key = key_table[key]
norm_dict[norm_key] = norm_arr
return norm_dict
[docs]
def filter_to_inv_years(
inv_years: np.ndarray, time_range: np.ndarray, interp_dict: dict
) -> dict:
"""Filters dictionary of interpolated arrays to items for inventory years only.
Args:
inv_years (numpy.ndarray): Array of inventory years
time_range (numpy.ndarray): time_range from config
interp_dict (dict): Dictionary of arrays, interpolated over time_range,
keys are e.g. evolution keys or species names
Returns:
dict: Dictionary of arrays, filtered to inventory years, same keys as
interp_dict
"""
mask = np.isin(time_range, inv_years)
filtered_dict = {key: arr[mask] for key, arr in interp_dict.items()}
return filtered_dict
[docs]
def norm_inv(inv_dict: dict, norm_dict: dict) -> dict:
"""Multiply data variables in emission inventories by normalization factors.
Depending on available keys in norm_dict, data variables are multiplied by
norm_fuel or (norm_fuel * norm_EI), with norm_EI = evo_EI / inv_EI.
Args:
inv_dict (dict): Dictionary of xarray Datasets, keys are years of
inventories
norm_dict (dict): Dictionary of normalization factors, keys are
"fuel" and species
Returns:
dict: Dictionary of xarray Datasets (normalized emission inventories),
keys are years of inventories
"""
# Initialize output inventory dictionary
out_inv_dict = {}
# Array index corresponding to inventory years
i = 0
for year, inv in inv_dict.items():
# Get global inventory attributes
inv_attrs = inv.attrs
# Initialize output inventory
out_inv = xr.Dataset()
# Create normalization sub dictionary for current inventory, with scalar values
norm_sub_dict = {}
for key, norm_arr in norm_dict.items():
norm_sub_dict[key] = norm_arr[i]
# Iterate over data variables in inventory
for data_key, data_arr in inv.items():
# Get attributes of data variable
data_attrs = data_arr.attrs
# Check if data_key is within norm_inv_dict
if data_key in norm_sub_dict:
# fuel: multiply with norm_fuel
if data_key == "fuel":
# Multiply inv.fuel by normalization factor norm_fuel
scaled_arr = data_arr * norm_sub_dict["fuel"]
# species: multiply by (norm_fuel * evo_EI / inv_EI)
else:
scaled_arr = data_arr * norm_sub_dict[data_key]
# lon, lat, plev: do NOT multiply
elif data_key in ["lon", "lat", "plev", "ac"]:
scaled_arr = data_arr
# species not in norm_inv_dict: multiply with norm_fuel
else:
scaled_arr = data_arr * norm_sub_dict["fuel"]
# Add data variable to output inventory
scaled_arr.attrs = data_attrs
out_inv = out_inv.merge({data_key: scaled_arr})
# Set global inventory attributes
out_inv.attrs = inv_attrs
out_inv_dict[year] = out_inv
i = i + 1
return out_inv_dict
[docs]
def norm_inventories(config: dict, inv_dict: dict) -> dict:
"""Applies normalization to a dictionary of emission inventories.
This function first interpolates evolution data variables to time_range
from config. It then gets inventory years, calculates inventory sums and
emission indices (dictionaries with spec keys and arrays over inventory
years). The emission indices dictionary is filtered to those species
specified in time evolution. The arrays in evolution data are filtered
to inventory years only. Next, the multipliers used for normalization
are calculated. Finally, the normalization is performed by multiplying
the inventory data variables with the normalization factors.
Args:
config (dict): Configuration dictionary
inv_dict (dict): Dictionary of xarray Datasets, keys are years of
inventories
Returns:
dict: Dictionary of normalized emission inventories, keys are years
of inventories
"""
# Interpolate evolution data variables to time_range from config
time_range, evo_interp_dict = interp_evolution(config)
# Get inventory years, calculate inventory sums and emission indices
# (dictionaries with spec keys and arrays over inventory years)
inv_years, _inv_sum_dict, ei_inv_dict = calc_inv_quantities(config, inv_dict)
# Filter emission indices dictionary to those species specified in time evolution
ei_inv_dict = filter_dict_to_evo_keys(config, ei_inv_dict)
# Filter arrays in evolution data to inventory years only
evo_filtered_dict = filter_to_inv_years(inv_years, time_range, evo_interp_dict)
# Calculate multipliers used for normalization (keys "fuel" and species)
norm_dict = calc_norm(evo_filtered_dict, ei_inv_dict)
# Perform actual normalization: multiply inventory data by norm factors
out_inv_dict = norm_inv(inv_dict, norm_dict)
return out_inv_dict
[docs]
def scale_inv(inv_dict: dict, scale_dict: dict) -> dict:
"""Multiply data variables in emission inventories by scaling factors.
Args:
inv_dict (dict): Dictionary of xarray Datasets, keys are years of inventories
scale_dict (dict): Dictionary of scaling factors {"scaling": np.ndarray}
Returns:
dict: Dictionary of xarray Datasets (scaled emission inventories),
keys are years of inventories
"""
# Get array with scaling multipliers
scale_arr = scale_dict["scaling"]
# Initialize output inventory dictionary
out_inv_dict = {}
# Array index corresponding to inventory years
i = 0
for year, inv in inv_dict.items():
# Get global inventory attributes
inv_attrs = inv.attrs
# Initialize output inventory
out_inv = xr.Dataset()
# Iterate over data variables in inventory
for data_key, data_arr in inv.items():
# Get attributes of data variable
data_attrs = data_arr.attrs
# lon, lat, plev: do NOT multiply
if data_key in ["lon", "lat", "plev", "ac"]:
scaled_arr = data_arr
# multiply fuel, species emissions, and distance by scaling multiplier
else:
scaled_arr = data_arr * scale_arr[i]
# Add data variable to output inventory
scaled_arr.attrs = data_attrs
out_inv = out_inv.merge({data_key: scaled_arr})
# Set global inventory attributes
out_inv.attrs = inv_attrs
out_inv_dict[year] = out_inv
i = i + 1
return out_inv_dict
[docs]
def scale_inventories(config: dict, inv_dict: dict) -> dict:
"""Applies scaling to a dictionary of emission inventories.
This function first interpolates evolution data variables to time_range
from config. It then gets inventory years. The array in evolution data
is filtered to inventory years only. Finally, the scaling is performed
by multiplying the inventory data variables with the scaling factors in
the filtered evolution data.
Args:
config (dict): Configuration dictionary
inv_dict (dict): Dictionary of xarray Datasets, keys are years of
inventories
Returns:
dict: Dictionary of scaled emission inventories, keys are years of
inventories
"""
# Interpolate evolution data variables to time_range from config
time_range, evo_interp_dict = interp_evolution(config)
# Get inventory years
inv_years = np.array(list(inv_dict))
# Filter scaling array in evolution data to inventory years only
evo_filtered_dict = filter_to_inv_years(inv_years, time_range, evo_interp_dict)
# Perform actual scaling: Multiply inventory data variables by scaling factors
out_inv_dict = scale_inv(inv_dict, evo_filtered_dict)
return out_inv_dict