Skip to content
Snippets Groups Projects
Commit 2d1e5893 authored by Aurelien Chauvigne's avatar Aurelien Chauvigne
Browse files

add level 1.5 support

parent 1e6c3477
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
ptrms_lib @ 27c3fd3e
Subproject commit 7057043e2685c422441ab228f0e670036cbf970a
Subproject commit 27c3fd3e513d28dc4df5915fc657d89d867c7852
src/converter/PtrmsConverter.py
+ 127
− 70
View file @ 2d1e5893
......@@ -101,11 +101,11 @@ class PtrmsConverter(Converter.Converter):
@param year year of the eventual yearly synthesis
"""
outdir = os.path.abspath(outdir)
# level 1 only
for lvl in levels:
for lvl in [levels]:
start_time = time.time()
if __DEBUG__:
print("*" * 20 + " LEVEL %d " % lvl + "*" * 20)
print("*" * 20 + " LEVEL %s " % lvl + "*" * 20)
try:
self.process_lvl(infiles, conf_file, cache_dir, outdir, lvl, hour)
except Exception as e:
......@@ -117,7 +117,7 @@ class PtrmsConverter(Converter.Converter):
if __DEBUG__:
print(
"*" * 20
+ " LEVEL %d" % lvl
+ " LEVEL %s" % lvl
+ " : "
+ time.strftime("%H:%M:%S ", time.gmtime(time.time() - start_time))
+ "*" * 20
......@@ -204,6 +204,8 @@ class PtrmsConverter(Converter.Converter):
prog_ptrms_lib.process(indata, year, levels=[0])
)
self.var_ids_std = [var_id + "_std" for var_id in self.lev1_mz]
# group outputs
indata = {}
outdata = {}
......@@ -212,24 +214,34 @@ class PtrmsConverter(Converter.Converter):
indata[Product.START_ACQ_ID] = raw_data_filtered[Product.START_ACQ_ID]
indata[Product.END_ACQ_ID] = raw_data_filtered[Product.END_ACQ_ID]
t_start = datetime(year,
TimeConverter.doy2dt(year, min(raw_data_filtered[Product.START_ACQ_ID])).month,
TimeConverter.doy2dt(year, min(raw_data_filtered[Product.START_ACQ_ID])).day,
TimeConverter.doy2dt(year, min(raw_data_filtered[Product.START_ACQ_ID])).hour,
TimeConverter.doy2dt(year, min(raw_data_filtered[Product.START_ACQ_ID])).minute,
TimeConverter.doy2dt(year, min(raw_data_filtered[Product.START_ACQ_ID])).second)
t_end = datetime(year,
TimeConverter.doy2dt(year, max(raw_data_filtered[Product.END_ACQ_ID])).month,
TimeConverter.doy2dt(year, max(raw_data_filtered[Product.END_ACQ_ID])).day,
TimeConverter.doy2dt(year, max(raw_data_filtered[Product.END_ACQ_ID])).hour,
TimeConverter.doy2dt(year, max(raw_data_filtered[Product.END_ACQ_ID])).minute,
TimeConverter.doy2dt(year, max(raw_data_filtered[Product.END_ACQ_ID])).second)
for var_id in self.var_ids_env:
indata[var_id] = raw_data_filtered[var_id]
if lvl == 2:
if lvl == "1.5":
# - set output time axis of level 2 products
dt = 3600 # hour by hour
d_start = datetime(year,
TimeConverter.doy2dt(year, min(raw_data_filtered[Product.START_ACQ_ID])).month,
TimeConverter.doy2dt(year, min(raw_data_filtered[Product.START_ACQ_ID])).day,
TimeConverter.doy2dt(year, min(raw_data_filtered[Product.START_ACQ_ID])).hour,
0,
0)
d_end = datetime(year,
TimeConverter.doy2dt(year, max(raw_data_filtered[Product.START_ACQ_ID])).month,
TimeConverter.doy2dt(year, max(raw_data_filtered[Product.START_ACQ_ID])).day,
TimeConverter.doy2dt(year, max(raw_data_filtered[Product.START_ACQ_ID])).hour,
0,
0) + timedelta(hours=1)
d_start = t_start(minutes=0, seconds=0)
d_end = t_end(minutes=0, seconds=0) + timedelta(hours=1)
# d_end = datetime(year,
# TimeConverter.doy2dt(year, max(raw_data_filtered[Product.START_ACQ_ID])).month,
# TimeConverter.doy2dt(year, max(raw_data_filtered[Product.START_ACQ_ID])).day,
# TimeConverter.doy2dt(year, max(raw_data_filtered[Product.START_ACQ_ID])).hour,
# 0,
# 0) + timedelta(hours=1)
ts_start = calendar.timegm(d_start.timetuple())
ts_end = calendar.timegm(d_end.timetuple())
......@@ -247,18 +259,17 @@ class PtrmsConverter(Converter.Converter):
shape = v_start_time.shape
# - init other output parametres buffers
for var_id in self.var_ids_env + self.lev1_mz + self.var_ids_flag:
for var_id in self.lev1_mz + self.var_ids_flag + self.var_ids_std:
outdata[var_id] = numpy.empty(
shape, dtype=numpy.float64) + self.get_fill(var_id)
else:
raise ValueError(
"ERROR: only level 0, 1 or 2 are requested for ACTRIS")
outdata = indata
# - acquisition time
v_start_acq_in = indata[Product.START_ACQ_ID]
v_start_acq_out = outdata[Product.START_ACQ_ID]
t_start = min(v_start_acq_in)
t_start = calendar.timegm(t_start.timetuple())
# get header
fname, sout, _ = self.get_header(t_start, lvl, v_start_acq_out, year)[:3]
......@@ -277,20 +288,27 @@ class PtrmsConverter(Converter.Converter):
v_start_acq_out = outdata["start_acq"]
v_end_acq_out = outdata["end_acq"]
sz = len(v_start_acq_out)
if lvl == 2:
if lvl == "1.5":
# - averaging
for i in range(sz):
# select records in out time range
ts_min = TimeConverter.ts2dt(v_start_acq_out[i])
ts_max = TimeConverter.ts2dt(v_end_acq_out[i])
# print(ts_min, ts_max)
# input()
s = ((v_start_acq_in.index >= ts_min) & (v_start_acq_in.index <= ts_max))
if numpy.any(s):
# -> computes means and percentiles for concentration
for var_id in self.lev1_mz:
self.set_stats_mean(
indata[var_id], indata_fill[var_id], s, i, outdata[var_id], self.get_fill(var_id))
# -> stddev
_var_id = "{0:s}_{1:s}".format(var_id, "std")
self.set_stats_std(
indata[var_id], indata_fill[var_id], s, i, outdata[_var_id])
if outdata[_var_id][i] != self.get_fill(_var_id):
outdata[_var_id][i] = outdata[_var_id][i] * 2
# -> set flags
for var_id, var_id_flag in zip(self.lev1_mz,self.var_ids_flag):
data = outdata[var_id][i]
......@@ -302,11 +320,17 @@ class PtrmsConverter(Converter.Converter):
outdata[var_id_flag][i] = 0.999
# ts to doy
dt_start_date = TimeConverter.dt_to_ts(datetime(year, 1, 1))
outdata["start_acq"] = numpy.array(
[(d - dt_start_date) / 86400.0 for d in outdata["start_acq"]])
outdata["end_acq"] = numpy.array(
[(d - dt_start_date) / 86400.0 for d in outdata["end_acq"]])
if lvl in ["1.5"]:
dt_start_date = TimeConverter.dt_to_ts(datetime(year, 1, 1))
outdata["start_acq"] = numpy.array(
[(d - dt_start_date) / 86400.0 for d in outdata["start_acq"]])
outdata["end_acq"] = numpy.array(
[(d - dt_start_date) / 86400.0 for d in outdata["end_acq"]])
# convert ppb to ppt
for var_id, var_id_flag in zip(self.lev1_mz, self.var_ids_flag):
if any(outdata[var_id][outdata[var_id_flag] == 0]):
outdata[var_id] *= 1000
# --- write data to output
if __DEBUG__:
......@@ -320,23 +344,37 @@ class PtrmsConverter(Converter.Converter):
outdata["start_acq"][i], outdata["end_acq"][i]
)
# environmental variables
if lvl in ["0", "1"]:
for var_id in self.var_ids_env:
data = outdata[var_id][i]
line += "{0:<20.3f}".format(data)
qa = 0.0
if data == self.get_fill(var_id):
qa = self.var_fills_qa
else:
qa = invalid_acq_val[var_id][i]
line += "{0:<20.3f}".format(qa)
# concentrations coef
for var_id in self.var_ids_env + self.lev1_mz:
for var_id in self.lev1_mz:
conc = outdata[var_id][i]
line += "{0:<20.3f}".format(conc)
if lvl == "1.5":
line += "{0:<20.3f}".format(outdata[var_id + '_std'][i])
qa = 0.0
if conc == self.var_fills_conc:
qa = self.var_fills_qa
if lvl < 2:
qa = invalid_acq_val[var_id][i]
if var_id in self.var_ids_env:
line += "{0:<20.3f}".format(conc)
else:
line += "{0:<20.3f}".format(conc)
if lvl in ["0", "1"]:
qa = invalid_acq_val[var_id][i]
else :
qa = outdata[var_id + '_numflag'][i]
line += "{0:<20.3f}".format(qa)
sout += line.strip() + os.linesep
# write out data
......@@ -363,7 +401,7 @@ class PtrmsConverter(Converter.Converter):
"""
if var_id in [self.QA_FLAGS_ID] + self.var_ids_flag:
return 9.999
elif var_id in self.lev1_mz:
elif var_id in self.lev1_mz + self.var_ids_std:
return self.var_fills_conc
elif var_id in self.var_ids_env:
return self.var_fills_env
......@@ -420,10 +458,13 @@ class PtrmsConverter(Converter.Converter):
@return the the number of output variables
"""
nvars = 1 # end time
nvars += 1 # status
nvars += len(self.var_ids_env) * 2 # T, P and QA
if lvl == "0":
nvars += 1 # status
if lvl in ["0", "1"]:
nvars += len(self.var_ids_env) * 2 # T, P and QA
nvars += len(self.lev1_mz) * 2 # m/e and QA flag
if lvl == "1.5":
nvars += len(self.lev1_mz) # std
return nvars
def get_time_tags(self, lvl, v_start_acq, year):
......@@ -434,7 +475,7 @@ class PtrmsConverter(Converter.Converter):
"""
dt_raw_tag = self.reader.attributes["Resolution code"]
dt_tag = self.reader.attributes["Sample duration"]
if lvl in [0, 1]:
if lvl in ["0", "1"]:
if dt_tag.endswith("mn"):
dt = float(dt_tag.replace("mn", "")) / float(
24 * 60
......@@ -442,12 +483,12 @@ class PtrmsConverter(Converter.Converter):
if dt_tag.endswith("s"):
dt = float(dt_tag.replace("s", "")) / float(24 * 60 * 60)
period_tag = "1d" # NRT mode
elif lvl in [2]:
elif lvl in ["1.5"]:
dt_tag = "1h"
dt = 1 / float(24) # 1h in decimal day unit
period_tag = "1y"
else:
raise ValueError("Invalid level %d" % lvl)
raise ValueError("Invalid level %s" % lvl)
return (period_tag, dt_tag, dt, dt_raw_tag)
def get_vars_fill(self, lvl):
......@@ -474,34 +515,40 @@ class PtrmsConverter(Converter.Converter):
@return (vars_long_name, vars_tags)
"""
vars_long_name = ""
vars_long_name += (
"{0:s}, {1:s}, Location=inlet, Matrix=instrument".format("temperature", "K")
+ os.linesep
)
vars_long_name += "numflag_temperature, no unit" + os.linesep
vars_long_name += (
"{0:s}, {1:s}, Location=inlet, Matrix=instrument".format("pressure", "hPa")
+ os.linesep
)
vars_long_name += "numflag_pressure, no unit" + os.linesep
if lvl in ["0", "1"]:
vars_long_name += (
"{0:s}, {1:s}, Location=inlet, Matrix=instrument".format("temperature", "K")
+ os.linesep
)
vars_long_name += "numflag temperature, no unit" + os.linesep
vars_long_name += (
"{0:s}, {1:s}, Location=inlet, Matrix=instrument".format("pressure", "hPa")
+ os.linesep
)
vars_long_name += "numflag pressure, no unit" + os.linesep
for var_id in self.lev1_mz_name:
vars_long_name += var_id + ", ppb" + os.linesep
vars_long_name += "numflag_" + var_id + ", no unit" + os.linesep
vars_long_name += var_id + ", pmol/mol" + os.linesep
if lvl == "1.5":
vars_long_name += var_id + ", pmol/mol, Statistics=expanded uncertainty 2sigma" + os.linesep
vars_long_name += "numflag " + var_id + ", no unit" + os.linesep
vars_long_name = vars_long_name.strip()
# tags
vars_tag = ""
for var_id in self.var_ids_env_tag:
var_id = var_id.replace(" ", "_")
vars_tag += "{0:<20s}".format(var_id)
vars_tag += "{0:<20s}".format(var_id + "_numflag")
if lvl in ["0", "1"]:
for var_id in self.var_ids_env_tag:
var_id = var_id.replace(" ", "_")
vars_tag += "{0:<20s}".format(var_id)
vars_tag += "{0:<20s}".format(var_id + "_numflag")
for var_id in self.lev1_mz:
vars_tag += "{0:<20s}".format(var_id)
if lvl == "1.5":
vars_tag += "{0:<20s}".format(var_id + "_std")
vars_tag += "{0:<20s}".format(var_id + "_numflag")
vars_tag = vars_tag.strip()
return (vars_long_name, vars_tag)
......@@ -549,8 +596,8 @@ class PtrmsConverter(Converter.Converter):
# set metdata for level 1
self.component = "ion_concentrations"
self.normal_comment_line = 58
self.unit = "ppb"
self.normal_comment_line = 59
self.unit = "pmol/mol"
# - output filename
fname = self.get_outfname(start_time, prod_time, lvl, v_start_acq, year)
......@@ -561,12 +608,12 @@ class PtrmsConverter(Converter.Converter):
# set type code + dt value : TU time regulary spaced, TI : irregulary
set_type_code = "TU" # should always be
if lvl in [0, 1]:
if lvl in ["0", "1"]:
diff = numpy.diff(v_start_acq)
if (diff.size > 0) and not numpy.all(diff == diff[0]):
set_type_code = "TI"
dt = 0
elif lvl in [2]:
elif lvl in ["1.5"]:
dt = 1 / float(24) # 1h in decimal day unit
period_tag = "1y"
......@@ -585,7 +632,7 @@ class PtrmsConverter(Converter.Converter):
# standard conditions of acquisition
std_temp = self.acq_temp
std_pres = self.acq_pres
if lvl in ["1", "2"]:
if lvl in ["1", "1.5"]:
std_temp = "%.2f" % Standardizer.T_STD
std_pres = "%.2f" % Standardizer.P_STD
......@@ -661,6 +708,16 @@ class PtrmsConverter(Converter.Converter):
"\t<conf-file> full path of the configuration file" + os.linesep
)
parser.add_option(
"-l",
"--levels",
action="store",
default=1,
type="str",
dest="levels",
help="Processed levels",
)
parser.add_option(
"-p",
"--prod-id",
......@@ -708,7 +765,7 @@ def main():
)
infiles, conf, (options, arg) = PtrmsConverter.get_parser(__APP__, SRC_DIR)
prog = PtrmsConverter(options.prod_id)
prog.process(infiles, conf, options.cache_dir, options.outdir, levels=[2])
prog.process(infiles, conf, options.cache_dir, options.outdir, levels=options.levels)
if __name__ == "__main__":
......
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