ENH: keep track of nsamples in timeavg flag

Author: wolfidan

src/pyrad_proc/pyrad/proc/process_intercomp.py

@@ -1053,7 +1053,13 @@ def process_weighted_time_avg(procstatus, dscfg, radar_list=None):
 def process_time_avg_flag(procstatus, dscfg, radar_list=None):
     """
     computes a flag field describing the conditions of the data used while
-    averaging
+    averaging. The flag is an integer that tracks up to 999 occurrences of the number of samples as well as
+    of three conditions during data accumulation: 𝜙𝑑𝑝 exceeding a threshold (𝜙𝑑𝑝ₘₐₓ), clutter,
+    and non-rain precipitation. It is a packed representation of 4 numbers into one large integer.
+    Flags are encoded by adding +1 for nsamples, +1E3 for 𝜙𝑑𝑝ₘₐₓ exceedance, + 1E6 for clutter,
+    and +1E9 for non-rain.
+    Inputs include the 𝜙𝑑𝑝 field, a dynamic clutter map, and either a hydrometeor
+    classification or temperature field to identify precipitation phase.
 
     Parameters
     ----------
@@ -1126,25 +1132,27 @@ def process_time_avg_flag(procstatus, dscfg, radar_list=None):
         time_avg_flag = pyart.config.get_metadata("time_avg_flag")
         time_avg_flag["data"] = np.ma.zeros((radar.nrays, radar.ngates), dtype=int)
 
+        time_avg_flag["data"] += 1  # number of samples
+
         if phidp_name not in radar.fields:
             warn("Missing PhiDP data")
-            time_avg_flag["data"] += 1
+            time_avg_flag["data"] += int(1e3)
         else:
             phidp_field = radar.fields[phidp_name]
-            time_avg_flag["data"][phidp_field["data"] > phidpmax] += 1
+            time_avg_flag["data"][phidp_field["data"] > phidpmax] += int(1e3)
 
         if echo_name is not None:
             if echo_name not in radar.fields:
                 warn("Missing echo ID data")
-                time_avg_flag["data"] += 100
+                time_avg_flag["data"] += int(1e6)
             else:
                 echo_field = radar.fields[echo_name]
-                time_avg_flag["data"][echo_field["data"] == 2] += 100
+                time_avg_flag["data"][echo_field["data"] == 2] += int(1e6)
 
         if hydro_name is not None and echo_name is not None:
             if (hydro_name not in radar.fields) or (echo_name not in radar.fields):
                 warn("Missing hydrometeor classification data")
-                time_avg_flag["data"] += 10000
+                time_avg_flag["data"] += int(1e9)
             else:
                 hydro_field = radar.fields[hydro_name]
                 # check where is no rain
@@ -1153,11 +1161,11 @@ def process_time_avg_flag(procstatus, dscfg, radar_list=None):
                 )
                 # where is no rain should be precip
                 is_not_rain = np.logical_and(is_not_rain, echo_field["data"] == 3)
-                time_avg_flag["data"][is_not_rain] += 10000
+                time_avg_flag["data"][is_not_rain] += int(1e9)
         elif temp_name is not None:
             if temp_name not in radar.fields:
                 warn("Missing temperature data")
-                time_avg_flag["data"] += 10000
+                time_avg_flag["data"] += int(1e9)
             else:
                 beamwidth = dscfg.get("beamwidth", None)
                 if beamwidth is None:
@@ -1185,11 +1193,11 @@ def process_time_avg_flag(procstatus, dscfg, radar_list=None):
                     iso0_field=iso0_name,
                     temp_ref="temperature",
                 )
-                time_avg_flag["data"][mask_fzl] += 10000
+                time_avg_flag["data"][mask_fzl] += int(1e6)
         elif iso0_name is not None:
             if iso0_name not in radar.fields:
                 warn("Missing height relative to iso0 data")
-                time_avg_flag["data"] += 10000
+                time_avg_flag["data"] += int(1e6)
             else:
                 beamwidth = dscfg.get("beamwidth", None)
                 if beamwidth is None:
@@ -1217,7 +1225,7 @@ def process_time_avg_flag(procstatus, dscfg, radar_list=None):
                     iso0_field=iso0_name,
                     temp_ref="height_over_iso0",
                 )
-                time_avg_flag["data"][mask_fzl] += 10000
+                time_avg_flag["data"][mask_fzl] += int(1e9)
 
         radar_aux = deepcopy(radar)
         radar_aux.fields = dict()
@@ -1882,7 +1890,6 @@ def process_intercomp(procstatus, dscfg, radar_list=None):
 
         fname = savedir + fname[0]
         coloc_data = read_colocated_data(fname)
-
         intercomp_dict = {
             "rad1_name": dscfg["global_data"]["rad1_name"],
             "rad1_time": coloc_data[0],
@@ -1947,13 +1954,13 @@ def process_intercomp_time_avg(procstatus, dscfg, radar_list=None):
         rng_tol : float. Dataset keyword
             range tolerance between the two radars. Default 50 m
         clt_max : int. Dataset keyword
-            maximum number of samples that can be clutter contaminated.
-            Default 100 i.e. all
+            maximum percentage of samples that can be clutter contaminated.
+            Default 100%& i.e. all
         phi_excess_max : int. Dataset keyword
-            maximum number of samples that can have excess instantaneous
-            PhiDP. Default 100 i.e. all
+            maximum percentage of samples that can have excess instantaneous
+            PhiDP. Default 100% i.e. all
         non_rain_max : int. Dataset keyword
-            maximum number of samples that can be no rain. Default 100 i.e. all
+            maximum percentage of samples that can be no rain. Default 100% i.e. all
         phi_avg_max : float. Dataset keyword
             maximum average PhiDP allowed. Default 600 deg i.e. any
 
@@ -2212,16 +2219,15 @@ def process_intercomp_time_avg(procstatus, dscfg, radar_list=None):
 
                 rad1_flag = flag1[rad1_ray_ind[i], ind_rng]
 
-                rad1_excess_phi = rad1_flag % 100
-                rad1_clt = ((rad1_flag - rad1_excess_phi) % 10000) / 100
-                rad1_prec = (
-                    (rad1_flag - rad1_clt * 100 - rad1_excess_phi) % 1000000
-                ) / 10000
-
-                flag1_vec[i] = int(
-                    10000 * np.max(rad1_prec)
-                    + 100 * np.max(rad1_clt)
-                    + np.max(rad1_excess_phi)
+                rad1_nsamples = rad1_flag % int(1e3)
+                rad1_excess_phi = (rad1_flag // int(1e3)) % int(1e3)
+                rad1_clt = (rad1_flag // int(1e6)) % int(1e3)
+                rad1_prec = (rad1_flag // int(1e9)) % int(1e3)
+                flag1_vec[i] = (
+                    int(1e9) * np.max(rad1_prec)
+                    + int(1e6) * np.max(rad1_clt)
+                    + int(1e3) * np.max(rad1_excess_phi)
+                    + np.max(rad1_nsamples)
                 )
                 is_valid_avg[i] = True
 
@@ -2270,16 +2276,16 @@ def process_intercomp_time_avg(procstatus, dscfg, radar_list=None):
 
                 rad2_flag = flag2[rad2_ray_ind[i], ind_rng]
 
-                rad2_excess_phi = rad2_flag % 100
-                rad2_clt = ((rad2_flag - rad2_excess_phi) % 10000) / 100
-                rad2_prec = (
-                    (rad2_flag - rad2_clt * 100 - rad2_excess_phi) % 1000000
-                ) / 10000
+                rad2_nsamples = rad1_flag % int(1e3)
+                rad2_excess_phi = (rad1_flag // int(1e3)) % int(1e3)
+                rad2_clt = (rad1_flag // int(1e6)) % int(1e3)
+                rad2_prec = (rad1_flag // int(1e9)) % int(1e3)
 
                 flag2_vec[i] = int(
-                    10000 * np.max(rad2_prec)
-                    + 100 * np.max(rad2_clt)
-                    + np.max(rad2_excess_phi)
+                    int(1e9) * np.max(rad2_prec)
+                    + int(1e6) * np.max(rad2_clt)
+                    + int(1e3) * np.max(rad2_excess_phi)
+                    + rad2_nsamples
                 )
                 is_valid_avg[i] = True
 
@@ -2396,17 +2402,15 @@ def process_intercomp_time_avg(procstatus, dscfg, radar_list=None):
         if rad1_time is None:
             return None, None
 
-        rad1_excess_phi = (rad1_flag % 100).astype(int)
-        rad2_excess_phi = (rad2_flag % 100).astype(int)
-        rad1_clt = (((rad1_flag - rad1_excess_phi) % 10000) / 100).astype(int)
-        rad2_clt = (((rad2_flag - rad2_excess_phi) % 10000) / 100).astype(int)
+        rad1_nsamples = rad1_flag % int(1e3)
+        rad1_excess_phi = (rad1_flag // int(1e3)) % int(1e3)
+        rad1_clt = (rad1_flag // int(1e6)) % int(1e3)
+        rad1_non_rain = (rad1_flag // int(1e9)) % int(1e3)
 
-        rad1_non_rain = (
-            ((rad1_flag - rad1_clt * 100 - rad1_excess_phi) % 1000000) / 10000
-        ).astype(int)
-        rad2_non_rain = (
-            ((rad2_flag - rad2_clt * 100 - rad2_excess_phi) % 1000000) / 10000
-        ).astype(int)
+        rad2_nsamples = rad2_flag % int(1e3)
+        rad2_excess_phi = (rad2_flag // int(1e3)) % int(1e3)
+        rad2_clt = (rad2_flag // int(1e6)) % int(1e3)
+        rad2_non_rain = (rad2_flag // int(1e9)) % int(1e3)
 
         clt_max = dscfg.get("clt_max", 100)
         phi_excess_max = dscfg.get("phi_excess_max", 100)
@@ -2417,14 +2421,14 @@ def process_intercomp_time_avg(procstatus, dscfg, radar_list=None):
         ind_val = np.where(
             np.logical_and.reduce(
                 (
-                    rad1_clt <= clt_max,
-                    rad2_clt <= clt_max,
-                    rad1_excess_phi <= phi_excess_max,
-                    rad2_excess_phi <= phi_excess_max,
-                    rad1_non_rain <= non_rain_max,
-                    rad2_non_rain <= non_rain_max,
-                    rad1_phi <= phi_avg_max,
-                    rad2_phi <= phi_avg_max,
+                    rad1_clt / rad2_nsamples <= clt_max,
+                    rad2_clt / rad2_nsamples <= clt_max,
+                    rad1_excess_phi / rad2_nsamples <= phi_excess_max,
+                    rad2_excess_phi / rad2_nsamples <= phi_excess_max,
+                    rad1_non_rain / rad1_nsamples <= non_rain_max,
+                    rad2_non_rain / rad1_nsamples <= non_rain_max,
+                    rad1_phi / rad1_nsamples <= phi_avg_max,
+                    rad2_phi / rad1_nsamples <= phi_avg_max,
                 )
             )
         )[0]

src/pyrad_proc/pyrad/prod/process_intercomp_products.py

@@ -153,7 +153,6 @@ def generate_intercomp_products(dataset, prdcfg):
         )
 
         fname = savedir + fname[0]
-
         write_colocated_data(dataset["intercomp_dict"], fname)
 
         print("saved colocated data file: " + fname)
@@ -194,6 +193,11 @@ def generate_intercomp_products(dataset, prdcfg):
         timeformat = "%Y%m%d"
         if scatter_type == "instant":
             timeformat = "%Y%m%d%H%M%S"
+            colname = "dBZavg"
+        else:
+            colname = "val"
+        if f"rad1_{colname}" not in dataset["intercomp_dict"]:
+            return
 
         field_name = get_fieldname_pyart(prdcfg["voltype"])
         savedir = get_save_dir(
@@ -247,8 +251,8 @@ def generate_intercomp_products(dataset, prdcfg):
                 continue
 
             hist_2d, bin_edges1, bin_edges2, stats = compute_2d_stats(
-                np.ma.asarray(dataset["intercomp_dict"]["rad1_val"][selection]),
-                np.ma.asarray(dataset["intercomp_dict"]["rad2_val"][selection]),
+                np.ma.asarray(dataset["intercomp_dict"][f"rad1_{colname}"][selection]),
+                np.ma.asarray(dataset["intercomp_dict"][f"rad2_{colname}"][selection]),
                 field_name,
                 field_name,
                 step1=step,

src/pyrad_proc/pyrad/util/__init__.py

@@ -98,5 +98,7 @@
 from .data_retrieval_utils import retrieve_CPCCV  # noqa
 from .data_retrieval_utils import retrieve_AQC_XLS  # noqa
 
+from .math_utils import bit_pack  # noqa
+from .math_utils import bit_unpack  # noqa
 
 __all__ = [s for s in dir() if not s.startswith("_")]

src/pyrad_proc/pyrad/util/math_utils.py

@@ -0,0 +1,71 @@
+import math
+
+
+def bit_pack(values, max_value, num_values):
+    """
+    Pack a list of integers into a single integer using bit manipulation.
+
+    Parameters
+    ----------
+    values : list of int
+        The values to pack. Each value must be <= `max_value`.
+        If fewer than `num_values` are given, values are right-padded with 0s.
+        If more than `num_values` are given, a ValueError is raised.
+    max_value : int
+        The maximum possible value any item in `values` can take.
+        Determines how many bits are allocated per value.
+    num_values : int
+        The exact number of values to pack. Controls total bit length.
+
+    Returns
+    -------
+    int
+        A single integer encoding all values.
+
+    Raises
+    ------
+    ValueError
+        If any value in `values` exceeds `max_value`, or if more than `num_values`
+        are provided.
+
+    """
+    bits = math.ceil(math.log2(max_value + 1))
+    if len(values) > num_values:
+        raise ValueError("Too many values provided")
+    padded_values = values + [0] * (num_values - len(values))
+
+    packed = 0
+    for v in padded_values:
+        if v > max_value:
+            raise ValueError(f"Value {v} exceeds max allowed {max_value}")
+        packed = (packed << bits) | v
+    return packed
+
+
+def bit_unpack(packed, max_value, num_values):
+    """
+    Unpack a single integer into a list of integers using bit manipulation.
+
+    Parameters
+    ----------
+    packed : int
+        The packed integer to decode.
+    max_value : int
+        The maximum possible value any original item could have taken.
+        Determines how many bits were allocated per value.
+    num_values : int
+        The number of values to extract.
+
+    Returns
+    -------
+    list of int
+        The unpacked values in their original order.
+
+    """
+    bits = math.ceil(math.log2(max_value + 1))
+    mask = (1 << bits) - 1
+    values = []
+    for _ in range(num_values):
+        values.append(packed & mask)
+        packed >>= bits
+    return list(reversed(values))