Relax obs requirements (#21)

* Up to 4 stations are allowed to be missing

* max_miss_stns configurable

* Update README.md

* update utils.py

* update utils.py

* small reformatting

---------

Co-authored-by: Andreas Pauling <paa@balfrin-ln003.cscs.ch>

Author: andreaspauling

README.md

@@ -110,6 +110,7 @@ const_file : <path>/CLON_CLAT_ICON-CH1.gb2
 station_obs_file : <path>/pollen_measured_2024020118.atab
 station_mod_file : <path>/pollen_modelled_2024020118.atab
 hour_incr : 1
+max_miss_stns : 4
 ```
 `pov_infile`: This GRIB2 file must include the fields `tthrs`, `tthre` (for POAC, `saisl` instead), `saisn` and `ctsum` if the module `update_phenology` is called. If the module `update_strength` is called `pov_infile` must include the fields `saisn` and `tune`. If at least one of these mandatory fields is missing the package exits with status 1 and tells the user. `pov_infile` is used as template for `pov_outfile`, i.e. the whole file is copied to `pov_outfile` with adapted values. Date and time information of `pov_infile` does not have to be correct, ICON just throws warnings.
 
@@ -119,6 +120,10 @@ hour_incr : 1
 `station_obs_file`: Observed hourly pollen concentrations (ATAB format) of the latest 120 hours relative to the target date of `pov_outfile`. The timestamps of the data in this file may vary depending on data availability, time of extraction etc. Missing values are allowed but at least 50% of each station must be there. If not, the package exits with status 1 and tells the user.
 `station_mod_file`: Modelled hourly pollen concentrations (ATAB format) of the latest 120 hours relative to the target date of `pov_outfile`. The timestamps of the data in this file may vary depending on data availability, time of extraction etc. In case of missing values the package exits with status 1 and tells the user. Same stations as in `station_obs_file` (only used if the module `update_strength` is called).
 `hour_incr`: Increment of the timestamp of the outfile relative to the infile in hours (defaults to 1; negative values also supported). This parameter should be adapted if the calibration is done for a subsequent run more than one hour ahead.
+`max_miss_stns`: Maximum number of stations with more than 50% missing observations. If
+there are less than 50% missing observation per station, they are replaced by the mean
+of the remaining observations; if there are more than 50% missing observations per station,
+the station is removed from pollen calibration. If there are more stations removed than mas_miss_stns, the pollen calibration is stopped and an alert is issued.
 
 
 ### How to run the package

pyproject.toml

@@ -7,7 +7,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "realtime-pollen-calibration"
-version = "1.1.0"
+version = "1.2.0"
 description = "Simon Adamov's Realtime Pollen Calibration, ICON extension by Andreas Pauling"
 readme = "README.md"
 keywords = [

src/realtime_pollen_calibration/set_up.py

@@ -33,6 +33,8 @@ def set_up_config(config_file: str):
 
     config.pov_outfile = data["pov_outfile"]
 
+    config.max_miss_stns = data["max_miss_stns"]
+
     config.hour_incr = data["hour_incr"]
 
     return config

src/realtime_pollen_calibration/update_phenology.py

@@ -133,7 +133,10 @@ def update_phenology_realtime(config_obj: utils.Config, verbose: bool = True):
     dict_fields = {}
     for pollen_type in utils.get_pollen_type(ds):
         obs_mod_data = utils.read_atab(
-            pollen_type, config_obj.station_obs_file, verbose=verbose
+            pollen_type,
+            config_obj.max_miss_stns,
+            config_obj.station_obs_file,
+            verbose=verbose,
         )
         change_phenology_fields = utils.get_change_phenol(
             pollen_type, obs_mod_data, ds, verbose

src/realtime_pollen_calibration/update_strength.py

@@ -102,6 +102,7 @@ def update_strength_realtime(config_obj: utils.Config, verbose: bool = True):
     for pollen_type in ptype_present:
         obs_mod_data = utils.read_atab(
             pollen_type,
+            config_obj.max_miss_stns,
             config_obj.station_obs_file,
             config_obj.station_mod_file,
             verbose=verbose,

src/realtime_pollen_calibration/utils.py

@@ -19,7 +19,7 @@
 
 
 @dataclass
-class Config:
+class Config:  # pylint: disable=too-many-instance-attributes
 
     pov_infile: str = ""
     """ICON GRIB2 file including path containing the pollen fields:
@@ -31,7 +31,7 @@ class Config:
     """ICON GRIB2 file including path of the desired output file."""
 
     t2m_file: str = ""
-    """"ICON GRIB2 file including path containing T_2M."""
+    """ICON GRIB2 file including path containing T_2M."""
 
     const_file: str = ""
     """ICON GRIB2 file including path containing Longitudes (clon)
@@ -48,7 +48,14 @@ class Config:
        pollen concentrations at the stations.
     """
 
+    max_miss_stns: int = 4
+    """Maximum number of stations with more than 50% missing data
+       allowed. If more than this number of stations are missing,
+       the pollen calibration is not performed.
+    """
+
     hour_incr: int = 1
+    """Hour increment (grib coding) between input and output POV file."""
 
 
 ObsModData = namedtuple(
@@ -125,7 +132,11 @@ def read_clon_clat(const_file):
 
 
 def read_atab(
-    pollen_type: str, file_obs_stns: str, file_mod_stns: str = "", verbose: bool = True
+    pollen_type: str,
+    max_miss_stns: int,
+    file_obs_stns: str,
+    file_mod_stns: str = "",
+    verbose: bool = True,
 ) -> ObsModData:
     # pylint: disable=too-many-locals
     """Read the pollen concentrations and the station locations from the ATAB files.
@@ -134,6 +145,7 @@ def read_atab(
         pollen_type: String describing the pollen type analysed.
         file_obs_stns: Location of the observation ATAB file.
         file_mod_stns: Location of the model ATAB file. (Optional)
+        max_miss_stns: Max. number of stations with more than 50% missing data
         verbose: Optional additional debug prints.
 
     Returns:
@@ -188,13 +200,13 @@ def read_obs_header(file_data: str):
         return HeaderData(coord_stns, missing_value, stn_indicators, n_header)
 
     headerdata = read_obs_header(file_obs_stns)
-    data = pd.read_csv(
+    data_obs = pd.read_csv(
         file_obs_stns,
         header=headerdata.n_header,
         sep=r"\s+",
         parse_dates=[[1, 2, 3, 4, 5]],
     )
-    data = data[data["PARAMETER"] == pollen_type].iloc[:, 2:].to_numpy()
+    data_obs = data_obs[data_obs["PARAMETER"] == pollen_type].iloc[:, 2:].to_numpy()
     if file_mod_stns != "":
         stn_indicators_mod = None
         missing_value = None
@@ -207,7 +219,6 @@ def read_obs_header(file_data: str):
                 if line.strip()[0:9] == "PARAMETER":
                     n_header_mod = n
                     break
-        istation_mod = get_mod_stn_index(headerdata.stn_indicators, stn_indicators_mod)
         data_mod = pd.read_csv(
             file_mod_stns,
             header=n_header_mod,
@@ -228,11 +239,24 @@ def read_obs_header(file_data: str):
     else:
         data_mod = 0
         istation_mod = 0
-    data = treat_missing(
-        data, headerdata.missing_value, headerdata.stn_indicators, verbose=verbose
+    data_obs, headerdata = treat_missing(
+        data_obs,
+        headerdata,
+        max_miss_stns,
+        headerdata.stn_indicators,
+        headerdata.missing_value,
+        verbose=verbose,
     )
+    # Calculating the station correspondence indices of obs/mod data.
+    if file_mod_stns != "":
+        istation_mod = get_mod_stn_index(headerdata.stn_indicators, stn_indicators_mod)
+
     return ObsModData(
-        data, headerdata.coord_stns, headerdata.missing_value, data_mod, istation_mod
+        data_obs,
+        headerdata.coord_stns,
+        headerdata.missing_value,
+        data_mod,
+        istation_mod,
     )
 
 
@@ -252,27 +276,36 @@ def create_data_arrays(cal_fields, clon, clat, time_values):
     return cal_fields_arrays
 
 
-def treat_missing(
+def treat_missing(  # pylint: disable=too-many-positional-arguments,too-many-arguments
     array,
+    headerdata,
+    max_miss_stns,
+    stn_indicators: np.ndarray,
     missing_value: float = -9999.0,
-    stn_indicators: str = "",
     verbose: bool = True,
 ):
     """Treat the missing values of the input array.
 
     Args:
         array: Array containing the concentration values.
+        headerdata: HeaderData containing ATAB header info.
         missing_value: Value considered as a missing value.
+        stn_indicators: Array of station indicators (PBS, PBU, ...).
+        max_miss_stns: Max. number of stations with more than 50% missing data
         verbose: Optional additional debug prints.
 
     Returns:
-        array: Modified array with removed missing values.
+        array: Modified array with replaced missing values or removed stations.
+        headerdata: Updated metadata on the stations matching the changes in array.
 
     """
     array_missing = array == missing_value
     nstns = array.shape[1]
     skip_missing_stn = np.zeros(nstns)
-    for istation in range(nstns):
+
+    stns_missing = 0
+    istation = 0
+    while istation < nstns:
         skip_missing_stn[istation] = np.count_nonzero(array_missing[:, istation])
         if verbose:
             print(
@@ -295,14 +328,48 @@ def treat_missing(
                 array[idx2, istation] = np.mean(array[idx1, istation])
             else:
                 print(
-                    f"Station {stn_indicators[istation]} has more than 50% missing ",
-                    "data. Please check the reason (Does jretrievedwh still work?).\n",
-                    "No pollen calibration update is performed until this is fixed! ",
-                    "Pollen in ICON will still work, but calibration fields get ",
-                    "more and more outdated.",
+                    f"Station {stn_indicators[istation]} has more than 50% missing\n",
+                    "data and is REMOVED from this pollen calibration run.\n",
+                    "Please check the reason (Does jretrievedwh still work?).\n",
+                    "If this occurs only at few stations the impact on the pollen\n",
+                    "calibration is minimal. However, if more stations are affected,\n",
+                    "switching off the pollen calibration should be considered,\n",
+                    "(depending on which stations, species and time of the year).",
                 )
-                sys.exit(1)
-    return array
+                stns_missing += 1
+                if stns_missing > max_miss_stns:
+                    print(
+                        f"\nALERT: More than {max_miss_stns} stations have more than\n",
+                        "50% missing data, no pollen calibration is performed!\n",
+                        "Pollen are still running but fix this asap by checking\n",
+                        "the reason for the missing observations.\n",
+                    )
+                    sys.exit(1)
+
+                # Remove the station from array, coord_stns, and stn_indicators
+                array = np.delete(array, istation, axis=1)
+
+                updated_coord_stns = (
+                    headerdata.coord_stns[:istation]
+                    + headerdata.coord_stns[istation + 1 :]
+                )
+                headerdata = HeaderData(
+                    coord_stns=updated_coord_stns,
+                    missing_value=headerdata.missing_value,
+                    stn_indicators=np.delete(headerdata.stn_indicators, istation),
+                    n_header=headerdata.n_header,
+                )
+                stn_indicators = np.delete(stn_indicators, istation)
+                array_missing = np.delete(array_missing, istation, axis=1)
+
+                # Update nstns since a station was removed
+                nstns -= 1
+                # Do not increment istation; check the next station at the same index
+                continue
+        # Increment istation to move to the next station
+        istation += 1
+
+    return array, headerdata
 
 
 def get_field_at(ds, field: str, coords: tuple):

tests/test_realtime_pollen_calibration/test_utils.py

@@ -38,8 +38,10 @@ def test_interpolation():
         encode_cf=("time", "geography", "vertical"),
     )
     obs_mod_data = utils.read_atab(
-        "ALNU",
-        str(here()) + "/RTcal_testdata/alnu_pollen_measured_values_2022022207.atab",
+        pollen_type="ALNU",
+        max_miss_stns=2,
+        file_obs_stns=str(here())
+        + "/RTcal_testdata/alnu_pollen_measured_values_2022022207.atab",
     )
     ######################
     # Specify the test

tools/get_data_and_test.sh

@@ -19,6 +19,7 @@ const_file : $localpath/RTcal_testdata/CLON_CLAT_ICON-CH1.gb2
 station_obs_file : $localpath/RTcal_testdata/pollen_measured_2024020118.atab
 station_mod_file : $localpath/RTcal_testdata/pollen_modelled_2024020118.atab
 hour_incr : 1
+max_miss_stns : 4
 EOF
 
 # activate conda env