Skip to content

Commit

Permalink
Fixing LWP offset correction and scan flagging
Browse files Browse the repository at this point in the history
  • Loading branch information
tobiasmarke committed Oct 2, 2024

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature.
1 parent 3c980bc commit bca9c98
Showing 3 changed files with 31 additions and 20 deletions.
15 changes: 8 additions & 7 deletions mwrpy/level1/quality_control.py
Original file line number Diff line number Diff line change
@@ -238,14 +238,16 @@ def spectral_consistency(
).mean()
tb_mean = tb_mean.reindex(tb_df.index, method="nearest")

fact = [5.0, 5.0] # factor for receiver retrieval uncertainty
fact = [2.0, 2.5] # factor for receiver retrieval uncertainty
# flag for individual channels based on channel retrieval uncertainty
sin_ele = np.sin(np.deg2rad(data["elevation_angle"][ele_ind]))
flag_ind[
np.where(
np.abs(tb_df["Tb"].values[:] - tb_mean["Tb"].values[:])
> ret_rm[:, ifreq]
* fact[data["receiver"][ifreq] - 1]
* (2.0 - np.sin(np.deg2rad(data["elevation_angle"][ele_ind])))
* 3.0
/ sin_ele
)[0],
ifreq,
] = True
@@ -332,14 +334,12 @@ def spectral_consistency(
np.abs(data["tb"][:, ifreq] - data["tb_spectrum"][:, ifreq])
)

th_rec = [1.5, 2.0] # threshold for receiver mean absolute difference
# receiver flag based on mean absolute difference
sin_ele = np.sin(np.deg2rad(data["elevation_angle"][:]))
for _, rec in enumerate(data["receiver_nb"]):
flag_ind[
np.ix_(
ma.mean(abs_diff[:, data["receiver"] == rec], axis=1)
> th_rec[rec - 1]
* (2.0 - np.sin(np.deg2rad(data["elevation_angle"][:]))),
ma.mean(abs_diff[:, data["receiver"] == rec], axis=1) > 1.5 / sin_ele,
data["receiver"] == rec,
)
] = True
@@ -366,4 +366,5 @@ def rm_retrieval(ele_obs: np.ndarray, coeff: dict, freq) -> np.ndarray:
ele_ret = coeff["AG"]
if ele_ret.shape == ():
ele_ret = np.array([ele_ret])
return np.array([rm_ret[freq_ind, (np.abs(ele_ret - v)).argmin()] for v in ele_obs])
# ind = np.argwhere(np.abs(ele_obs - ele_ret[:, np.newaxis]) < 0.5)[:, 0]
return np.array([rm_ret[freq_ind, np.abs(v - ele_ret) < 0.5] for v in ele_obs])
8 changes: 5 additions & 3 deletions mwrpy/level2/lwp_offset.py
Original file line number Diff line number Diff line change
@@ -46,9 +46,11 @@ def correct_lwp_offset(
pd.tseries.frequencies.to_offset("20min"), center=True, min_periods=100
).max()

lwp[(lwcl_i != 0) | (lwp > 0.06) | (lwp_max["Lwp"] > (tb_max["Tb"] * 0.0025))] = (
np.nan
)
lwp[
((lwcl_i != 0) & (lwp_max["Lwp"] > (tb_max["Tb"] * 0.0025 / 1.5)))
| (lwp > 0.06)
| (lwp_max["Lwp"] > (tb_max["Tb"] * 0.0025))
] = np.nan
lwp_df = pd.DataFrame({"Lwp": lwp}, index=ind)
lwp_offset = lwp_df.rolling(
pd.tseries.frequencies.to_offset("20min"), center=True, min_periods=100
28 changes: 18 additions & 10 deletions mwrpy/level2/write_lev2_nc.py
Original file line number Diff line number Diff line change
@@ -203,7 +203,7 @@ def get_products(
else:
rpg_dat["iwv"] = ret_product

_get_qf(rpg_dat, lev1, coeff, index, product)
_get_qf(rpg_dat, lev1, coeff, index, index_ret, product)

elif data_type in ("2P01", "2P03"):
if data_type == "2P01":
@@ -288,8 +288,6 @@ def get_products(
if product == "absolute_humidity":
tmp_dat = tmp_dat / 1000.0

_get_qf(rpg_dat, lev1, coeff, index, product)

index_ret = np.where(
np.any(
np.abs(
@@ -311,6 +309,8 @@ def get_products(
)
rpg_dat[product][index_ret, :] = tmp_dat[index_ret, :]

_get_qf(rpg_dat, lev1, coeff, index, index_ret, product)

elif data_type == "2P02":
coeff = get_mvr_coeff(site, "tpb", lev1["frequency"][:], coeff_files)
if coeff[0]["RT"] < 2:
@@ -428,7 +428,7 @@ def get_products(
+ coeff["output_offset"][:, np.newaxis]
)

_get_qf(rpg_dat, lev1, coeff, index, "temperature")
_get_qf(rpg_dat, lev1, coeff, index, np.array(range(len(index))), "temperature")

elif data_type in ("2P04", "2P07", "2P08"):
assert temp_file is not None
@@ -497,16 +497,24 @@ def get_products(


def _get_qf(
rpg_dat: dict, lev1: nc.Dataset, coeff: dict, index: np.ndarray, product: str
rpg_dat: dict,
lev1: nc.Dataset,
coeff: dict,
index: np.ndarray,
index_ret: np.ndarray,
product: str,
) -> None:
rpg_dat[product + "_quality_flag"] = ma.masked_all((len(index)), np.int32)
rpg_dat[product + "_quality_flag_status"] = ma.masked_all((len(index)), np.int32)

_, freq_ind, _ = np.intersect1d(
lev1["frequency"][:],
coeff["FR"][:],
assume_unique=False,
return_indices=True,
)
rpg_dat[product + "_quality_flag"] = np.bitwise_or.reduce(
lev1["quality_flag"][:, freq_ind][index], axis=1
rpg_dat[product + "_quality_flag"][index_ret] = np.bitwise_or.reduce(
lev1["quality_flag"][:, freq_ind][index[index_ret]], axis=1
)
seqs_all = [
(key, len(list(val))) for key, val in groupby(lev1["pointing_flag"][:] & 1 > 0)
@@ -527,9 +535,9 @@ def _get_qf(
flg = np.bitwise_or.reduce(lev1["quality_flag"][scan, freq_ind], axis=1)
rpg_dat[product + "_quality_flag"][ind] = np.bitwise_or.reduce(flg)

rpg_dat[product + "_quality_flag_status"] = lev1["quality_flag_status"][
rpg_dat[product + "_quality_flag_status"][index_ret] = lev1["quality_flag_status"][
:, freq_ind[0]
][index]
][index[index_ret]]


def _combine_lev1(
@@ -588,7 +596,7 @@ def ele_retrieval(ele_obs: np.ndarray, coeff: dict) -> np.ndarray:
ele_ret = coeff["AG"]
if ele_ret.shape == ():
ele_ret = np.array([ele_ret])
ind = np.argmin(np.abs(ele_obs - ele_ret[:, np.newaxis]), axis=0)
ind = np.argwhere(np.abs(ele_obs - ele_ret[:, np.newaxis]) < 0.5)[:, 0]
return ele_ret[ind]


0 comments on commit bca9c98

Please sign in to comment.