import datetime
import os
from pathlib import Path
import numpy as np
import xarray as xr
import pandas as pd
import string
import calendar
import pytz
from logging import info, debug
from .....utils.datastores.empty import init_empty
from .....utils.hdf5 import _hdf5_lock
[docs]
def read(
self,
name,
varnames,
dates,
files,
interpol_flx=False,
comp_type=None,
tracer=None,
**kwargs
):
"""Get TNO fluxes and apply GNFR time (and, for point sources, vertical) profiles.
Loads the month/day/hour (and, for point sources, height) GNFR profile
CSVs from `tracer.dir_profiles`. Then, for each requested date/file
(re-reading category indices, source types and emission values only
when the file changes), loops over the active emission categories
(`tracer.cat_select`, or all categories found in the file), matches
each to its GNFR temporal profile, and applies the corresponding
month-of-year, day-of-week and hour-of-day scaling coefficients to the
raw category emission (and, for point sources, the height-distribution
coefficients to spread across vertical levels). Results are
accumulated into a ``(lev, lat, lon)`` array per date, either scattered
onto point-source domain cells or gridded area cells (area-normalized
using the file's ``area`` variable).
Args:
self: The flux tracer plugin instance (unused directly; `tracer` is
used instead).
name: Unused directly, kept for interface consistency with other
flux plugins.
varnames: Name of the emission variable to read in the file.
dates: list of ``[start, end]`` date intervals to extract; only the
start of each interval is used to compute the local
month/weekday/hour.
files: list of files matching `dates`.
interpol_flx (bool): Unused, kept for interface consistency.
comp_type: Unused, kept for interface consistency.
tracer: The flux tracer plugin, providing ``dir_profiles``,
``point_sources``, ``cat_select`` and ``domain``.
Returns:
xr.DataArray: the flux data with dimensions
``(time, lev, lat, lon)``.
"""
# Time profiles
Pdir = tracer.dir_profiles
info(Pdir)
TPmonth_file = Path(Pdir, 'timeprofiles-month-in-year_GNFR.csv')
TPday_file = Path(Pdir, 'timeprofiles-day-in-week_GNFR.csv')
TPhour_file = Path(Pdir, 'timeprofiles-hour-in-day_GNFR.csv')
VerticalP_file = Path(Pdir, 'TNO_height-distribution_GNFR.csv')
info(TPmonth_file)
coef_dict = {}
for key, f in zip(['month', 'day', 'hour', 'height'],
[TPmonth_file, TPday_file, TPhour_file, VerticalP_file]):
coef = pd.read_csv(f, sep=';', comment='#')
l = list(coef.columns)
l2 = l[l.index('GNFR_Category_Name') + 1:]
coef_dict[key] = coef[l2].values
profil_cat_list = coef['TNO GNFR sectors Sept 2018'].values
nlevel = len(coef_dict['height'][0])
if not tracer.point_sources:
nlevel = 1
debug(f"List of categories: {profil_cat_list}")
debug(f"Number of levels: {nlevel}")
# list of the various fields read
data = []
outdate = []
nc = None
opened_file = ""
idate = 0
for ddi, ff in zip(dates, files):
# dd_UTC = ddi.tz_localize('UTC')
# dd_CET = dd_UTC.tz_convert('Europe/Berlin')
dd_CET = ddi[0]
yr = dd_CET.year
mm = dd_CET.month
dd = dd_CET.weekday()
hh = dd_CET.hour
if ff != opened_file or nc is None:
debug(f'Reading of {[varnames]} in {ff} for {ddi}')
opened_file = ff
ds = pd.DataFrame({})
with _hdf5_lock:
nc = xr.open_dataset(ff[0], decode_times=False)
list_cat_name = nc['emis_cat_code'].values
list_cat_name = [b.decode("utf-8") for b in list_cat_name]
debug(f"List of category names: {list_cat_name}")
ds['cat_index'] = nc['emission_category_index'].values
ds['source_type'] = nc['source_type_index'].values
ds['emis'] = nc[varnames].values
nlon = nc.dims["longitude"]
nlat = nc.dims["latitude"]
if tracer.point_sources:
ds["lon"] = nc["longitude_source"]
ds["lat"] = nc["latitude_source"]
else:
ds['ilon'] = nc['longitude_index'].values
ds['ilat'] = nc['latitude_index'].values
areas = nc['area'].values
ds = ds.loc[ds['ilat'] > 0]
ds = ds.loc[ds['ilon'] > 0]
if tracer.cat_select:
cat_list = tracer.cat_select
else:
cat_list = np.arange(
ds['cat_index'].min(), ds['cat_index'].max() + 1)
# Mask of locations according to file and categories
TNO_array = np.zeros((nlevel, nlat, nlon))
if tracer.point_sources:
mask_domain = tracer.domain.value_file == ff[0]
ds = ds[(ds.cat_index.isin(cat_list)) & (ds.source_type == 2)]
mask_nonzero = ds['emis'] != 0
ds = ds.loc[mask_nonzero]
TNO_array = np.zeros((nlevel, tracer.domain.nlat, tracer.domain.nlon))
# Looping over categories
for cat in cat_list:
debug(f'Processing category: {cat}')
catinprofiles = np.where(profil_cat_list==list_cat_name[cat-1])[0][0]
debug(f"Category names: {list_cat_name[cat-1]}")
debug(f"Category profiles: {catinprofiles}")
local_month_coef = coef_dict['month'][catinprofiles][mm - 1]
local_nbdays = \
calendar.mdays[mm] + (mm == 2 and calendar.isleap(yr))
debug(f"Daily profiles: {coef_dict['day'][catinprofiles]}")
local_wday_coef = coef_dict['day'][catinprofiles][dd]
local_hour_coef = np.array(coef_dict['hour'][catinprofiles])[hh]
local_height_coef = np.array(coef_dict['height'][catinprofiles])
if tracer.point_sources:
# Filtering category emissions
mask_loc = ds["cat_index"] == cat
ds_cat = ds.loc[mask_loc]
target_mask = np.where(mask_domain)[0][mask_loc]
# Computing emissions for this category
emis = ds_cat['emis'].values
emis_hour = (emis / 12 * local_month_coef
/ local_nbdays * local_wday_coef
/ 24 * local_hour_coef)
# Factorizing by levels
emis_level = emis_hour[np.newaxis] * local_height_coef[:, np.newaxis]
TNO_array[:, :, target_mask] += emis_level[:, np.newaxis, :]
else:
mask_cat = (ds.cat_index == cat) & (ds.source_type == 1)
ds_cat = ds.loc[mask_cat]
emis = ds_cat['emis'].values
emis_hour = (emis / 12 * local_month_coef
/ local_nbdays * local_wday_coef
/ 24 * local_hour_coef)
ilons = ds_cat['ilon']
ilats = ds_cat['ilat']
np.add.at(TNO_array[0],
(ilats - 1, ilons - 1),
emis_hour / areas[ilats - 1, ilons - 1] )
data.append(TNO_array)
outdate.append(ddi[0])
idate += 1
return xr.DataArray(
data,
coords={"time": outdate},
dims=("time", "lev", "lat", "lon"),
)
