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 TNO domain, dispatching to a gridded or point-source layout.
Delegates to `get_area_domain` (regular grid) or `get_point_domain`
(unstructured, one cell per emission point) depending on
`tracer.point_sources`.
Args:
--------
ref_dir: directory where the original files are found (unused
directly; the reference file is taken from `tracer.input_files`)
ref_file: (template) name of the original files (unused directly)
input_interval: list of the periods to simulate, each item is the list of the dates of
the period (unused directly)
target_dir: directory where the links to the orginal files are created
(unused directly)
tracer: the flux tracer plugin, providing ``point_sources`` and
``input_files``
Ouputs:
--------
setup of the domain in section "Initializes domain"
"""
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, surface-only domain from a reference TNO file.
Reads ``longitude``/``latitude``/``longitude_bounds``/``latitude_bounds``
from the first fetched reference file, and builds a `Domain` with a
single dummy vertical level.
Args:
tracer: The flux tracer plugin, providing ``input_files``.
Returns:
Domain: a regular gridded, single-level `Domain`.
Raises:
CifError: If no reference file is available.
"""
# 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 all fetched TNO files.
Reads the vertical height-distribution CSV (`TNO_height-distribution_GNFR.csv`)
to get the level tops, then loops over all files listed in
`tracer.input_files`, collecting the (category- and source-type-
filtered) non-zero point-source emission locations from each. Builds an
unstructured `Domain` with one "cell" per point and a ``value_file``
attribute mapping each point to its originating file.
Args:
tracer: The flux tracer plugin, providing ``dir_profiles``,
``input_dates``, ``input_files``, ``varname`` and
``cat_select``.
Returns:
Domain: an unstructured domain with one cell per point source.
"""
# 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)
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
