import copy
import datetime
import glob
import os
from logging import debug
import numpy as np
import pandas as pd
from ......utils.parallel import thread
from ...utils.read import read_flexpart_gridinit
from ......utils.check.errclass import CifError
[docs]
def inicond_contribution(
self, mode, dataobs,
fp_header_init, spec, ddi,
batch_computation=False
):
"""Compute the initial-condition–concentration contribution for one species and period.
Convolves the FLEXPART initial-condition sensitivity fields
(``grid_initial_*``) with the stored initial-condition array to produce
the background concentration increment for each observation in *dataobs*.
Args:
self: Lagrangian model plugin instance.
mode (str): ``'fwd'`` or ``'tl'``.
dataobs: CIF observation data-store for the period.
fp_header_init: FLEXPART header for the initial-condition sensitivity.
spec (str): species name.
ddi (datetime): sub-simulation period start.
batch_computation (bool): ensemble batch mode flag.
Returns:
dataobs updated with ``'spec'`` contributions from initial conditions.
"""
inicond = {
spec: {"spec": self.datainicond[("inicond", spec)][ddi]["spec"]}}
if mode == "tl" and "incr" in self.datainicond[("inicond", spec)][ddi]:
inicond[spec]["incr"] = \
self.datainicond[("inicond", spec)][ddi]["incr"]
# Initialize debug columns if needed
if self.dump_debug:
dataobs[spec][("flexpart", "background")] = 0.
dataobs[spec][("flexpart", "background_tl")] = 0.
# Specific treatment for batch computing
sample_species = [spec]
if batch_computation:
ref_species = self.perturbed_species[spec]
sample_species = [
s for s in self.chemistry.acspecies.attributes
if self.perturbed_species[s] == ref_species
]
inicond = {
s: {"spec": self.datainicond[("inicond", s)][ddi]["spec"]}
for s in sample_species
}
if mode == "tl" and "incr" in self.datainicond[("inicond", spec)][ddi]:
for s in sample_species:
inicond[s]["incr"] = self.datainicond[(
"inicond", s)][ddi]["incr"]
if self.dump_debug:
for s in sample_species:
dataobs[s][("flexpart", "background")] = 0.
dataobs[s][("flexpart", "background_tl")] = 0.
# Execute parallel threads
nthreads = self.nthreads
nobs = len(dataobs[sample_species[0]])
thread_intervals = np.linspace(0, nobs, nthreads + 1).astype(int)
@thread
def thread_function(ithread):
for obs_i in range(thread_intervals[ithread], thread_intervals[ithread + 1]):
process_obs_row(
self, dataobs, ithread, fp_header_init, inicond, obs_i,
sample_species, batch_computation
)
thread_function(range(nthreads))
# Flush inicond
for s in sample_species:
self.datainicond[("inicond", s)][ddi]["spec"] = None
self.datainicond[("inicond", s)][ddi]["incr"] = None
[docs]
def process_obs_row(self, dataobs, ithread,
fp_header_init,
inicond, obs_i,
sample_species, batch_computation):
"""Compute the initial-condition contribution for a single observation row.
Reads the initial-condition sensitivity field for the observation at
*obs_i* and convolves it with *inicond* to produce the background
concentration contribution.
Args:
self: Lagrangian model plugin instance.
dataobs (dict): species-keyed observation data-stores.
ithread (int): thread index.
fp_header_init: FLEXPART initial-condition sensitivity header.
inicond (dict): initial-condition arrays keyed by species.
obs_i (int): row index into the observation data-store.
sample_species (list): species names to process.
batch_computation (bool): ensemble batch mode flag.
"""
ref_spec = sample_species[0]
row = dataobs[ref_spec]["metadata"].iloc[obs_i]
station = row.station
network = row.network
subdir = row.date.strftime(self.footprint_dir_format)
# Translate station name if needed
if hasattr(self, "station_name_dict"):
station = self.dict_station_name[station.upper()].upper()
# Infer folder structure
runsubdir_init = os.path.join(
self.run_dir_bg,
self.footprint_stat_subdir_format.format(
stat=station.upper(), network=network),
subdir
)
release_date = row.date - pd.to_timedelta(self.release_shift)
file_date = release_date.strftime('%Y%m%d%H%M%S')
file_name = self.file_ini_format.format(
date=file_date, stat=station.upper(), network=network)
list_valid = glob.glob(os.path.join(runsubdir_init, file_name))
if list_valid == []:
debug(f"WARNING: file not found: {os.path.join(runsubdir_init, file_name)}")
return
elif len(list_valid) > 1:
raise CifError(
f"Multiple files fit the specified format {self.file_ini_format}. "
f"This can be related to the use of a wildcard... "
f"Please check your yml"
)
file_name = os.path.basename(list_valid[0])
if self.preloaded_footprints.get("loaded_file", "") == list_valid[0] \
and self.reload_footprints:
debug("Using preloaded concentrations from "
f"{file_name} for station {station}")
for s in sample_species:
sample_inicond = inicond[s]
for data_id in sample_inicond:
dataobs[s].iloc[
obs_i, dataobs[s].columns.get_loc(("maindata", data_id))
] = self.preloaded_footprints[data_id][s]
else:
debug(f"Thread #{ithread}: Reading {file_name} for station {station}")
grid_init = read_flexpart_gridinit(
runsubdir_init, file_name, fp_header_init)
if self.reload_footprints:
self.preloaded_footprints = {
"loaded_file": list_valid[0],
"data": copy.deepcopy(grid_init),
"spec": {},
"incr": {}
}
# Normalize grid_init to make sure that total sensitivity is 1
grid_init /= grid_init.sum()
# Multiply 3-D sensitivity to background concentrations
# WARNING: do not deal with temporal and vertical dimension yet
nz = (fp_header_init.outheight != 0.).sum()
ini_sensit = grid_init.T.reshape(nz, -1)
for s in sample_species:
sample_inicond = inicond[s]
for data_id in sample_inicond:
dataini = sample_inicond[data_id]
ini_dates = pd.DatetimeIndex(
dataini.time.values).to_pydatetime()
inds_inicond = np.argmin(
(np.array(
[row.date - pd.to_timedelta(self.backward_trajdays)]
)[:, np.newaxis] - ini_dates[np.newaxis, :]
) >= datetime.timedelta(0), axis=1) - 1
istartsensit = (
0 if self.domain.zlon.size == self.domain.zlon_in.size
else self.domain.zlon_in.size
)
ini_sim = (
dataini[inds_inicond, :, 0,
istartsensit:].values * ini_sensit
).sum()
# Filling simulation
dataobs[s].iloc[
obs_i, dataobs[s].columns.get_loc(("maindata", data_id))
] += ini_sim
# Dump debug
if self.dump_debug:
if data_id == "spec":
sim_col = "background"
else:
sim_col = "background_tl"
dataobs[s].iloc[
obs_i, dataobs[s].columns.get_loc(
("flexpart", sim_col))
] = ini_sim
# Save simulation for later if needed
if self.reload_footprints:
self.preloaded_footprints[data_id][s] = ini_sim
