Include flag in LWP offset correction and adapt for frequent scanning

mwrpy/atmos.py

@@ -164,19 +164,21 @@ def find_lwcl_free(lev1: dict) -> tuple[np.ndarray, np.ndarray]:
         )
         ind = utils.time_to_datetime_index(lev1["time"][:])
         tb_df = pd.DataFrame({"Tb": tb}, index=ind)
+        offset = "3min" if np.nanmean(np.diff(lev1["time"])) < 1.8 else "10min"
         tb_std = tb_df.rolling(
-            pd.tseries.frequencies.to_offset("2min"), center=True, min_periods=40
+            pd.tseries.frequencies.to_offset(offset), center=True, min_periods=50
         ).std()
+        offset = "20min" if np.nanmean(np.diff(lev1["time"])) < 1.8 else "60min"
         tb_mx = tb_std.rolling(
-            pd.tseries.frequencies.to_offset("20min"), center=True, min_periods=100
+            pd.tseries.frequencies.to_offset(offset), center=True, min_periods=100
         ).max()
 
         tb_wv = np.squeeze(lev1["tb"][:, freq_22])
         tb_rat = pd.DataFrame({"Tb": tb_wv / tb}, index=ind)
         tb_rat = tb_rat.rolling(
-            pd.tseries.frequencies.to_offset("20min"), center=True, min_periods=100
+            pd.tseries.frequencies.to_offset(offset), center=True, min_periods=100
         ).max()
-        index[tb_mx["Tb"] < np.nanmedian(tb_rat["Tb"]) * 0.1] = 0
+        index[tb_mx["Tb"] < tb_rat["Tb"] * 0.075] = 0
 
         df = pd.DataFrame({"index": index}, index=ind)
         df = df.bfill(limit=120)

mwrpy/level2/lwp_offset.py

@@ -1,13 +1,18 @@
 """Module for LWP offset correction."""
 
+from itertools import groupby
+
 import numpy as np
 import pandas as pd
 
 from mwrpy import utils
 
 
 def correct_lwp_offset(
-    lev1: dict, lwp_org: np.ndarray, index: np.ndarray
+    lev1: dict,
+    lwp_org: np.ndarray,
+    index: np.ndarray,
+    qf: np.ndarray,
 ) -> tuple[np.ndarray, np.ndarray]:
     """This function corrects Lwp offset using the
     2min Lwp standard deviation and the water vapor
@@ -17,6 +22,7 @@ def correct_lwp_offset(
         lev1: Level 1 data.
         lwp_org: Lwp array.
         index: Index to use.
+        qf: Quality flag
     """
     if "elevation_angle" in lev1:
         elevation_angle = lev1["elevation_angle"][:]
