import copy
import os
from logging import debug
import numpy as np
import pandas as pd
from ......utils.check.errclass import CifFileNotFoundError, CifRuntimeError
[docs]
def read_binary(path, ncols):
"""Read a Fortran-order binary float file into a 2-D NumPy array.
Args:
path (str): path to the binary file.
ncols (int): number of columns in the output array.
Returns:
np.ndarray: shape ``(nrows, ncols)`` float64 array.
"""
data = np.fromfile(path, dtype="float")
data = data.reshape((-1, ncols), order="F")
return data
[docs]
def read_obs(runsubdir):
"""Read LMDZ-ACC simulated observation binary output.
Reads ``obs.bin`` from *runsubdir* (6-column Fortran-order float binary).
Returns a zero-row array when the file does not exist.
Args:
runsubdir (str): path to the period run directory.
Returns:
np.ndarray: shape ``(nobs, 6)`` observation array.
"""
obs_file = os.path.join(runsubdir, "obs.bin")
if os.path.isfile(obs_file):
obs = read_binary(obs_file, 6)
else:
obs = np.zeros((0, 6)) # Array of size 0
return obs
[docs]
def read_sim(self, data2dump, runsubdir, mode, ddi, ddf):
"""Read the mod.txt file as simulated by LMDz"""
# Read observations
obs = read_obs(runsubdir)
# Read simulation
obs_file = os.path.join(runsubdir, "obs.bin")
sim_file = os.path.join(runsubdir, "obs_out.bin")
if os.path.isfile(sim_file):
sim = read_binary(sim_file, 6)
else:
if not os.path.isfile(obs_file):
sim = np.zeros((0, 6)) # Array of size 0
else:
raise CifFileNotFoundError(
f"LMDZ did not produce a obs_out.bin file in '{runsubdir}'. "
"Please check the directory. "
"If using the option autorestart in the observation operator, "
"always remove the first non-valid sub-directory"
)
# Observations that were not extracted by LMDZ are set to NaN
sim[sim == 0] = np.nan
col_names = ("spec", "incr", "pressure", "dpressure", "hlay", "airm")
sim = pd.DataFrame(data=sim, columns=col_names)
dataout = {}
for trid in data2dump:
component, spec = trid
dataloc = data2dump[trid][ddi]
if len(dataloc) == 0:
continue
# Putting values to the local data store
ind_spec = self.chemistry.acspecies.attributes.index(
spec.replace("__sample#", "_")) + 1
mask = np.floor(obs[:, -1]) == ind_spec
sim_spec = sim.loc[mask]
# Put simulated value into correct column
# Different case if concs, or other parameters such as pressure
# Put pressure and other auxiliary data into spec column for later
# interpolation
if component != "concs":
if len(self.chunk_indexes[ddi]["concs"][spec]) > 0:
out_transforms = np.unique(
self.chunk_indexes[ddi]["concs"][spec])
valid_transforms = [
i for i in out_transforms
if (self.chunk_indexes[ddi]["concs"][spec] == i).sum()
== len(dataloc)
]
if len(valid_transforms) == 1:
sim_spec = sim_spec.loc[self.chunk_indexes[ddi]["concs"][spec]
== valid_transforms[0]]
else:
raise CifRuntimeError(
"I could not determine what chunk should be extracted. "
"This is unlikely to happen, but if does, please "
"report to developers"
)
inds = np.concatenate([
[0], np.cumsum(dataloc.loc[:, ('metadata', "dtstep")].values[:-1])])
avg_inds = pd.Series(np.nan, index=np.arange(len(sim_spec)))
avg_inds.loc[inds] = np.arange(len(dataloc))
avg_inds = avg_inds.ffill()
dataavg = sim_spec.groupby(avg_inds.values).sum()
# Put values in dataloc
dataloc.loc[:, ("maindata", "spec")] = dataavg.loc[:, "spec"].values
if mode == "tl":
dataloc.loc[:, ("maindata", "incr")
] = dataavg.loc[:, "incr"].values
# Put simulated value into correct column
# Different case if concs, or other parameters such as pressure
# Put pressure and other auxiliary data into spec column for later
# interpolation
if trid[0] != "concs":
# Column name
col = trid[0]
dataloc.loc[:, ("maindata", "spec")] = copy.deepcopy(
dataavg.loc[:, col].values)
dataloc.loc[:, ("maindata", "incr")] = 0.
dataout[trid] = copy.deepcopy(dataloc)
return dataout
