import calendar
from logging import debug,info
from pathlib import Path
import numpy as np
import pandas as pd
import xarray as xr
[docs]
def read(
self,
name,
varnames,
dates,
files,
interpol_flx=False,
comp_type=None,
tracer=None,
**kwargs,
):
"""Get fluxes from raw files and load them into a pyCIF
variables
"""
# Time profiles
Pdir = tracer.dir_profiles
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')
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('Reading of {} in {} for {}'.format([varnames], ff, ddi))
opened_file = ff
ds = pd.DataFrame({})
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}')
idx = np.where(profil_cat_list==list_cat_name[cat-1])
if idx[0].size == 0:
debug(f" !!! use default profile {tracer.default_profil} : to check !!!")
catinprofiles = tracer.default_profil
else:
catinprofiles = idx[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] / 1e6) #from /km2 to /m2
data.append(TNO_array)
outdate.append(ddi[0])
idate += 1
return xr.DataArray(
data,
coords={"time": outdate},
dims=("time", "lev", "lat", "lon"),
)
