Commit 2cfc978b authored by Antoine Berchet's avatar Antoine Berchet
Browse files

Fix issue with FLEXPART flux read

parent 3c493cd8
......@@ -53,12 +53,12 @@ def read(
mesh_network, mesh_ts = np.meshgrid(network, ts)
mesh_agl, mesh_ts = np.meshgrid(agl, ts)
ds = init_empty()
ds["spec"] = data.flatten() * getattr(tracer, "numscale", 1)
ds["date"] = mesh_ts.flatten()
ds["station"] = mesh_stat.flatten()
ds["network"] = mesh_network.flatten()
ds["parameter"] = name
ds["duration"] = 1
ds["obserror"] = 0
ds[("metadata", "date")] = mesh_ts.flatten()
ds[("metadata", "station")] = mesh_stat.flatten()
ds[("metadata", "network")] = mesh_network.flatten()
ds[("metadata", "parameter")] = name
ds[("metadata", "duration")] = 1
ds[("maindata", "obserror")] = 0
ds[("maindata", "spec")] = data.flatten() * getattr(tracer, "numscale", 1)
return ds
......@@ -35,5 +35,5 @@ def fetch(ref_dir, ref_file, date_interval, target_dir, tracer=None, **kwargs):
target_file = '{}/{}'.format(target_dir,
os.path.basename(file_glob))
path.link(file_glob, target_file)
return list_files, list_dates
......@@ -5,6 +5,7 @@ import numpy as np
from ....utils.netcdf import readnc
from ....utils.dataarrays.reindex import reindex
from logging import debug
import pandas as pd
# from ....utils.dates import j2d
......@@ -43,10 +44,11 @@ def read(self,
[self.varname_flx, self.latname_flx,
self.lonname_flx, self.timename_flx])
file_times_in = [
datetime.datetime(1900, 1, 1) + datetime.timedelta(int(t))
for t in time_jd_in]
idata_in = file_times_in.index(dd) if dd in file_times_in else 0
ddi = datetime.datetime(year=dd[0].year, month=1, day=1)
ddf = datetime.datetime(year=ddi.year + 1, month=1, day=1)
file_dates = pd.date_range(ddi, ddf, freq="1MS").to_list()
idata_in = file_dates.index(dd[0]) if dd in file_dates else 0
data_in = data_in[idata_in]
# Convert julian day (since 1-1-1900) to datetime
......
......@@ -24,11 +24,26 @@ def dump_debug(transform, transf_mapper, tmp_datastore, runsubdir, ddi):
)
else:
if "spec" not in tmp_datastore["outputs"][trid]:
continue
if "spec" in tmp_datastore["outputs"][trid]:
debug_file = "{}/dataarray_debug.nc".format(debug_dir)
pathlib.Path(debug_file).unlink(missing_ok=True)
# Turn data to dataset in case it is a dictionary of dataarrays
xr.Dataset(tmp_datastore["outputs"][trid]).to_netcdf(debug_file, mode="w")
debug_file = "{}/dataarray_debug.nc".format(debug_dir)
pathlib.Path(debug_file).unlink(missing_ok=True)
else:
for d in tmp_datastore["outputs"][trid]:
if "spec" not in tmp_datastore["outputs"][trid][d]:
continue
debug_file = d.strftime(
"{}/dataarray_debug.nc_%Y%m%d%H%M".format(debug_dir))
pathlib.Path(debug_file).unlink(missing_ok=True)
# Turn data to dataset in case it is a dictionary of dataarrays
xr.Dataset(tmp_datastore["outputs"][trid][d]).to_netcdf(
debug_file, mode="w")
# Turn data to dataset in case it is a dictionary of dataarrays
xr.Dataset(tmp_datastore["outputs"][trid]).to_netcdf(debug_file, mode="w")
\ No newline at end of file
......@@ -20,17 +20,21 @@ def forward(
if type(xmod_out[("concs", transform.spec)][di]) == dict:
return
concs = xmod_in[("concs", transform.spec)][di].set_index(["date", "station"])
bkg = xmod_in[("background", transform.spec)][di].set_index(["date", "station"])
concs = xmod_in[("concs", transform.spec)][di].set_index(
[("metadata", "date"), ("metadata", "station")])
bkg = xmod_in[("background", transform.spec)][di].set_index(
[("metadata", "date"), ("metadata", "station")])
# Drop duplicate stations
bkg = bkg[~bkg.index.duplicated(keep='first')]
# Addind background
xmod_out[("concs", transform.spec)][di].loc[:, "spec"] = \
bkg.reindex(concs.index)["spec"].values + concs["spec"].values
xmod_out[("concs", transform.spec)][di].loc[:, ("maindata", "spec")] = \
bkg.reindex(concs.index)["maindata"]["spec"].values \
+ concs["maindata"]["spec"].values
if mode == "tl":
xmod_out[("concs", transform.spec)][di].loc[:, "incr"] = \
bkg.reindex(concs.index)["incr"].values + concs["incr"].values
xmod_out[("concs", transform.spec)][di].loc[:, ("maindata", "incr")] = \
bkg.reindex(concs.index)["maindata"]["incr"].values \
+ concs["maindata"]["incr"].values
\ No newline at end of file
......@@ -61,7 +61,7 @@ def forward(
# Directly return what was read for sparse data
if getattr(tracer, "sparse_data", False):
inputs["incr"] = 0.
inputs[("maindata", "incr")] = 0.
xmod[trid][ddi] = inputs
continue
......@@ -99,12 +99,6 @@ def forward(
# Deals with different resolutions
xmod_out[x] = scale2map(tmp, tracer, cdates, tracer.domain)
if trcr == "CH4":
print(__file__)
import code
code.interact(local=dict(locals(), **globals()))
# Now deals with scalars and physical variables
if getattr(tracer, "type", "scalar") == "scalar":
# Read the tracer array and apply the present control vector
......
......@@ -199,17 +199,31 @@ def dump_datastore(
# Expand lists objects to dump in netCDF
for var in ds.variables:
if type(ds[var].values[0]) == list:
if ds[var].values.dtype != "O":
continue
# Check if non NaNs values are lists
mask = [type(d) == list for d in ds[var].values]
if np.any(mask):
expanded = np.array(list(itertools.zip_longest(
*ds[var].values, fillvalue=np.nan))).T
*ds[var].values[mask], fillvalue=np.nan))).T
max_len = expanded.shape[1]
if max_len == 1:
ds[var] = ("index", [d[0] for d in ds[var].values])
ds[var] = ("index", [d[0] if m else d
for d, m in zip(ds[var].values, mask)])
else:
for ivar in range(max_len):
out_var = []
ind_mask = 0
for k, m in enumerate(mask):
out_var.append(expanded[ind_mask][ivar]
if m else ds[var].values[k])
if m:
ind_mask += 1
ds["{}_{:02d}".format(var, ivar)] = \
("index", [e[ivar] for e in expanded])
("index", out_var)
del ds[var]
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment