Merge pull request #6 from MeteoSwiss/meteoswiss-dev

Minor release Sep. 2025

Author: leric2

.github/workflows/CL_test_metoffice.yaml

@@ -0,0 +1,44 @@
+# This workflow will install Python dependencies and run tests on different python versions
+# more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions
+
+name: CI
+
+on:
+  push:
+    branches:
+      - meteoswiss-metoffice
+  pull_request:
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: ['3.7']
+    steps:
+      - uses: actions/checkout@v3
+      - name: Set up python ${{ matrix.python-version }}
+        uses: actions/setup-python@v4
+        with:
+          python-version: ${{ matrix.python-version }}
+      - name: Display Python version
+        run: python -c "import sys; print(sys.version)"
+#      - name: Cache PIP and Poetry
+#        uses: actions/cache@v2.1.6
+#        with:
+#          path: |
+#            ~/.cache/pip/
+#            ~/.cache/pypoetry/
+#          key: pip&poetry-${{ runner.os }}-${{ hashFiles('pyproject.toml') }}
+      - name: Install poetry
+        run: |
+          python -m pip install --upgrade pip
+          pip install poetry
+      - name: Install package
+        run: poetry install --only main
+      - name: Display installed packages
+        run: |
+          poetry env list
+          poetry show
+      - name: Run tests
+        run: poetry run python -m unittest discover tests/

.gitignore

@@ -5,6 +5,7 @@ __pycache__/
 *log.txt
 .idea/
 .~lock.*
+.vscode/
 
 offline/
 dist/

mwr_raw2l1/config/L1_format.yaml

@@ -35,7 +35,8 @@ variables:
     _FillValue: null
     optional: False
     attributes:
-      long_name: Time interval endpoints
+      long_name: Start and end time (UTC) of the measurement
+      units: seconds since 1970-01-01 00:00:00
 
   lat:
     name: station_latitude
@@ -72,6 +73,7 @@ variables:
       long_name: Altitude above mean sea level of measurement station
       standard_name: altitude
       units: m
+      positive: up
 
   frequency:
     name: frequency
@@ -458,7 +460,7 @@ attributes:
   license: Closed-Use Non-Commercial General Licence 1.0 (CUNCGL)
   network_name: E-PROFILE
   campaign_name: ''
-  comment: ''
+  comment: ' '
   source: Ground Based Remote Sensing
   dependencies: None
   # history is directly set in NetCDF writer

mwr_raw2l1/config/config_0-20000-0-06610_A.yaml

@@ -19,6 +19,10 @@ station_latitude: 46.81
 station_longitude: 6.94
 station_altitude: 491.
 
+# liquid cloud check parameters
+liquid_cloud_check:
+  do_check: False
+  multiplying_factor: 0.075
 
 # instrument parameters (frequency-dependent parameters go are associated to channels sorted by increasing frequency)
 # ---------------------

mwr_raw2l1/config/qc_config.yaml

@@ -13,6 +13,3 @@ check_Tb_offset: False
 delta_azi_sun: 7  # minimum accepted absolute azimuth offset between instrument pointing and sun
 delta_ele_sun: 7  # minimum accepted absolute elevation offset between instrument pointing and sun
 Tb_threshold: [2.7, 330.0]   # Threshold for min and max Tb
-
-
-# potential future extension for liquid_cloud_flag

mwr_raw2l1/main.py

