import numpy as np
import xarray as xr
import glob
import datetime
import os
import pandas as pd
import itertools
from .....utils.classes.setup import Setup
from .....utils.classes.domains import Domain
from logging import info, debug
from .....utils.check.errclass import CifError
from .....utils.hdf5 import _hdf5_lock
[docs]
def get_domain(ref_dir, ref_file, input_interval, target_dir, tracer=None):
"""Build the horizontal (and, for point sources, vertical) domain.
Dispatches to :func:`get_area_domain` for the default gridded area-source
mode, or to :func:`get_point_domain` when ``tracer.point_sources`` is set.
Args:
ref_dir (str): directory where the original files are found.
ref_file (str): (template) name of the original files.
input_interval (list): simulation interval, as a list of the two
bounding dates.
target_dir (str): directory where links to the original files are
created.
tracer: the tracer Plugin, corresponding to the paragraph
:bash:`datavect/components/fluxes/parameters/my_species` in the
configuration yaml.
Returns:
Domain: a domain class object, gridded (area mode) or unstructured
(point-source mode).
"""
if not tracer.point_sources:
domain = get_area_domain(tracer)
else:
domain = get_point_domain(tracer)
return domain
[docs]
def get_area_domain(tracer):
"""Build a regular gridded domain from a TNO/CAMS-REG reference file.
Reads ``longitude``/``latitude`` and their ``_bounds`` variables from the
first fetched input file to build the grid centers and corners, with a
single dummy vertical level (surface only).
Args:
tracer: the tracer Plugin holding the list of fetched input files
(``tracer.input_files``).
Returns:
Domain: a gridded, single-level domain.
Raises:
CifError: if no reference file could be found among the fetched
input files.
"""
# Looking for a reference file to read lon/lat in
domain_file = None
ref_file = list(itertools.chain.from_iterable(tracer.input_files.values()))
if len(ref_file) != 0:
domain_file = ref_file[0][0]
if domain_file is None:
raise CifError(
"TNO domain could not be initialized as no file was found"
)
# Read lon/lat in
with _hdf5_lock:
nc = xr.open_dataset(domain_file, decode_times=False)
llon = nc['longitude'].values
llat = nc['latitude'].values
llonb = nc['longitude_bounds'].values
llatb = nc['latitude_bounds'].values
# compute the corner matrix
llonc = np.append(llonb[:, 0], llonb[-1, 1])
llatc = np.append(llatb[:, 0], llatb[-1, 1])
lon, lat = np.meshgrid(llon, llat)
lonc, latc = np.meshgrid(llonc, llatc)
nlat, nlon = lat.shape[0], lat.shape[1]
# If no vertical dimension for emissions, provide dummy vertical
punit = "Pa"
nlevs = 1
sigma_a_mid = np.array([0])
sigma_b_mid = np.array([1])
# Put it to a domain Plugin
domain = Domain(nlon=nlon, nlat=nlat,
zlon=lon, zlat=lat,
zlonc=lonc, zlatc=latc,
nlev=nlevs, pressure_unit=punit,
sigma_b_mid=sigma_b_mid, sigma_a_mid=sigma_a_mid)
return domain
[docs]
def get_point_domain(tracer):
"""Build an unstructured point-source domain from TNO/CAMS-REG files.
Reads the TNO vertical height-distribution profile CSV to determine the
vertical levels, then loops over all fetched input files to extract
non-zero point-source emission locations (filtered by emission category
and source type), building an unstructured domain where each point is
associated with the file it came from (``value_file``).
Args:
tracer: the tracer Plugin, giving access to ``dir_profiles``,
``cat_select``, and the fetched ``input_dates``/``input_files``.
Returns:
Domain: an unstructured domain with one horizontal "cell" per
point source and ``nlev`` vertical height bins.
"""
# Get vertical extent of point sources
VerticalP_file = tracer.dir_profiles + '/TNO_height-distribution_GNFR.csv'
coef = pd.read_csv(VerticalP_file, sep=';', comment='#')
l = list(coef.columns)
info(l)
height = l[l.index('GNFR_Category_Name') + 1:]
height_down = np.array([float(x.split('-')[0].replace(" ", ""))
for x in height])
height_top = np.array(
[float(x.split('-')[1].replace("m", "").replace(" ", ""))
for x in height])
nlevs = len(height_top)
info(f'Top heighs: {height_top}')
# Now loop over files to get lon/lat of point sources
dates = tracer.input_dates
files = tracer.input_files
varnames = tracer.varname
TNO_array = []
nc = None
opened_file = ""
for dd_ref in dates:
for ddi, ff in zip(dates[dd_ref], files[dd_ref]):
# Open file if not already processed
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"]
ds["lon"] = nc["longitude_source"]
ds["lat"] = nc["latitude_source"]
if tracer.cat_select:
cat_list = tracer.cat_select
else:
cat_list = np.arange(
ds['cat_index'].min(), ds['cat_index'].max() + 1)
else:
continue
# Loop over point sources
ds_cat = ds[(ds.cat_index.isin(cat_list)) & (ds.source_type == 2)]
mask_nonzero = ds_cat['emis'] != 0
ds_cat = ds_cat.loc[mask_nonzero]
lons = ds_cat["lon"]
lats = ds_cat["lat"]
tmp_ds = pd.DataFrame({
"lon": lons, "lat": lats,
"file": len(lats) * [ff[0]],
})
tmp_ds["file"] = tmp_ds["file"].astype("category")
TNO_array.append(tmp_ds)
TNO_array = pd.concat(TNO_array)
# Now prepare variables for the domain
zlon = TNO_array["lon"].values[np.newaxis, :]
zlat = TNO_array["lat"].values[np.newaxis, :]
zlonc = zlon[:, :]
zlatc = zlat[:, :]
nlon = len(TNO_array["lon"])
nlat = 1
# Put it to a domain Plugin
domain = Domain(
nlon=nlon, nlat=nlat,
zlon=zlon, zlat=zlat,
zlonc=zlonc, zlatc=zlatc,
nlev=nlevs,
heights=height_top,
height_unit="m",
unstructured_domain=True,
value_file=TNO_array["file"],
)
return domain
