diff --git a/mwrpy/level1/quality_control.py b/mwrpy/level1/quality_control.py index c222229..8934379 100644 --- a/mwrpy/level1/quality_control.py +++ b/mwrpy/level1/quality_control.py @@ -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]) diff --git a/mwrpy/level2/lwp_offset.py b/mwrpy/level2/lwp_offset.py index 1495b24..f9d2d6f 100644 --- a/mwrpy/level2/lwp_offset.py +++ b/mwrpy/level2/lwp_offset.py @@ -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 diff --git a/mwrpy/level2/write_lev2_nc.py b/mwrpy/level2/write_lev2_nc.py index 478ce8e..1da4e4b 100644 --- a/mwrpy/level2/write_lev2_nc.py +++ b/mwrpy/level2/write_lev2_nc.py @@ -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]