@@ -28,12 +34,17 @@ def correct_lwp_offset(
     lwp = np.copy(lwp_org)
     lwp[elevation_angle[index] < 89.0] = np.nan
     lwp_df = pd.DataFrame({"Lwp": lwp}, index=ind)
+    offset = "3min" if np.nanmean(np.diff(lev1["time"])) < 1.8 else "10min"
     lwp_std = lwp_df.rolling(
-        pd.tseries.frequencies.to_offset("2min"), center=True, min_periods=40
+        pd.tseries.frequencies.to_offset(offset), center=True, min_periods=50
     ).std()
+    offset = "20min" if np.nanmean(np.diff(lev1["time"])) < 1.8 else "60min"
     lwp_max = lwp_std.rolling(
-        pd.tseries.frequencies.to_offset("20min"), center=True, min_periods=100
+        pd.tseries.frequencies.to_offset(offset), center=True, min_periods=100
     ).max()
+    lwp_min = lwp_df.rolling(
+        pd.tseries.frequencies.to_offset(offset), center=True, min_periods=100
+    ).min()
     freq_31 = np.where(np.round(lev1["frequency"][:].data, 1) == 31.4)[0]
     freq_22 = np.where(np.round(lev1["frequency"][:].data, 1) == 22.2)[0]
     tb = lev1["tb"][:][index, :]
@@ -43,33 +54,43 @@ def correct_lwp_offset(
         index=ind,
     )
     tb_max = tb_rat.rolling(
-        pd.tseries.frequencies.to_offset("20min"), center=True, min_periods=100
+        pd.tseries.frequencies.to_offset(offset), center=True, min_periods=100
     ).max()
 
     lwp[
-        ((lwcl_i != 0) & (lwp_max["Lwp"] > (tb_max["Tb"] * 0.0025 / 1.5)))
+        (
+            (lwcl_i != 0)
+            & (lwp_max["Lwp"] > (tb_max["Tb"] * 0.0025 / 1.5))
+            & (lwp_min["Lwp"] > -0.025)
+        )
         | (lwp > 0.06)
         | (lwp_max["Lwp"] > (tb_max["Tb"] * 0.0025))
     ] = np.nan
+
+    seqs_all = [(key, len(list(val))) for key, val in groupby(np.isfinite(lwp))]
+    seqs = np.array(
+        [
+            (key, sum(s[1] for s in seqs_all[:i]), length)
+            for i, (key, length) in enumerate(seqs_all)
+            if bool(key) is True
+        ]
+    )
+    if len(seqs) > 0:
+        for sec in range(len(seqs)):
+            if (
+                lev1["time"][seqs[sec, 1] + seqs[sec, 2] - 1]
+                - lev1["time"][seqs[sec, 1]]
+                < 180
+            ):
+                lwp[seqs[sec, 1] : seqs[sec, 1] + seqs[sec, 2] - 1] = 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
+        pd.tseries.frequencies.to_offset("60min"), center=True, min_periods=10
     ).mean()
     lwp_offset = lwp_offset.interpolate(method="linear")
     lwp_offset = lwp_offset.bfill()
 
-    lwp_mx = lwp_offset.rolling(
-        pd.tseries.frequencies.to_offset("60min"), center=True, min_periods=100
-    ).max()
-    lwp_mn = lwp_offset.rolling(
-        pd.tseries.frequencies.to_offset("60min"), center=True, min_periods=100
-    ).min()
-    lwp_offset.loc[
-        (lwp_mx["Lwp"][:] - lwp_mn["Lwp"][:]) / tb_max["Tb"] > (tb_max["Tb"] * 0.002),
-        "Lwp",
-    ] = np.nan
-    lwp_offset = lwp_offset.interpolate(method="linear")
-    lwp_offset = lwp_offset.bfill()
     lwp_offset.loc[(np.isnan(lwp_offset["Lwp"])) | (lwp_org == -999.0), "Lwp"] = 0.0
     lwp_org -= lwp_offset["Lwp"].values
     return lwp_org, lwp_offset["Lwp"].values

mwrpy/level2/write_lev2_nc.py

@@ -95,8 +95,8 @@ def get_products(
     else:
         elevation_angle = 90 - lev1["zenith_angle"][:]
 
-    rpg_dat, coeff, index, scan_time = (
-        {},
+    rpg_dat: dict = {}
+    coeff, index, scan_time = (
         {},
         np.empty([0], np.int32),
         np.empty([0], np.int32),
@@ -206,6 +206,7 @@ def get_products(
         if product in ("lwp", "iwv"):
             ret_product = ma.masked_all(len(index), np.float32)
             ret_product[index_ret] = tmp_product[index_ret]
+            _get_qf(rpg_dat, lev1, coeff, index, index_ret, product)
             if product == "lwp":
                 freq_win = np.where((np.round(lev1["frequency"][:].data, 1) == 31.4))[0]
                 rpg_dat["lwp"], rpg_dat["lwp_offset"] = (
@@ -217,13 +218,14 @@ def get_products(
                         rpg_dat["lwp"][index_ret],
                         rpg_dat["lwp_offset"][index_ret],
                     ) = correct_lwp_offset(
-                        lev1.variables, ret_product[index_ret], index[index_ret]
+                        lev1.variables,
+                        ret_product[index_ret],
+                        index[index_ret],
+                        rpg_dat["lwp_quality_flag"][index_ret],
                     )
             else:
                 rpg_dat[product] = ret_product
 
-            _get_qf(rpg_dat, lev1, coeff, index, index_ret, product)
-
         else:
             product_list = (
                 "k_index",