@@ -58,6 +58,15 @@ def run(inst_config_file, nc_format_config_file=None, qc_config_file=None, conca
     conf_inst = get_inst_config(inst_config_file)
     conf_nc = get_nc_format_config(nc_format_config_file)
     conf_qc = get_qc_config(qc_config_file)
+    try:
+        if conf_inst['lwcl_check'] and 'do_check' in conf_inst['lwcl_check']:
+            logger.info('Liquid cloud check activated for this instrument.')
+            conf_qc['lwcl_check'] = conf_inst['lwcl_check']['do_check']
+            conf_qc['lwcl_multiplying_factor'] = conf_inst['lwcl_check']['multiplying_factor']
+    except KeyError:
+        conf_qc['lwcl_check'] = False
+        conf_qc['lwcl_multiplying_factor'] = None
+        logger.info('No liquid cloud check configured in instrument config file.')
 
     reader = get_reader(conf_inst['reader'])
     meas_constructor = get_meas_constructor(conf_inst['meas_constructor'])

mwr_raw2l1/measurement/measurement.py

@@ -5,7 +5,7 @@
 from mwr_raw2l1.log import logger
 from mwr_raw2l1.measurement.measurement_constructors import MeasurementConstructors
 from mwr_raw2l1.measurement.measurement_helpers import channels2receiver, get_receiver_vars, is_var_in_data
-from mwr_raw2l1.measurement.measurement_qc_helpers import check_rain, check_receiver_sanity, check_sun
+from mwr_raw2l1.measurement.measurement_qc_helpers import check_rain, check_receiver_sanity, check_sun, find_lwcl_from_mwr
 from mwr_raw2l1.utils.num_utils import setbit, timedelta2s, unsetbit
 
 
@@ -17,7 +17,7 @@ def run(self, conf_qc):
         Args:
             conf_qc: configuration dictionary of the quality control. For defaults use mwr_raw2l1/config/qc_config.yaml
         """
-        self.set_coords()
+        self.set_coords(primary_src='conf')
         self.set_wavelength()
         self.set_receivers()
         self.set_inst_params()
@@ -252,6 +252,11 @@ def apply_quality_control(self, conf_qc):
             qc_thresholds = qc_thresholds[:-1]
         self.data['qc_thresholds'] = qc_thresholds
 
+        # Compute the liquid cloud flag using MWRpy threshold method:
+        if conf_qc['lwcl_check']:
+            self.data = find_lwcl_from_mwr(self.data, multiplying_factor=conf_qc['lwcl_multiplying_factor'])
+            self.data['liquid_cloud_flag_status'] = xr.ones_like(self.data['liquid_cloud_flag'], dtype=np.int32)
+
     def _setbits_qc(self, bit_nb, channel, mask_fail, mask_applied=None):
         """set values for quality_flag and quality_flag status for executed checks"""
         if mask_applied is None:

mwr_raw2l1/measurement/measurement_construct_helpers.py

@@ -1,5 +1,6 @@
 import numpy as np
 import xarray as xr
+import pandas as pd
 
 from mwr_raw2l1.errors import DimensionError, MissingInputArgument, TimeMismatch
 from mwr_raw2l1.log import logger
@@ -109,7 +110,10 @@ def rpg_to_si(all_data):
         all_data['met']['windspeed'] = all_data['met']['windspeed'] / 3.6  # km/h -> m/s
     except KeyError:  # KeyError will only occur if quantity not in data, what can well happen. Do nothing in this case
         pass
-
+    try:
+        all_data['irt']['IRT'] = all_data['irt']['IRT'] + 273.15  # °C -> K
+    except KeyError:  # KeyError will only occur if quantity not in data, what can well happen. Do nothing in this case
+        pass
     return all_data
 
 
@@ -149,7 +153,7 @@ def make_dataset(data, dims, vars, vars_opt=None, multidim_vars=None, time_vecto
     if data is None or not data:
         if time_vector is None:
             raise MissingInputArgument('if data is empty or None the input argument time_vector must be specified')
-        data = {'time': time_vector}  # start overwriting empty data variable
+        data = {'time': pd.to_datetime(time_vector)}  # start overwriting empty data variable
         for dim in dims[1:]:  # assume first dimension to be 'time'
             data[dim] = np.array([missing_val])  # other dimensions all one-element
         for var in all_vars:
@@ -174,7 +178,11 @@ def make_dataset(data, dims, vars, vars_opt=None, multidim_vars=None, time_vecto
             raise DimensionError(dims, var, nd)
         spec[var] = dict(dims=dims[0:nd], data=data[var])
 
-    return xr.Dataset.from_dict(spec)
+    ds = xr.Dataset.from_dict(spec)
+    # For some reason, this does not keep the formatting of the time coordinates so we overwrite it again
+    if not isinstance(ds['time'].data[0], np.datetime64):
+        ds['time'] = spec['time']['data'].values
+    return ds
 
 
 def to_single_dataset(data_dicts, *args, **kwargs):
@@ -221,6 +229,8 @@ def merge_aux_data(mwr_data, all_data, srcs_to_ignore=None):
                 all_data[src] = all_data[src].rename(varname_map)
 
         # interp to same time grid (time grid from blb now stems from some interp) and merge into out
+        # Note that this does not do any extrapolation which leaves some values (e.g. IRT) to NaN
+        # in case of a file starting with a scan
         srcdat_interp = all_data[src].interp(time=out['time'], method='nearest')  # nearest: flags stay integer
         out = out.merge(srcdat_interp, join='left')
 
@@ -264,7 +274,9 @@ def merge_brt_blb(all_data):
                 logger.warning(
                     'Skipping {} of {} scanning observations due to identical timestamp with zenith obs for {}'.format(
                         len(duplicate_times), len(blb_ts.time), duplicate_times))
-                out = out.merge(blb_ts, join='outer', compat='override')
+                # remove duplicate times from BRT and merge
+                out = out.sel(time=~out.time.isin(duplicate_times))
+                out = out.merge(blb_ts, join='outer')
         else:
             out = scan_to_timeseries_from_aux(all_data['blb'], hkd=all_data['hkd'])
 

mwr_raw2l1/measurement/measurement_constructors.py

@@ -60,7 +60,7 @@ def from_radiometrics(cls, readin_data, conf_inst=None):
         # dimensions and variable names for usage with make_dataset
         dims = {'mwr': ['time', 'frequency'],
                 'aux': ['time']}
-        vars = {'mwr': ['Tb', 'ele', 'azi', 'quality'],
+        vars = {'mwr': ['Tb', 'ele', 'azi', 'quality', 'T_amb'],
                 'aux': ['IRT', 'p', 'T', 'RH', 'rainflag', 'quality']}
         vars_opt = {'mwr': [],
                     'aux': []}
@@ -72,6 +72,9 @@ def from_radiometrics(cls, readin_data, conf_inst=None):
         all_data['mwr']['scanflag'] = ('time', flags_here)
         data = merge_aux_data(all_data['mwr'], all_data)
 
+        # adapt the dimensions for the T_amb variable as only 1 temperature is given for 2 receivers
+        data['T_amb'] = data['T_amb'].expand_dims(dim={'receiver_nb':2}, axis=1)
+
         data['mfr'] = 'radiometrics'  # manufacturer (lowercase)
 
         return cls(data, conf_inst)

mwr_raw2l1/measurement/measurement_qc_helpers.py

@@ -1,5 +1,6 @@
 import ephem
 import numpy as np
+import xarray as xr
 
 from mwr_raw2l1.errors import UnknownManufacturer
 from mwr_raw2l1.log import logger
@@ -20,12 +21,17 @@ def check_receiver_sanity(data, channel):
     if data['mfr'] == 'attex':
         logger.info('Cannot check receiver sanity for Attex as no status variable is in data file')
         return None, False
-    elif data['mfr'] == 'radiometrics':  # quality good if quality=0
-        return flag_check(data, 'quality', 1, channel=None)
+    elif data['mfr'] == 'radiometrics':  # quality good if quality=1 (confirmed by Radiometrics to CHMI in March 2025) 
+        return flag_check(data, 'quality', 0, channel=None)
     elif data['mfr'] == 'rpg':  # quality good if channel_quality_ok=1 and alarm=0
         masks_and_checks = []  # collect all output tuples from flag_check here
         masks_and_checks.append(flag_check(data, 'channel_quality_ok', 0, channel))
         masks_and_checks.append(flag_check(data, 'alarm', 1, channel=None))
+        masks_and_checks.append(flag_check(data, 'noisediode_ok_hum', 0, channel=None))
+        masks_and_checks.append(flag_check(data, 'noisediode_ok_temp', 0, channel=None))
+        masks_and_checks.append(flag_check(data, 'Tstab_ok_hum', 0, channel=None))
+        masks_and_checks.append(flag_check(data, 'Tstab_ok_temp', 0, channel=None))
+        masks_and_checks.append(flag_check(data, 'Tstab_ok_amb', 0, channel=None))
         # TODO: could add checks for noisediode_ok_hum, noisediode_ok_temp, Tstab_ok_hum, Tstab_ok_temp, Tstab_ok_amb
         check_applied_all = [m[1] for m in masks_and_checks]
         if any(check_applied_all):
@@ -164,3 +170,73 @@ def flag_check(data, varname, value, channel=None):
     else:
         logger.info("Cannot apply check for '{}' during quality control as variable does not exist".format(varname))
         return None, False
+
+def find_lwcl_from_mwr(data, multiplying_factor=0.075):
+    """
+    This is a copy of the MWRpy function to find liquid water cloud free periods using 31.4 GHz TB variability.
+    It uses water vapor channel as proxy for a humidity dependent threshold.
+
+    For now, it works only for HATPRO instrument as this includes some empirically derived parameters. 
+
+    Refactored to work directly with xarray data instead of dict
+
+    Args:
+        data: dataset, commonly Measurement.data
+        multiplying_factor: factor to multiply the threshold with, empirically derived
+
+    Returns:
+        dataset with liquid cloud flag set
+    """
+    # Different frequencies for window and water vapor channels depending on instrument type
+    freq_win = np.where(
+        (np.isclose(data["frequency"].values, 31.4, atol=2))
+        | (np.isclose(data["frequency"].values, 190.8, atol=1))
+    )[0]
+    freq_win = np.array([freq_win[0]]) if len(freq_win) > 1 else freq_win
+    freq_wv = np.where(
+        (np.isclose(np.round(data["frequency"][:], 1), 22.2))
+        | (np.isclose(np.round(data["frequency"][:], 1), 183.9))
+    )[0]
+
+    if len(freq_win) == 1 and len(freq_wv) == 1:
+        tb = data["Tb"].isel(frequency=freq_win)
+        tb = tb.squeeze(dim='frequency', drop=True)
+        tb_zenith = tb.where(data["scanflag"] == 0, drop=True).where((data["ele"] > 89.0) & (data["ele"] < 91.0), drop=True)
+        mean_diff_t = np.nanmean(tb.time.diff(dim='time').dt.seconds)
+        number_of_samples = 180/mean_diff_t.round() if mean_diff_t < 1.8 else 180/mean_diff_t.round()
+        # tb_std = tb_df.rolling(
+        #     pd.tseries.frequencies.to_offset(offset), center=True, min_periods=50
+        # ).std()
+        tb_std = tb_zenith.rolling(
+            time=int(number_of_samples), center=True
+        ).std()
+        number_of_samples = 600/mean_diff_t.round() if mean_diff_t < 1.8 else 600/mean_diff_t.round()
+        # tb_mx = tb_std.rolling(
+        #     pd.tseries.frequencies.to_offset(offset), center=True, min_periods=100
+        # ).max()
+        tb_mx = tb_std.rolling(
+            time=int(number_of_samples), center=True
+        ).max()
+        #tb_wv = np.squeeze(ds["tb"][:, freq_wv])
+        tb_wv = data["Tb"].isel(frequency=freq_wv)
+        tb_wv = tb_wv.squeeze(dim='frequency', drop=True)
+        # In order to compute the ratio, we need to get rid of the frequency coordinates
+
+        tb_rat = tb_wv / tb
+        tb_rat = tb_rat.rolling(
+            time=int(number_of_samples), center=True
+        ).max()
+
+        threshold_rat = tb_rat * multiplying_factor
+        
+        
+        data['liquid_cloud_flag'] = xr.where(
+            tb_mx < threshold_rat,
+            0,
+            1
+        )
+    data['liquid_cloud_flag'] = xr.where((data["ele"] > 89.0) & (data["ele"] < 91.0), data['liquid_cloud_flag'], 2)
+    # also fill nans with 2
+    data['liquid_cloud_flag'] = data['liquid_cloud_flag'].fillna(2)
+
+    return data