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
[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
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('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}')
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"),
)
