Skip to content

Commit 6c9a231

Browse files
tobiasmarketobiasmarke
authored
Implement support for covariance data (#31)
1 parent f48eb39 commit 6c9a231

12 files changed

Lines changed: 867 additions & 146 deletions

‎mwrpy/level1/lev1_meta_nc.py‎

Lines changed: 16 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -23,16 +23,6 @@ def get_data_attributes(rpg_variables: dict, data_type: str) -> dict:
2323
from level1.lev1_meta_nc import get_data_attributes
2424
att = get_data_attributes('data','data_type')
2525
"""
26-
if data_type not in (
27-
"1B01",
28-
"1B11",
29-
"1B21",
30-
"1C01",
31-
):
32-
raise RuntimeError(
33-
["Data type " + data_type + " not supported for file writing."]
34-
)
35-
3626
if data_type in ("1B01", "1B11", "1B21"):
3727
read_att = att_reader[data_type]
3828
attributes = dict(ATTRIBUTES_COM, **read_att)
@@ -41,6 +31,10 @@ def get_data_attributes(rpg_variables: dict, data_type: str) -> dict:
4131
attributes = dict(
4232
ATTRIBUTES_COM, **ATTRIBUTES_1B01, **ATTRIBUTES_1B11, **ATTRIBUTES_1B21
4333
)
34+
else:
35+
raise RuntimeError(
36+
["Data type " + data_type + " not supported for file writing."]
37+
)
4438

4539
for key in list(rpg_variables):
4640
if key in attributes:
@@ -185,25 +179,18 @@ def get_data_attributes(rpg_variables: dict, data_type: str) -> dict:
185179
comment="0=horizon, 90=zenith",
186180
dimensions=("time",),
187181
),
188-
# "tb_accuracy": MetaData(
189-
# long_name="Total absolute calibration uncertainty of brightness temperature,\n"
190-
# "one standard deviation",
191-
# units="K",
192-
# comment="specify here source of this variable, e.g. literature value,\n"
193-
# "specified by manufacturer, result of validation effort\n"
194-
# "(updated irregularily) For RDX systems, derived from analysis\n"
195-
# "performed by Tim Hewsion (Tim J. Hewison, 2006: Profiling Temperature\n"
196-
# "and Humidity by Ground-based Microwave Radiometers, PhD Thesis,\n"
197-
# "University of Reading.) Derived from sensitivity analysis of LN2\n"
198-
# "calibration plus instrument noise levels (ACTRIS work), \n"
199-
# "currently literature values (Maschwitz et al. for HATPRO, ? for radiometrics)",
200-
# ),
201-
# "tb_cov": MetaData(
202-
# long_name="Error covariance matrix of brightness temperature channels",
203-
# units="K*K",
204-
# comment="the covariance matrix has been determined using the xxx method\n"
205-
# "from observations at a blackbody target of temperature t_amb",
206-
# ),
182+
"tb_cov_amb": MetaData(
183+
long_name="Error covariance matrix of brightness temperature channels on ambient target.",
184+
units="K*K",
185+
comment="Brightness temperature error covariance matrix determined on ambient target.",
186+
dimensions=("frequency", "frequency"),
187+
),
188+
"tb_cov_ln2": MetaData(
189+
long_name="Error covariance matrix of brightness temperature channels on LN2 target.",
190+
units="K*K",
191+
comment="Brightness temperature error covariance matrix determined during absolute calibration on LN2 target.",
192+
dimensions=("frequency", "frequency"),
193+
),
207194
"quality_flag": MetaData(
208195
long_name="Quality flag",
209196
units="1",

‎mwrpy/level1/met_quality_control.py‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -45,3 +45,5 @@ def apply_met_qc(data: dict, params: dict) -> None:
4545
threshold_low, threshold_high = params["met_thresholds"][bit]
4646
ind = (data[name][:] < threshold_low) | (data[name][:] > threshold_high)
4747
data["met_quality_flag"][ind] = setbit(data["met_quality_flag"][ind], bit)
48+
if name == "rainfall_rate" and np.all(data[name][:] * 1000.0 * 3600.0 > 60.0):
49+
data["met_quality_flag"][:] = setbit(data["met_quality_flag"][:], bit)

‎mwrpy/level1/rpg_bin.py‎

Lines changed: 255 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -24,7 +24,9 @@ def _stack_data(source: dict, target: dict, fun: Callable):
2424
value = ma.array(value)
2525
if value.ndim > 0 and name in target:
2626
if target[name].ndim == value.ndim:
27-
if (
27+
if name == "covariance_matrix":
28+
target[name] = np.array([target[name], value])
29+
elif (
2830
value.ndim > 1
2931
and value.shape[1] != target[name].shape[1]
3032
and name in ("irt", "tb")
@@ -34,6 +36,10 @@ def _stack_data(source: dict, target: dict, fun: Callable):
3436
)
3537
else:
3638
target[name] = fun((target[name], value))
39+
elif target[name].ndim != value.ndim and name == "covariance_matrix":
40+
target[name] = np.concatenate(
41+
(target[name], value[np.newaxis, :, :])
42+
)
3743
elif value.ndim > 0 and name not in target:
3844
target[name] = value
3945

@@ -72,12 +78,13 @@ def __init__(
7278
self, file_list: list[str], time_offset: datetime.timedelta | None = None
7379
):
7480
self.header, self.raw_data = stack_files(file_list)
75-
self.raw_data["time"] = utils.epoch2unix(
76-
self.raw_data["time"], self.header["_time_ref"], time_offset=time_offset
77-
)
81+
if str(file_list[0][-3:]).lower() not in ("log", "his"):
82+
self.raw_data["time"] = utils.epoch2unix(
83+
self.raw_data["time"], self.header["_time_ref"], time_offset=time_offset
84+
)
7885
self.data: dict = {}
7986
self._init_data()
80-
if str(file_list[0][-3:]).lower() != "his":
87+
if str(file_list[0][-3:]).lower() not in ("log", "his"):
8188
try:
8289
self.find_valid_times()
8390
except ValueError as err:
@@ -389,6 +396,247 @@ def read_met(file_name: str) -> tuple[dict, dict]:
389396
return header, data
390397

391398

399+
def read_log(file_name: str) -> tuple[dict, dict]:
400+
"""Reads LOG files and returns header and data as dictionary."""
401+
header: dict = {}
402+
data: dict = {}
403+
# Reads covariance matrix
404+
if "covmatrix" in file_name.lower():
405+
with open(file_name, "r", encoding="utf-8") as file:
406+
lines = file.readlines()
407+
data["cal_date"] = np.array(
408+
[
409+
datetime.datetime.strptime(
410+
str(lines[0].split(" ")[2] + lines[0].split(" ")[4]).strip(),
411+
"%d.%m.%Y%H:%M:%S",
412+
).timestamp()
413+
]
414+
)
415+
for lineno, line in enumerate(lines):
416+
if line.startswith("Frequencies [GHz]"):
417+
da = line.split(":")[1].split()
418+
header["_f"] = np.array([float(f) for f in da], dtype=float)
419+
if line.startswith("Calibrated"):
420+
da = line.split(":")[1].split()
421+
data["status"] = np.array([int(f) for f in da], dtype=int, ndmin=2)
422+
if line.startswith(("Alpha", "Gain", "Trec", "Tnoise")):
423+
da = line.split(":")[1].split()
424+
data[line.split("[")[0].strip()] = np.array(
425+
[float(f) for f in da], dtype=float, ndmin=2
426+
)
427+
if line.startswith("Covariance Matrix [K^2]"):
428+
data["covariance_matrix"] = np.ndarray(
429+
(len(header["_f"]), len(header["_f"])), np.float32
430+
)
431+
for nf in range(len(header["_f"])):
432+
line_a = lines[lineno + nf + 1].split(":")[1].split()
433+
data["covariance_matrix"][nf, :] = np.array(
434+
[float(f) for f in line_a]
435+
)
436+
# Reads absolute calibration log files
437+
elif "abscal" in file_name.lower():
438+
with open(file_name, "r", encoding="utf-8") as file:
439+
text = file.read()
440+
if (
441+
"cancelled" in text
442+
or "failed" in text
443+
or "Rec1" not in text
444+
or "Rec2" not in text
445+
or "New" not in text
446+
):
447+
return header, data
448+
with open(file_name, "r", encoding="utf-8") as file:
449+
lines = file.readlines()
450+
data["cal_date"] = np.array(
451+
[
452+
datetime.datetime.strptime(
453+
lines[1].split("_")[1] + lines[1].split("_")[2].split(".")[0],
454+
"%y%m%d%H%M%S",
455+
).timestamp()
456+
]
457+
)
458+
for lineno, line in enumerate(lines):
459+
if line.strip().startswith("Number of calibration cycles"):
460+
header["n"] = int(line.split(":")[1].strip())
461+
if line.strip().startswith("New"):
462+
ll = 1
463+
while lines[lineno + ll].strip() != "":
464+
for varo, var in enumerate(("Trec", "Tnoise", "Gain", "Alpha")):
465+
if var in data:
466+
data[var] = np.append(
467+
data[var],
468+
np.array(
469+
[
470+
float(
471+
lines[lineno + ll].split(" ")[
472+
varo * 4 + 1
473+
]
474+
)
475+
],
476+
dtype=float,
477+
ndmin=2,
478+
),
479+
axis=1,
480+
)
481+
else:
482+
data[var] = np.array(
483+
[
484+
float(
485+
lines[lineno + ll].split(" ")[varo * 4 + 1]
486+
)
487+
],
488+
dtype=float,
489+
ndmin=2,
490+
)
491+
ll += 1
492+
# Assign frequencies based on number of channels (currently not working for LHATPRO)
493+
header["_f"] = (
494+
np.array([183.91, 184.81, 185.81, 186.81, 188.31, 190.81, 90.0])
495+
if len(data["Trec"]) == 7
496+
else np.array(
497+
[
498+
22.24,
499+
23.04,
500+
23.84,
501+
25.44,
502+
26.24,
503+
27.84,
504+
31.4,
505+
51.26,
506+
52.28,
507+
53.86,
508+
54.94,
509+
56.66,
510+
57.3,
511+
58.0,
512+
]
513+
)
514+
)
515+
data["Gain"] *= 1e3 # Convert to mV/K
516+
else:
517+
raise InvalidFileError(f"LOG file type not supported")
518+
519+
return header, data
520+
521+
522+
def read_his(file_name: str) -> tuple[dict, dict]:
523+
"""This function reads RPG MWR ABSCAL.HIS binary files."""
524+
Fill_Value_Float = -999.0
525+
Fill_Value_Int = -99
526+
527+
with open(file_name, "rb") as file:
528+
code = np.fromfile(file, np.int32, 1)
529+
if code != 39583209:
530+
raise RuntimeError(["Error: CAL file code " + str(code) + " not supported"])
531+
532+
def _get_header():
533+
"""Read header info."""
534+
n = int(np.fromfile(file, np.uint32, 1)[0])
535+
time_ref = 1
536+
header_names = ["_code", "n", "_time_ref"]
537+
header_values = [code, n, time_ref]
538+
return dict(zip(header_names, header_values))
539+
540+
def _create_variables():
541+
"""Initialize data arrays."""
542+
vrs = {
543+
"len": np.ones(header["n"], np.int32) * Fill_Value_Int,
544+
"rad_id": np.ones(header["n"], np.int32) * Fill_Value_Int,
545+
"cal1_t": np.ones(header["n"], np.int32) * Fill_Value_Int,
546+
"cal2_t": np.ones(header["n"], np.int32) * Fill_Value_Int,
547+
"t1": np.ones(header["n"], np.int32) * Fill_Value_Int,
548+
"cal_date": np.ones(header["n"], np.int32) * Fill_Value_Int,
549+
"t2": np.ones(header["n"], np.int32) * Fill_Value_Int,
550+
"a_temp1": np.ones(header["n"], np.float32) * Fill_Value_Float,
551+
"a_temp2": np.ones(header["n"], np.float32) * Fill_Value_Float,
552+
"p1": np.ones(header["n"], np.float32) * Fill_Value_Float,
553+
"p2": np.ones(header["n"], np.float32) * Fill_Value_Float,
554+
"hl_temp1": np.ones(header["n"], np.float32) * Fill_Value_Float,
555+
"hl_temp2": np.ones(header["n"], np.float32) * Fill_Value_Float,
556+
"cl_temp1": np.ones(header["n"], np.float32) * Fill_Value_Float,
557+
"cl_temp2": np.ones(header["n"], np.float32) * Fill_Value_Float,
558+
"spare": np.ones([header["n"], 5], np.float32) * Fill_Value_Float,
559+
"n_ch1": np.ones(header["n"], np.int32) * Fill_Value_Int,
560+
"freq1": np.ones([header["n"], 7], np.float32) * Fill_Value_Float,
561+
"n_ch2": np.ones(header["n"], np.int32) * Fill_Value_Int,
562+
"freq2": np.ones([header["n"], 7], np.float32) * Fill_Value_Float,
563+
"cal_flag": np.ones([header["n"], 14], np.int32) * Fill_Value_Int,
564+
"Gain": np.ones([header["n"], 14], np.float32) * Fill_Value_Float,
565+
"Tnoise": np.ones([header["n"], 14], np.float32) * Fill_Value_Float,
566+
"Trec": np.ones([header["n"], 14], np.float32) * Fill_Value_Float,
567+
"Alpha": np.ones([header["n"], 14], np.float32) * Fill_Value_Float,
568+
}
569+
return vrs
570+
571+
def _get_data():
572+
"""Loop over file to read data."""
573+
data = _create_variables()
574+
for sample in range(header["n"]):
575+
data["len"][sample] = np.fromfile(file, np.int32, 1)[0]
576+
data["rad_id"][sample] = np.fromfile(file, np.int32, 1)[0]
577+
data["cal1_t"][sample] = np.fromfile(file, np.int32, 1)[0]
578+
data["cal2_t"][sample] = np.fromfile(file, np.int32, 1)[0]
579+
data["t1"][sample] = np.fromfile(file, np.int32, 1)[0]
580+
data["cal_date"][sample] = utils.epoch2unix(data["t1"][sample], 1)
581+
data["t2"][sample] = np.fromfile(file, np.int32, 1)[0]
582+
data["a_temp1"][sample] = np.fromfile(file, np.float32, 1)[0]
583+
data["a_temp2"][sample] = np.fromfile(file, np.float32, 1)[0]
584+
data["p1"][sample] = np.fromfile(file, np.float32, 1)[0]
585+
data["p2"][sample] = np.fromfile(file, np.float32, 1)[0]
586+
data["hl_temp1"][sample] = np.fromfile(file, np.float32, 1)[0]
587+
data["hl_temp2"][sample] = np.fromfile(file, np.float32, 1)[0]
588+
data["cl_temp1"][sample] = np.fromfile(file, np.float32, 1)[0]
589+
data["cl_temp2"][sample] = np.fromfile(file, np.float32, 1)[0]
590+
data["spare"][sample,] = np.fromfile(file, np.float32, 5)
591+
data["n_ch1"][sample] = np.fromfile(file, np.int32, 1)[0]
592+
data["n_ch1"][sample] = data["n_ch1"][sample]
593+
data["freq1"][sample, 0 : data["n_ch1"][sample]] = np.fromfile(
594+
file, np.float32, int(data["n_ch1"][sample])
595+
)
596+
data["n_ch2"][sample] = np.fromfile(file, np.int32, 1)[0]
597+
data["freq2"][sample, 0 : int(data["n_ch2"][sample])] = np.fromfile(
598+
file, np.float32, int(data["n_ch2"][sample])
599+
)
600+
data["cal_flag"][
601+
sample, 0 : int(data["n_ch1"][sample] + data["n_ch2"][sample])
602+
] = np.fromfile(
603+
file, np.int32, int(data["n_ch1"][sample] + data["n_ch2"][sample])
604+
)
605+
data["Gain"][
606+
sample, 0 : int(data["n_ch1"][sample] + data["n_ch2"][sample])
607+
] = (
608+
np.fromfile(
609+
file,
610+
np.float32,
611+
int(data["n_ch1"][sample] + data["n_ch2"][sample]),
612+
)
613+
* 1000.0
614+
)
615+
data["Tnoise"][
616+
sample, 0 : int(data["n_ch1"][sample] + data["n_ch2"][sample])
617+
] = np.fromfile(
618+
file, np.float32, int(data["n_ch1"][sample] + data["n_ch2"][sample])
619+
)
620+
data["Trec"][
621+
sample, 0 : int(data["n_ch1"][sample] + data["n_ch2"][sample])
622+
] = np.fromfile(
623+
file, np.float32, int(data["n_ch1"][sample] + data["n_ch2"][sample])
624+
)
625+
data["Alpha"][
626+
sample, 0 : int(data["n_ch1"][sample] + data["n_ch2"][sample])
627+
] = np.fromfile(
628+
file, np.float32, int(data["n_ch1"][sample] + data["n_ch2"][sample])
629+
)
630+
631+
file.close()
632+
return data
633+
634+
header = _get_header()
635+
data = _get_data()
636+
header["_f"] = np.concatenate((data["freq1"][0, :], data["freq2"][0, :]))
637+
return header, data
638+
639+
392640
def _read(file: BinaryIO, fields: list[Field], count: int) -> np.ndarray:
393641
arr = np.fromfile(file, np.dtype(fields), count)
394642
if (read := len(arr)) != count:
@@ -474,4 +722,6 @@ def _fix_header(header: dict) -> dict:
474722
"hkd": read_hkd,
475723
"blb": read_blb,
476724
"bls": read_bls,
725+
"log": read_log,
726+
"his": read_his,
477727
}

0 commit comments

Comments
 (0)