Masked arrays and other

mwrpy/atmos.py

@@ -91,10 +91,10 @@ def calc_p_baro(
         (~ma.getmaskarray(q).any(axis=1)) & (~ma.getmaskarray(T).any(axis=1)), 0
     ] = p[(~ma.getmaskarray(q).any(axis=1)) & (~ma.getmaskarray(T).any(axis=1))]
     for ialt in np.arange(len(z) - 1) + 1:
-        p_baro[:, ialt] = p_baro[:, ialt - 1] * np.exp(
+        p_baro[:, ialt] = p_baro[:, ialt - 1] * ma.exp(
             -scipy.constants.g
             * (z[ialt] - z[ialt - 1])
-            / (con.RS * np.mean([Tv[:, ialt], Tv[:, ialt - 1]]))
+            / (con.RS * (Tv[:, ialt] + Tv[:, ialt - 1]) / 2.0)
         )
 
     return p_baro
@@ -179,7 +179,7 @@ def find_lwcl_free(lev1: dict) -> tuple[np.ndarray, np.ndarray]:
         tb_rat = tb_rat.rolling(
             pd.tseries.frequencies.to_offset("20min"), center=True, min_periods=100
         ).max()
-        index[tb_mx["Tb"] < np.median(tb_rat["Tb"]) * 0.1] = 0
+        index[tb_mx["Tb"] < np.nanmedian(tb_rat["Tb"]) * 0.1] = 0
 
         df = pd.DataFrame({"index": index}, index=ind)
         df = df.bfill(limit=120)

mwrpy/level1/quality_control.py

@@ -12,9 +12,6 @@
 from mwrpy.level2.write_lev2_nc import retrieval_input
 from mwrpy.utils import get_coeff_list, setbit
 
-Fill_Value_Float = -999.0
-Fill_Value_Int = -99
-
 
 def apply_qc(
     site: str | None, data_in: RpgBin, params: dict, coeff_files: list | None
@@ -54,7 +51,7 @@ def apply_qc(
                 data["quality_flag_status"][:, freq], 0
             )
         else:
-            ind = np.where(data["tb"][:, freq] == Fill_Value_Float)
+            ind = np.where(ma.getmaskarray(data["tb"][:, freq]))
             data["quality_flag"][ind, freq] = setbit(data["quality_flag"][ind, freq], 0)
 
         # Bit 2: TB threshold (lower range)
@@ -126,12 +123,12 @@ def orbpos(data: dict, params: dict) -> np.ndarray:
     and returns index for observations in the direction of the sun"""
 
     sun: dict = {
-        "azimuth_angle": np.zeros(data["time"].shape) * Fill_Value_Float,
-        "elevation_angle": np.zeros(data["time"].shape) * Fill_Value_Float,
+        "azimuth_angle": ma.masked_all(data["time"].shape),
+        "elevation_angle": ma.masked_all(data["time"].shape),
     }
     moon: dict = {
-        "azimuth_angle": np.zeros(data["time"].shape) * Fill_Value_Float,
-        "elevation_angle": np.zeros(data["time"].shape) * Fill_Value_Float,
+        "azimuth_angle": ma.masked_all(data["time"].shape),
+        "elevation_angle": ma.masked_all(data["time"].shape),
     }
 
     sol = ephem.Sun()
@@ -170,23 +167,39 @@ def orbpos(data: dict, params: dict) -> np.ndarray:
 
     flag_ind = np.where(
         (
-            (data["elevation_angle"] != Fill_Value_Float)
+            (~ma.getmaskarray(data["elevation_angle"]))
             & (data["elevation_angle"] <= np.max(sun["elevation_angle"]) + 10.0)
             & (data["time"] >= sun["rise"])
             & (data["time"] <= sun["set"])
             & (data["elevation_angle"] >= sun["elevation_angle"] - params["saf"])
             & (data["elevation_angle"] <= sun["elevation_angle"] + params["saf"])
-            & (data["azimuth_angle"] >= sun["azimuth_angle"] - params["saf"])
-            & (data["azimuth_angle"] <= sun["azimuth_angle"] + params["saf"])
+            & (
+                data["azimuth_angle"]
+                >= sun["azimuth_angle"]
+                - params["saf"] / np.cos(np.deg2rad(data["elevation_angle"]))
+            )
+            & (
+                data["azimuth_angle"]
+                <= sun["azimuth_angle"]
+                + params["saf"] / np.cos(np.deg2rad(data["elevation_angle"]))
+            )
         )
         | (
             (data["elevation_angle"] <= np.max(moon["elevation_angle"]) + 10.0)
             & (data["time"] >= moon["rise"])
             & (data["time"] <= moon["set"])
             & (data["elevation_angle"] >= moon["elevation_angle"] - params["saf"])
             & (data["elevation_angle"] <= moon["elevation_angle"] + params["saf"])
-            & (data["azimuth_angle"] >= moon["azimuth_angle"] - params["saf"])
-            & (data["azimuth_angle"] <= moon["azimuth_angle"] + params["saf"])
+            & (
+                data["azimuth_angle"]
+                >= moon["azimuth_angle"]
+                - params["saf"] / np.cos(np.deg2rad(data["elevation_angle"]))
+            )
+            & (
+                data["azimuth_angle"]
+                <= moon["azimuth_angle"]
+                + params["saf"] / np.cos(np.deg2rad(data["elevation_angle"]))
+            )
         )
     )[0]
 
@@ -201,7 +214,7 @@ def spectral_consistency(
 
     flag_ind = np.zeros(data["tb"].shape, dtype=np.int32)
     abs_diff = ma.masked_all(data["tb"].shape, dtype=np.float32)
-    data["tb_spectrum"] = np.ones(data["tb"].shape) * np.nan
+    data["tb_spectrum"] = ma.masked_all(data["tb"].shape)
 
     c_list = get_coeff_list(site, "spc", coeff_files)
 

mwrpy/level1/rpg_bin.py

@@ -6,13 +6,11 @@
 from typing import Any, BinaryIO, Literal, TypeAlias
 
 import numpy as np
+from numpy import ma
 
 from mwrpy import utils
 from mwrpy.exceptions import InvalidFileError, MissingInputData
 
-Fill_Value_Float = -999.0
-Fill_Value_Int = -99
-
 Dim = int | tuple[int, ...]
 Field = tuple[str, str] | tuple[str, str, Dim]
 
@@ -22,17 +20,15 @@ def stack_files(file_list: list[str]) -> tuple[dict, dict]:
 
     def _stack_data(source: dict, target: dict, fun: Callable):
         for name, value in source.items():
+            value = ma.array(value)
             if value.ndim > 0 and name in target:
                 if target[name].ndim == value.ndim:
                     if (
                         value.ndim > 1
                         and value.shape[1] != target[name].shape[1]
                         and name == "irt"
                     ):
-                        value = np.hstack(
-                            (value, np.ones((len(value), 1)) * Fill_Value_Float)
-                        )
-                        target[name] = fun((target[name], value))
+                        raise NotImplementedError("Inconsistent number of IRT channels")
                     else:
                         target[name] = fun((target[name], value))
             elif value.ndim > 0 and name not in target:

mwrpy/level1/write_lev1_nc.py

@@ -21,8 +21,6 @@
     update_lev1_attributes,
 )
 
-Fill_Value_Float = -999.0
-Fill_Value_Int = -99
 FuncType: TypeAlias = Callable[[str], np.ndarray]
 
 
@@ -151,10 +149,10 @@ def prepare_data(
                 rpg_blb = RpgBin(file_list_blb)
                 _add_blb(rpg_bin, rpg_blb, rpg_hkd, params)
 
-        if params["azi_cor"] != Fill_Value_Float:
+        if params["azi_cor"] != -999.0:
             _azi_correction(rpg_bin.data, params)
 
-        if params["const_azi"] != Fill_Value_Float:
+        if params["const_azi"] != -999.0:
             rpg_bin.data["azimuth_angle"] = (
                 rpg_bin.data["azimuth_angle"] + params["const_azi"]
             ) % 360
@@ -164,7 +162,7 @@ def prepare_data(
                 file_list_irt = get_file_list(path_to_files, "IRT")
                 if len(file_list_irt) > 0:
                     rpg_irt = RpgBin(file_list_irt)
-                    rpg_irt.data["irt"][rpg_irt.data["irt"] <= 125.5] = Fill_Value_Float
+                    rpg_irt.data["irt"][rpg_irt.data["irt"] <= 125.5] = ma.masked
                     rpg_bin.data["ir_wavelength"] = rpg_irt.header["_f"] * 1e-6
                     rpg_bin.data["ir_bandwidth"] = params["ir_bandwidth"] * 1e-6
                     rpg_bin.data["ir_beamwidth"] = params["ir_beamwidth"]
@@ -277,7 +275,7 @@ def _append_hkd(
             "latitude",
         )
     else:
-        idx = np.where(hkd.data["latitude"] != Fill_Value_Float)[0]
+        idx = np.where(~ma.getmaskarray(hkd.data["latitude"]))[0]
         add_interpol1d(
             rpg_bin.data,
             hkd.data["latitude"][idx],
@@ -292,7 +290,7 @@ def _append_hkd(
             "longitude",
         )
     else:
-        idx = np.where(hkd.data["longitude"] != Fill_Value_Float)[0]
+        idx = np.where(~ma.getmaskarray(hkd.data["longitude"]))[0]
         add_interpol1d(
             rpg_bin.data,
             hkd.data["longitude"][idx],
@@ -301,7 +299,7 @@ def _append_hkd(
         )
 
     if data_type in ("1B01", "1C01"):
-        hkd.data["temp"][hkd.data["temp"] >= 350.0] = Fill_Value_Float
+        hkd.data["temp"][hkd.data["temp"] >= 350.0] = ma.masked
         add_interpol1d(
             rpg_bin.data, hkd.data["temp"][:, 0:2], hkd.data["time"], "t_amb"
         )
@@ -356,8 +354,8 @@ def _add_bls(brt: RpgBin, bls: RpgBin, hkd: RpgBin, params: dict) -> None:
 
     bls.data["pointing_flag"] = np.ones(len(bls.data["time"]), np.int32)
     bls.data["liquid_cloud_flag"] = np.ones(len(bls.data["time"]), np.int32) * 2
-    bls.data["liquid_cloud_flag_status"] = (
-        np.ones(len(bls.data["time"]), np.int32) * Fill_Value_Int
+    bls.data["liquid_cloud_flag_status"] = ma.masked_all(
+        len(bls.data["time"]), np.int32
     )
     brt.data["time"] = np.concatenate((brt.data["time"], bls.data["time"]))
     ind = np.argsort(brt.data["time"])
@@ -507,7 +505,7 @@ def _add_blb(brt: RpgBin, blb: RpgBin, hkd: RpgBin, params: dict) -> None:
     if len(time_add) > 0:
         pointing_flag_add = np.ones(len(time_add), np.int32)
         liquid_cloud_flag_add = np.ones(len(time_add), np.int32) * 2
-        liquid_cloud_flag_status_add = np.ones(len(time_add), np.int32) * Fill_Value_Int
+        liquid_cloud_flag_status_add = ma.masked_all(len(time_add), np.int32)
         brt.data["time"] = np.concatenate((brt.data["time"], time_add))
         ind = np.argsort(brt.data["time"])
         brt.data["time"] = brt.data["time"][ind]

mwrpy/level2/get_ret_coeff.py

@@ -1,12 +1,10 @@
 """Module to load in retrieval coefficient files"""
 import netCDF4 as nc
 import numpy as np
+from numpy import ma
 
 from mwrpy.utils import get_coeff_list
 
-Fill_Value_Float = -999.0
-Fill_Value_Int = -99
-
 
 def get_mvr_coeff(
     site: str | None, prefix: str, freq: np.ndarray, coeff_files: list | None
@@ -55,12 +53,12 @@ def get_mvr_coeff(
         coeff["FR_BL"] = coeff["FR"]
 
     elif (str(c_list[0][-2:]).lower() == "nc") and (len(c_list) > 0):
-        coeff["RT"] = Fill_Value_Int
+        coeff["RT"] = -1
         N = len(c_list)
 
         if prefix in ("lwp", "iwv"):
-            coeff["AG"] = np.ones(N) * Fill_Value_Float
-            coeff["FR"] = np.ones([len(freq), N]) * Fill_Value_Float
+            coeff["AG"] = ma.masked_all(N)
+            coeff["FR"] = ma.masked_all([len(freq), N])
             coeff["TL"] = np.zeros([N, len(freq)])
             coeff["TQ"] = np.zeros([N, len(freq)])
             coeff["OS"] = np.zeros(N)
@@ -89,8 +87,8 @@ def get_mvr_coeff(
             c_file = nc.Dataset(c_list[0])
             n_height_grid = c_file.dimensions["n_height_grid"].size
 
-            coeff["AG"] = np.ones(N) * Fill_Value_Float
-            coeff["FR"] = np.ones([len(freq), N]) * Fill_Value_Float
+            coeff["AG"] = ma.masked_all(N)
+            coeff["FR"] = ma.masked_all([len(freq), N])
             coeff["TL"] = np.zeros([N, n_height_grid, len(freq)])
             coeff["TQ"] = np.zeros([N, n_height_grid, len(freq)])
             coeff["OS"] = np.zeros([n_height_grid, N])
@@ -145,7 +143,7 @@ def get_mvr_coeff(
                 ]
             )
 
-    if (coeff["RT"] < 2) & (len(coeff["AL"]) == 1):
+    if (coeff["RT"] < 2) and (len(coeff["AL"]) == 1):
 
         def f_offset(x):
             return np.array(
@@ -160,7 +158,7 @@ def f_lin(x):
                 [coeff["TL"][(np.abs(coeff["AG"] - v)).argmin(), :] for v in x]
             )
 
-        if coeff["RT"] in (1, Fill_Value_Int):
+        if coeff["RT"] in (1, -1):
             if coeff["RT"] == 1:
                 coeff["TQ"] = coeff["TQ"][np.newaxis, :]
 
@@ -169,7 +167,7 @@ def f_quad(x):
                     [coeff["TQ"][(np.abs(coeff["AG"] - v)).argmin(), :] for v in x]
                 )
 
-    elif (coeff["RT"] < 2) & (len(coeff["AL"]) > 1) & (prefix != "tpb"):
+    elif (coeff["RT"] < 2) and (len(coeff["AL"]) > 1) and (prefix != "tpb"):
 
         def f_offset(x):
             return np.array(
@@ -184,7 +182,7 @@ def f_lin(x):
                 [coeff["TL"][(np.abs(coeff["AG"] - v)).argmin(), :, :] for v in x]
             )
 
-        if coeff["RT"] in (1, Fill_Value_Int):
+        if coeff["RT"] in (1, -1):
             if coeff["RT"] == 1:
                 coeff["TQ"] = coeff["TQ"][np.newaxis, :, :]
 
@@ -193,7 +191,7 @@ def f_quad(x):
                     [coeff["TQ"][(np.abs(coeff["AG"] - v)).argmin(), :, :] for v in x]
                 )
 
-    elif (coeff["RT"] < 2) & (len(coeff["AL"]) > 1) & (prefix == "tpb"):
+    elif (coeff["RT"] < 2) and (len(coeff["AL"]) > 1) and (prefix == "tpb"):
 
         def f_offset(_x):
             return coeff["OS"]
@@ -262,8 +260,8 @@ def factor(x):
     if str(c_list[0][-3:]).lower() == "ret":
         retrieval_type = ["linear regression", "quadratic regression", "neural network"]
         coeff["retrieval_type"] = retrieval_type[int(coeff["RT"][0])]
-        coeff["retrieval_elevation_angles"] = str(coeff["AG"])
-        coeff["retrieval_frequencies"] = str(coeff["FR"][:])
+        coeff["retrieval_elevation_angles"] = coeff["AG"]
+        coeff["retrieval_frequencies"] = coeff["FR"]
         if coeff["TS"] == 0:
             coeff["retrieval_auxiliary_input"] = "no_surface"
         else:
@@ -272,8 +270,8 @@ def factor(x):
 
     elif str(c_list[0][-2:]).lower() == "nc":
         coeff["retrieval_type"] = c_file.regression_type
-        coeff["retrieval_elevation_angles"] = str(coeff["AG"])
-        coeff["retrieval_frequencies"] = str(c_file["freq"][:])
+        coeff["retrieval_elevation_angles"] = coeff["AG"]
+        coeff["retrieval_frequencies"] = c_file["freq"]
         coeff["retrieval_auxiliary_input"] = c_file.surface_mode
         coeff["retrieval_description"] = c_file.retrieval_version