from logging import warning
import xarray as xr
import pandas as pd
import numpy as np
import os
from ......utils.hdf5 import _hdf5_lock
[docs]
def create_oem_tv_scaling_factors(self, ddi, ddf, data, oem_dir, tfactors_oem_group):
"""Write the OEM temporal scaling-factor file for one sub-period.
When ``self.use_hourofyear=True``, extracts the hourly temporal profiles
for the period ``[ddi, ddf]`` from the full-year profile file and writes
a reduced ``hour_of_year.nc`` to *oem_dir* that ICON-ART reads at
runtime.
Args:
self: ICON-ART model plugin instance.
ddi (datetime): period start date.
ddf (datetime): period end date.
data (xr.Dataset): CIF flux data for the period.
oem_dir (str): OEM output directory.
tfactors_oem_group: OEM temporal-factor group object.
"""
# Create and update the OEM hour-of-year scaling factors
if self.use_hourofyear:
t_profiles_file = os.path.join(oem_dir, f'hour_of_year.nc')
hourly_time = pd.date_range(
data.time[0].values,
data.time[-1].values + pd.Timedelta(self.input_resolution),
freq="1h")
with _hdf5_lock:
ds_tf = xr.open_dataset(t_profiles_file) \
if os.path.exists(t_profiles_file) else xr.Dataset()
ds_tf[tfactors_oem_group] = xr.DataArray(
np.ones((hourly_time.shape[0], data.shape[1])),
coords=[np.arange(hourly_time.shape[0]), np.arange(data.shape[1])],
dims=['hourofyear', 'country']
)
def relative_hour_of_year(current_date):
return int((current_date - ddi).total_seconds() // 3600)
hourofyear_idxs = [relative_hour_of_year(
dd) for dd in hourly_time] # TODO: check its fine
data_hourly = data.interp(time=hourly_time, kwargs={
"fill_value": "extrapolate"})
tf = data_hourly / data.mean(dim='time')
tf = tf.fillna(1) # TODO: Is there really any nan left ?
ds_tf[tfactors_oem_group][hourofyear_idxs, :] = tf.data
if os.path.exists(t_profiles_file):
os.remove(t_profiles_file)
with _hdf5_lock:
ds_tf.to_netcdf(t_profiles_file)
else:
warning("Careful, deriving hourofday, dayofweek and "
"monthofyear temporal files is only relevant "
"if the fluxes follow such a temporal periodicity.\n"
"Prefer the 'use_hourofyear' option if you think it "
"is not the case.")
for time_type, time_len, file_name, time_variable in zip(
['hourofday', 'dayofweek', 'monthofyear'],
[24, 7, 12],
['hour_of_day', 'day_of_week', 'month_of_year'],
['hour', 'dayofweek', 'month']
):
t_profiles_file = os.path.join(oem_dir, f'{file_name}.nc')
with _hdf5_lock:
ds_tf = xr.open_dataset(t_profiles_file) \
if os.path.exists(t_profiles_file) else xr.Dataset()
if tfactors_oem_group in ds_tf:
continue
ds_tf[tfactors_oem_group] = xr.DataArray(
np.ones((time_len, data.shape[1])),
coords=[np.arange(time_len), np.arange(data.shape[1])],
dims=[time_type, 'country']
)
tf = data.groupby(f'time.{time_variable}').mean()
if len(tf[time_variable]) == time_len:
tf = tf / tf.mean(dim=time_variable)
tf = tf.fillna(1)
ds_tf[tfactors_oem_group][:] = tf.data
if os.path.exists(t_profiles_file):
os.remove(t_profiles_file)
with _hdf5_lock:
ds_tf.to_netcdf(t_profiles_file)
# Create the vertical scaling factors files
v_profiles_file = os.path.join(oem_dir, 'vertical_profile.nc')
ds_v = xr.Dataset(
data_vars={
'layer_bot': ('level', [0., 20., 92., 184., 324., 522., 781.]),
'layer_mid': ('level', [10., 56, 138, 254, 423, 651.5, 943.5]),
'layer_top': ('level', [20., 92., 184., 324., 522., 781., 1106.]),
'ground_emissions': ('level', [1., 0., 0., 0., 0., 0., 0.]),
# Ground emissions are all emitted on the smallest level
},
)
if os.path.exists(v_profiles_file):
os.remove(v_profiles_file)
with _hdf5_lock:
ds_v.to_netcdf(v_profiles_file)
