import os
import glob
import copy
import xarray as xr
import numpy as np
import pandas as pd
import time as tm
import multiprocessing as mp
from logging import debug
from itertools import repeat
from ..utils import (
OUTPUT_DIRNAME,
OUTPUT_PATTERN,
)
from ......utils.hdf5 import _hdf5_lock
try:
from sklearn.neighbors import BallTree
except ImportError:
pass
from logging import debug, warning, info
# ---------------------------------------------------------
# -- CONSTANTS
# ---------------------------------------------------------
R_EARTH = 6371 # Radius of earth [km]
N_AVO = 6.022e23 # Avogadro number
M_DRYAIR = 28.97e-3 # Molecular weight of dry air [kg/mol]
# ---------------------------------------------------------
# -- SIDE FUNCTIONS
# ---------------------------------------------------------
[docs]
def reduce_size(file, output_path, cells, levs, nlev, ncells, model_timestep):
""" Reduce the size of the file by selecting relevant cells
and levels only.
Args:
file (str): Path of the file
output_path (str): Path of the output directory
cells (list): List of indexes to cells
levs (list): List of indexes to levels
nlev (int): Number of levels
ncells (int): Number of cells
model_timestep (float): Timestep of the model
"""
# Open the dataset
with _hdf5_lock:
ds = xr.open_dataset(file)
# Find name of relevant dimensions
dim_height = [dim for (dim, ndim) in ds.dims.items() if ndim == nlev][0]
dim_cell = [dim for (dim, ndim) in ds.dims.items() if ndim == ncells][0]
# Build hlay [m]
dim_height_ifc = [dim for dim in ds["z_ifc"].dims if "height" in dim][0]
hlay = -ds["z_ifc"].diff(dim=dim_height_ifc)
hlay = hlay.rename({dim_height_ifc: dim_height})
hlay = hlay.expand_dims({"time": ds["time"]})
hlay[dim_height] = ds[dim_height]
ds["hlay"] = hlay
# Build dpressure [Pa]
dim_pres_ifc = [dim for dim in ds["pres_ifc"].dims if "height" in dim][0]
dpressure = ds["pres_ifc"].diff(dim=dim_pres_ifc)
dpressure = dpressure.rename({dim_pres_ifc: dim_height})
dpressure[dim_height] = ds[dim_height]
ds["dpressure"] = dpressure
# Build airm [molecules/cm3]
ds["airm"] = ds["rho"] * (1 - ds["qv"]) * 1e-6 * N_AVO / M_DRYAIR
# Build pressure [Pa]
ds = ds.rename({"pres": "pressure"})
# Add time dimension to z_mc
ds["z_mc"] = ds["z_mc"].expand_dims({"time": ds["time"]})
# Select relevant indexes
dim_height = [dim for (dim, ndim) in ds.dims.items() if ndim == nlev][0]
dim_cell = [dim for (dim, ndim) in ds.dims.items() if ndim == ncells][0]
vars2fetch = [var for var in list(ds.data_vars) if ds.data_vars[var].dims == ("time", dim_height, "ncells")]
ds = ds[vars2fetch].sel({dim_height: levs, dim_cell: cells})
ds["ncells"] = cells
# Offset the time [trick to get first timestep AFTER restart]
# ds["time"] = ds.time - pd.Timedelta(seconds=model_timestep)
new_date = pd.to_datetime(ds["time"][0].values).strftime("__%Y%m%dT%H%M%SZ")
filename_out = os.path.join(output_path, "reduced", os.path.basename(file).split("__")[0] + new_date + "_short.nc")
# Save the dataset
with _hdf5_lock:
ds.to_netcdf(filename_out, format="NETCDF4_CLASSIC")
[docs]
def concatenate_byday(output_path, files2concatenate):
"""
Concatenate the files by day to make the subsequent processing faster.
"""
with _hdf5_lock:
ds_var = xr.open_mfdataset(files2concatenate, concat_dim="time", combine="nested")
ds_var["time"] = pd.to_datetime(ds_var["time"])
filename_begin = os.path.basename(files2concatenate[0])
filename_out = os.path.join(output_path, "concatenated_byday", f"{filename_begin[:-3]}.nc")
if os.path.isfile(filename_out):
os.remove(filename_out)
with _hdf5_lock:
ds_var.to_netcdf(filename_out, format="NETCDF4_CLASSIC")
# TODO: remove occurences of "height"", not supposed to know the name
# If all observations are processed
[docs]
def apply_interpolation_by_chunk_full(segment_idx_obs, ds_icon, ds_interp, all_trcrs):
"""Apply full 3-D horizontal+vertical interpolation to a chunk of observations.
Identical in semantics to
:func:`apply_interpolation.apply_interpolation_by_chunk_full` but
operates on a reduced ICON output with dimension names ``height`` and
``ncells``.
Args:
segment_idx_obs (tuple[int, int]): (start, end) observation slice.
ds_icon (xr.Dataset): reduced ICON output on the icosahedral grid.
ds_interp (xr.Dataset): pre-computed interpolation metadata.
all_trcrs (list[str]): tracer variable names.
Returns:
xr.Dataset: interpolated concentrations for the observation chunk.
"""
debug(f"Segment {segment_idx_obs} : started...")
idx_obs_start, idx_obs_end = segment_idx_obs
# Get interpolation data for the right indexes
ds_interp_red = ds_interp.isel(index=slice(idx_obs_start, idx_obs_end))
# Fetch interpolation data
ilev_below = ds_interp_red["ilev_below"].astype(int)
ilev_above = ds_interp_red["ilev_above"].astype(int)
iadjcells = ds_interp_red["iadjcells"]
vert_scaling_fact = ds_interp_red["vert_scale_fact"]
inv_dist_adjcells = ds_interp_red["inv_dist_adjcells"]
tsteps = ds_interp_red["tstep"]
parameters = ds_interp_red["parameter"]
# Convert interpolation data to DataArrays
vals_lev_above = ilev_above.values + 1
vals_lev_above = xr.DataArray(vals_lev_above, dims=("obs", "adjcell"))
vals_lev_below = ilev_below.values + 1
vals_lev_below = xr.DataArray(vals_lev_below, dims=("obs", "adjcell"))
vals_cell = iadjcells.values
vals_cell = xr.DataArray(vals_cell, dims=("obs", "adjcell"))
vals_time = ds_icon["time"][tsteps].values
vals_time = xr.DataArray(vals_time, dims=("obs"))
vals_vert_scal = vert_scaling_fact.rename({"index": "obs"})
vals_inv_dist = inv_dist_adjcells.rename({"index": "obs"})
# Apply interpolation to model values
ds_above_adj = ds_icon[all_trcrs].sel(time=vals_time, height=vals_lev_above, ncells=vals_cell)
ds_below_adj = ds_icon[all_trcrs].sel(time=vals_time, height=vals_lev_below, ncells=vals_cell)
ds_adj = ds_below_adj + vals_vert_scal * (ds_above_adj - ds_below_adj)
ds_out = (ds_adj * vals_inv_dist).sum(dim="adjcell") / vals_inv_dist.sum(dim="adjcell")
# Change the index
new_index = ds_interp_red.index.values
ds_out["obs"] = new_index
# Add parameter information
ds_out["parameter"] = xr.DataArray(data=parameters.values, dims="obs", coords={"obs": ds_out.obs})
debug(f"Segment {segment_idx_obs} : completed!")
return ds_out.load()
# If observations are grouped by index
[docs]
def apply_interpolation_by_chunk_reduced(segment_idx_obs, ds_icon, ds_interp, all_trcrs):
"""Apply reduced (level-fixed) horizontal interpolation to a chunk of observations.
Identical in semantics to
:func:`apply_interpolation.apply_interpolation_by_chunk_reduced` but
operates on a pre-reduced ICON output dataset.
Args:
segment_idx_obs (tuple[int, int]): (start, end) observation slice.
ds_icon (xr.Dataset): reduced ICON output on the icosahedral grid.
ds_interp (xr.Dataset): pre-computed interpolation metadata.
all_trcrs (list[str]): tracer variable names.
Returns:
xr.Dataset: interpolated concentrations for the observation chunk.
"""
idx_obs_start, idx_obs_end = segment_idx_obs
ds_interp_red = ds_interp.isel(index=slice(idx_obs_start, idx_obs_end))
# Process data that do not need to be interpolated vertically
mask_alt_nan = np.isnan(ds_interp_red["ilev_below"][:, 0].values)
iadjcells = ds_interp_red["iadjcells"][mask_alt_nan]
inv_dist_adjcells = ds_interp_red["inv_dist_adjcells"][mask_alt_nan]
tsteps = ds_interp_red["tstep"][mask_alt_nan]
parameters = ds_interp_red["parameter"][mask_alt_nan]
vals_cell = iadjcells.values
vals_cell = xr.DataArray(vals_cell, dims=("obs_red", "adjcell"))
vals_time = ds_icon["time"][tsteps].values
vals_time = xr.DataArray(vals_time, dims=("obs_red"))
vals_inv_dist = inv_dist_adjcells.rename({"index": "obs_red"})
# Unroll model values using all levels
min_lev = int(ds_icon["height"].min().values)
max_lev = int(ds_icon["height"].max().values)
nlev = max_lev + 1 - min_lev
ds_adj = ds_icon[all_trcrs].sel(
time=vals_time,
height=slice(min_lev, max_lev + 1),
ncells=vals_cell,
)
ds_adj["obs_red"] = vals_inv_dist.obs_red
# Perform horizontal interpolation
ds_out = (ds_adj * vals_inv_dist).sum(dim="adjcell") / vals_inv_dist.sum(dim="adjcell")
ds_out = ds_out.reindex(height=list(reversed(ds_out.height)))
ds_out = ds_out.stack(obs=("obs_red", "height"), create_index=False)
ds_out = ds_out.drop_vars("height")
# Change the index
new_index = ds_interp_red.index[mask_alt_nan].values
new_index = np.repeat(new_index, (max_lev - min_lev + 1))
ds_out["obs"] = new_index
# Solve duplicates due to time averaging
u, c = np.unique(ds_interp_red.index, return_counts=True)
index_dup = u[c > 1]
if len(index_dup) != 0:
mask_dup = np.isin(new_index, index_dup)
da_out = ds_out.to_array().load()
index_in = np.arange(len(new_index))[mask_dup]
array_in = ds_out.isel(obs=index_in).to_array().values
indices_in = np.vstack(
(
np.hstack(np.split(np.arange(len(index_in)), len(index_in) // nlev)[::2]),
np.hstack(np.split(np.arange(len(index_in)), len(index_in) // nlev)[1::2]),
)
).ravel("F")
da_out[:, index_in] = copy.deepcopy(array_in[:, indices_in])
ds_out = da_out.to_dataset(dim="variable")
# Add parameter information
ds_out["parameter"] = xr.DataArray(
data=np.repeat(parameters.values, (max_lev - min_lev + 1)), dims="obs", coords={"obs": new_index}
)
# Process data that need to be interpolated vertically
if np.sum(~mask_alt_nan):
# Fetch interpolation data
ilev_below = ds_interp_red["ilev_below"][~mask_alt_nan]
ilev_above = ds_interp_red["ilev_above"][~mask_alt_nan]
iadjcells = ds_interp_red["iadjcells"][~mask_alt_nan]
vert_scaling_fact = ds_interp_red["vert_scale_fact"][~mask_alt_nan]
inv_dist_adjcells = ds_interp_red["inv_dist_adjcells"][~mask_alt_nan]
tsteps = ds_interp_red["tstep"][~mask_alt_nan]
parameters = ds_interp_red["parameter"][~mask_alt_nan]
# Convert interpolation data to DataArrays
vals_ilev_above = ilev_above.values
vals_ilev_above = xr.DataArray(vals_ilev_above, dims=("obs", "adjcell"))
vals_ilev_below = ilev_below.values
vals_ilev_below = xr.DataArray(vals_ilev_below, dims=("obs", "adjcell"))
vals_icell = xr.DataArray(iadjcells.values, dims=("obs", "adjcell"))
vals_itime = xr.DataArray(tsteps.values, dims=("obs"))
vals_vert_scal = vert_scaling_fact.rename({"index": "obs"})
vals_inv_dist = inv_dist_adjcells.rename({"index": "obs"})
# Apply interpolation to model values
ds_above_adj = ds_icon[all_trcrs].isel(time=vals_itime, height=vals_ilev_above, ncells=vals_icell)
ds_below_adj = ds_icon[all_trcrs].isel(time=vals_itime, height=vals_ilev_below, ncells=vals_icell)
ds_adj = ds_below_adj + vals_vert_scal * (ds_above_adj - ds_below_adj)
ds_out_with_alt = (ds_adj * vals_inv_dist).sum(dim="adjcell") / vals_inv_dist.sum(dim="adjcell")
ds_out_with_alt["obs"] = ds_interp_red.index[~mask_alt_nan].values
# Add parameter information
ds_out_with_alt["parameter"] = xr.DataArray(data=parameters.values, dims="obs", coords={"obs": ds_out.obs})
# Concatenate the two datasets and sort it with obs index
ds_out = xr.concat([ds_out, ds_out_with_alt], dim="obs")
ds_out = ds_out.sortby("obs")
return ds_out.load()
# ---------------------------------------------------------
# -- MAIN FUNCTION
# ---------------------------------------------------------
[docs]
def process_output(self, runsubdir, ddi):
"""Post-process ICON-ART output: reduce size, concatenate, and interpolate.
1. Reduces each output NetCDF to the observation-relevant cells and
levels using :func:`reduce_size` to limit disk and memory usage.
2. Concatenates daily files using :func:`concatenate_byday`.
3. Applies horizontal + vertical interpolation to all observation
locations using :func:`apply_interpolation_by_chunk_full` or
:func:`apply_interpolation_by_chunk_reduced`.
Results are stored in ``self.sim_data`` for later extraction by
``read_sim``.
Args:
self: ICON-ART model plugin instance.
runsubdir (str): path to the period run directory.
ddi (datetime): sub-simulation period start.
"""
NADJ = self.output_interp_neighbors
output_path = os.path.join(runsubdir, OUTPUT_DIRNAME)
# ---------------------------------------------------------
# -- Initialize interpolation data
# ---------------------------------------------------------
# Get outputs files
list_files = glob.glob(os.path.join(output_path, OUTPUT_PATTERN))
list_files = sorted(list_files)
# Find index of maximum height (indices start at 0)
with _hdf5_lock:
f0 = xr.open_dataset(list_files[0])
da_alt = f0["z_mc"]
# Check obs.csv exists
if not os.path.isfile(f"{runsubdir}/obs.csv"):
warning("No obs.csv file found")
return
# Extract station information
data_obs = pd.read_csv(f"{runsubdir}/obs.csv", na_values=-999, index_col=0)
data_obs["i"] = data_obs["i"].astype(int)
data_obs["tstep"] = data_obs["tstep"].astype(int)
# min_ilevel = int(data_obs['icon_ilevel'].min()) + 1 # height dim is 1..nlev+1
# max_ilevel = int(data_obs['icon_ilevel'].max()) + 1 # height dim is 1..nlev+1
nobs = len(data_obs)
icells_obs = data_obs["i"].values
lons_obs = data_obs["lon"].values
lats_obs = data_obs["lat"].values
tsteps_obs = data_obs["tstep"].values
tracers_obs = data_obs["parameter"].values
index_obs = data_obs.index.values
# Create a Dataset to store all the information.
# icell and ilev starts at 0
# BUT, coord cell starts at 0 and coord height starts at 1
# so icell ~= coord cell but coord height ~= ilev + 1
# initialization to -1
ds_interp = xr.Dataset(
data_vars=dict(
icell=(["index"], icells_obs),
tstep=(["index"], tsteps_obs),
parameter=(["index"], tracers_obs),
iadjcells=(["index", "adjcell"], np.nan * np.ones((nobs, NADJ))),
alt_adjcells=(["index", "adjcell"], np.nan * np.ones((nobs, NADJ))),
ilev_below=(["index", "adjcell"], np.nan * np.ones((nobs, NADJ))),
ilev_above=(["index", "adjcell"], np.nan * np.ones((nobs, NADJ))),
vert_scale_fact=(["index", "adjcell"], np.nan * np.ones((nobs, NADJ))),
inv_dist_adjcells=(["index", "adjcell"], np.nan * np.ones((nobs, NADJ))),
),
coords=dict(index=index_obs, icon_ilevel=range(1, self.domain.nlev + 1)),
)
ds_interp["icell"] = ds_interp["icell"].astype(int)
# ------------------------------------------------------------
# Horizontal interpolation parameters
# ------------------------------------------------------------
info("Calculating horizontal interpolation parameters...")
# Get grid data
with _hdf5_lock:
ds_grid = xr.open_dataset(self.domain.dynamics_grid)
# Create a BallTree for speeding up interpolation
grid_centers_coords_rad = np.vstack([ds_grid["clat"].values, ds_grid["clon"].values]).T
tree = BallTree(grid_centers_coords_rad, metric="haversine")
# Format inputs for the tree query
obs_coords_rad = np.deg2rad(np.vstack([lats_obs, lons_obs]).T)
# Find the closest cells using the tree
query = tree.query(obs_coords_rad, k=NADJ)
distances_sorted = query[0] * R_EARTH
iadjcells_sorted = query[1]
# Store the data
ds_interp["iadjcells"][:] = iadjcells_sorted
ds_interp["iadjcells"] = ds_interp["iadjcells"].astype(int)
ds_interp["inv_dist_adjcells"][:] = 1 / distances_sorted
# ------------------------------------------------------------
# Vertical interpolation parameters
# ------------------------------------------------------------
info("Calculating vertical interpolation parameters...")
# Get relevant data from obs.csv
alt_obs = data_obs["alt"].values
ilevels_obs = data_obs["icon_ilevel"].values
# Extract height coordinate
height_dimname = da_alt.dims[0]
# If full interpolation, using ilevel to fill missing altitude
if self.full_interpolation:
info("Using full interpolation...")
mask_nan_alt = np.isnan(alt_obs)
if np.sum(mask_nan_alt):
vals_levels = ilevels_obs[mask_nan_alt] + 1
vals_levels = xr.DataArray(vals_levels.astype(int), dims="index")
vals_cells = icells_obs[mask_nan_alt]
vals_cells = xr.DataArray(vals_cells, dims="index")
# TODO: get altitude of level using an inverse-weighted average of all adjacent cells
alt_obs[mask_nan_alt] = da_alt.loc[vals_levels, vals_cells].values
# Duplicate and store the altitude of observations for the adjacent cells
mask_nan_alt = np.isnan(alt_obs)
ds_interp["alt_adjcells"][:] = np.tile(alt_obs.reshape(-1, 1), NADJ)
# If altitude data exists, find the adjacent levels and
# calculate the interpolation factors
# TODO: accelerate this part like apply_interpolation_by_chunk_full
nobs_with_alt = np.sum(~mask_nan_alt)
if nobs_with_alt:
# Finding the heights levels of each observation for timestep 0
iadjcells_obs_flatten = ds_interp["iadjcells"][~mask_nan_alt].values.flatten()
alt_adjcells = ds_interp["alt_adjcells"][~mask_nan_alt].values.flatten()
alt_model_adjcells = da_alt.sel(ncells=iadjcells_obs_flatten)
# Check if there is NaN values in output heights
nan_in_height = np.any(np.isnan(alt_model_adjcells))
negative_in_height = np.any(alt_model_adjcells < 0)
if nan_in_height:
warning("WARNING: NaNs detected in z_mc at cells" " matching the observations!")
if negative_in_height:
warning("WARNING: Negative values detected in z_mc" " at cells matching the observations!")
mask_below = alt_model_adjcells <= alt_adjcells
ilev_adjcells_below = mask_below.argmax(dim=height_dimname).values
ilev_adjcells_above = ilev_adjcells_below - 1
mask_nolevelbelow = np.sum(mask_below, axis=0) == 0
mask_nolevelabove = np.sum(mask_below, axis=0) == self.domain.nlev
ilev_adjcells_below[mask_nolevelbelow] = ilev_adjcells_above[mask_nolevelbelow] = self.domain.nlev - 1
ilev_adjcells_below[mask_nolevelabove] = ilev_adjcells_above[mask_nolevelabove] = 0
ilev_adjcells_both = np.vstack([ilev_adjcells_above, ilev_adjcells_below])
vertical_distances = np.take_along_axis((alt_model_adjcells - alt_adjcells).values, ilev_adjcells_both, axis=0)
vertical_distances = np.abs(vertical_distances)
vert_scaling_fact = vertical_distances[1] / vertical_distances.sum(axis=0)
vert_scaling_fact[mask_nolevelbelow] = 0
vert_scaling_fact[mask_nolevelabove] = 0
ds_interp["vert_scale_fact"][~mask_nan_alt] = vert_scaling_fact.reshape((nobs_with_alt, NADJ))
ds_interp["ilev_below"][~mask_nan_alt] = ilev_adjcells_below.reshape((nobs_with_alt, NADJ))
ds_interp["ilev_above"][~mask_nan_alt] = ilev_adjcells_above.reshape((nobs_with_alt, NADJ))
min_ilevel = int(ilev_adjcells_above.min())
max_ilevel = int(ilev_adjcells_below.max())
# Save the dataset
with _hdf5_lock:
ds_interp.to_netcdf(os.path.join(output_path, "interpolation/ds_data_interpolation.nc"))
# ---------------------------------------------------------
# -- Reduce the sizes of the files
# ---------------------------------------------------------
info("Reducing the sizes...")
# Get indexes to save
cells2save = np.unique(iadjcells_sorted)
levs2save = range(min_ilevel + 1, max_ilevel + 2)
# Reduce
func_arguments = zip(
list_files,
repeat(output_path),
repeat(cells2save),
repeat(levs2save),
repeat(self.domain.nlev),
repeat(self.domain.nlon),
repeat(self.timestep),
)
with mp.Pool(min(len(list_files), 36)) as pool:
pool.starmap(reduce_size, func_arguments)
# ---------------------------------------------------------
# -- Concatenate the reduced files by day
# ---------------------------------------------------------
info("Concatenating the reduced files by day...")
file_path_reduced = os.path.join(output_path, "reduced", OUTPUT_PATTERN)
list_files_reduced = sorted(glob.glob(str(file_path_reduced)))
with _hdf5_lock:
f0 = xr.open_dataset(list_files_reduced[0])
fn = xr.open_dataset(list_files_reduced[-1])
time = pd.date_range(f0.time.values[0], fn.time.values[-1], freq=self.output_resolution)
df = pd.DataFrame(list_files_reduced, columns=["file"], index=time)
list_files2concatenate = df["file"].groupby(pd.Grouper(freq="D")).apply(list).tolist()
# Process
func_arguments = zip(repeat(output_path), list_files2concatenate)
with mp.Pool(min(len(list_files2concatenate), 36)) as pool:
pool.starmap(concatenate_byday, func_arguments)
# ------------------------------------------------------------
# Concatenate the files over the full period
# ------------------------------------------------------------
info("Concatenating the files over the full period...")
file_path_day = os.path.join(output_path, "concatenated_byday", OUTPUT_PATTERN)
list_files_day = sorted(glob.glob(str(file_path_day)))
with _hdf5_lock:
ds_icon = xr.open_mfdataset(list_files_day, concat_dim="time", combine="nested")
# Extract variable information
dim_height = ds_icon["z_mc"].dims[1]
all_trcrs = [tr for tr in list(ds_icon.data_vars) if ds_icon.data_vars[tr].dims == ("time", dim_height, "ncells")]
# Apply the interpolation for every chunk
info(
f"Applying horizontal and vertical interpolation with {self.interpolation_apply_nchunks} chunks..."
)
nobs = len(ds_interp.index)
if self.interpolation_apply_nchunks > 1:
nobs_per_chunk = int(np.ceil(nobs / self.interpolation_apply_nchunks))
bounds_idx_obs = list(range(0, nobs, nobs_per_chunk)) + [nobs]
list_segments_idx_obs = [
[idx_start, idx_end] for (idx_start, idx_end) in zip(bounds_idx_obs[:-1], bounds_idx_obs[1:])
]
func_arguments = zip(list_segments_idx_obs, repeat(ds_icon), repeat(ds_interp), repeat(all_trcrs))
with mp.Pool(min(self.interpolation_apply_nchunks, 36)) as pool:
if self.full_interpolation:
list_ds_out = pool.starmap(apply_interpolation_by_chunk_full, func_arguments)
else:
list_ds_out = pool.starmap(apply_interpolation_by_chunk_reduced, func_arguments)
# Concatenate the data
info("Concatenating outputs from multiprocessing...")
ds_out_concat = xr.concat(list_ds_out, dim="obs")
else:
if self.full_interpolation:
ds_out_concat = apply_interpolation_by_chunk_full([0, nobs], ds_icon, ds_interp, all_trcrs)
else:
ds_out_concat = apply_interpolation_by_chunk_reduced([0, nobs], ds_icon, ds_interp, all_trcrs)
# Fill dataout with interpolated data
info("Filling output DataFrame with calculated data...")
dataout = pd.DataFrame(index=ds_out_concat.obs, columns=all_trcrs + ["parameter"])
dataout.loc[:, "parameter"] = ds_out_concat["parameter"].values
ds_out_concat = ds_out_concat.drop(["parameter"])
dataout.loc[:, all_trcrs] = ds_out_concat.to_array().values.T
# Save the dataset
info("Saving the DataFrame as a NetCDF file...")
filepath_dataout = os.path.join(output_path, "dataout/dataout.nc")
with _hdf5_lock:
dataout.to_xarray().to_netcdf(filepath_dataout)
info("Output processing and interpolation are complete!")
# # ---------------------------------------------------------
# # -- MAIN
# # ---------------------------------------------------------
# if __name__ == "__main__":
# start_time = tm.time()
# # PARAMETERS
# # runsubdir = sys.argv[1]
# # LOGFILE = sys.argv[2]
# # self.domain.dynamics_grid = sys.argv[3]
# # OUTPUT_DIRNAME = sys.argv[4]
# # OUTPUT_PATTERN = sys.argv[5]
# # self.output_resolution = sys.argv[6]
# # NADJ = int(sys.argv[7])
# # self.domain.nlev = int(sys.argv[8])
# # self.domain.nlon = int(sys.argv[9])
# # self.timestep = int(sys.argv[10])
# # self.full_interpolation = bool(int(sys.argv[11]))
# runsubdir = "/scratch/snx3000/jthanwer/cif_workdirs/cif_workdir_ensrf_co2_intercomp/ensemble/20190101-00:00_to_20190121-00:00/batch_sampling/obsoperator/fwd_0000/2019-01-01_00-00"
# LOGFILE = "/scratch/snx3000/jthanwer/cif_workdirs/cif_workdir_ensrf_co2_intercomp/ensemble/20190101-00:00_to_20190121-00:00/batch_sampling/obsoperator/fwd_0000/2019-01-01_00-00/outputs_log.std"
# self.domain.dynamics_grid = "/scratch/snx3000/jthanwer/cif_workdirs/cif_workdir_ensrf_co2_intercomp/ensemble/20190101-00:00_to_20190121-00:00/batch_sampling/obsoperator/fwd_0000/2019-01-01_00-00/dyn_grid.nc"
# OUTPUT_DIRNAME = 'outputs'
# OUTPUT_PATTERN = 'OUTPUT_*.nc'
# self.output_resolution = '1H'
# NADJ = 5
# self.domain.nlev = 60
# self.domain.nlon = 21344
# self.timestep = 120
# self.full_interpolation = True
# self.interpolation_apply_nchunks = 10
# # Process
# process_output(runsubdir, self.domain.dynamics_grid, OUTPUT_RESOLUTION,
# NADJ, self.domain.nlev, self.domain.nlon, self.timestep,
# self.full_interpolation, self.interpolation_apply_nchunks)
