Icon support

data/grib2_files_ICON-CH1/iaf2023042500

@@ -0,0 +1 @@
+/scratch/mch/paa/wd/23042500_icon_BETU_POAC_testchain/lm_coarse/iaf2023042500

data/grib2_files_ICON-CH1/lfff00000000c

@@ -0,0 +1 @@
+/oprusers/osm/opr.aare/data/ICON_INPUT/ICON-CH1-EPS/lfff00000000c

src/realtime_pollen_calibration/main.py

@@ -16,19 +16,20 @@
 
 cwd = Path(os.getcwd())
 data_path = str(cwd)
-# data_path = str(cwd.parent.parent.absolute())
 
-file_data = data_path + "/data/atabs/alnu_pollen_measured_values_2022021314.atab"
-file_grib = data_path + "/data/grib2_files_cosmo1e/laf2022021314_ALNUtthrs_tthre"
-file_out = data_path + "/data/2022021315_ALNUtthrs_tthre"
+file_obs   = data_path + "/data/atabs/alnu_pollen_measured_values_2022021314.atab"
+file_POV   = data_path + "/data/grib2_files_ICON-CH1/ART_POV_iconR19B08-grid_0001_all_specs_values"
+file_T_2M  = data_path + "/data/grib2_files_ICON-CH1/iaf2023042500"
+file_Const = data_path + "/data/grib2_files_ICON-CH1/lfff00000000c"
+file_out   = data_path + "/data/grib2_files_ICON-CH1/POV_out"
 update_phenology_realtime.update_phenology_realtime(
-    file_data, file_grib, file_out, False
+    file_obs, file_POV, file_T_2M, file_Const, file_out, True
 )
 
-file_data = data_path + "/data/atabs/alnu_pollen_measured_values_2022022207.atab"
-file_data_mod = data_path + "/data/atabs/alnu_pollen_modelled_values_2022022207.atab"
-file_grib = data_path + "/data/grib2_files_cosmo1e/laf2022022207_ALNUtune"
+file_obs = data_path + "/data/atabs/alnu_pollen_measured_values_2022022207.atab"
+file_obs_mod = data_path + "/data/atabs/alnu_pollen_modelled_values_2022022207.atab"
+file_POV = data_path + "/data/grib2_files_cosmo1e/laf2022022207_ALNUtune"
 file_out = data_path + "/data/2022022208_ALNUtune"
 update_strength_realtime.update_strength_realtime(
-    file_data, file_data_mod, file_grib, file_out, False
+    file_obs, file_obs_mod, file_POV, file_out, False
 )

src/realtime_pollen_calibration/update_phenology_realtime.py

@@ -7,28 +7,139 @@
 """A module for the update of the pollen start and end of season in real time."""
 
 # Third-party
-import cfgrib  # type: ignore
 import numpy as np
+import xarray as xr
+from datetime import datetime
+from eccodes import (
+    codes_get,
+    codes_get_array,
+    codes_grib_new_from_file,
+    codes_release,
+    codes_set_array,
+    codes_write,
+)
 
 # First-party
 from realtime_pollen_calibration import utils
 
 
-def update_phenology_realtime(file_obs, file_in, file_out, verbose=False):
+def update_phenology_realtime(file_obs, file_POV, file_T_2M, file_Const, file_out, verbose=False):
     """Advance the temperature threshold fields by one hour.
 
     Args:
         file_obs: Location of ATAB file containing the pollen
                 concentration information at the stations.
-        file_in: Location of GRIB2 file containing the following fields:
-                'T_2M', 'tthrs', 'tthre' (for POAC, 'saisl' instead),
+        file_POV: Location of ICON GRIB2 file containing the pollen fields:
+                'tthrs', 'tthre' (for POAC, 'saisl' instead),
                 'saisn' and 'ctsum'. Lat-lon information of the grid must be
                 present in the file.
+        file_T_2M: Location of GRIB2 file containing T_2M.
+        file_Const: Location of GRIB2 file containing Longitudes and Latitudes of the 
+                unstructured ICON grid.
         file_out: Location of the desired output file.
         verbose: Optional additional debug prints.
 
     """
-    ds = cfgrib.open_dataset(file_in, encode_cf=("time", "geography", "vertical"))
+
+    #file_POV = "data/grib2_files_ICON-CH1/ART_POV_iconR19B08-grid_0001_all_specs_values"
+    #file_obs= "data/atabs/alnu_pollen_measured_values_2022021314.atab"
+
+    fh_POV = open(file_POV, "rb")
+    fh_Const = open(file_Const, "rb")
+    fh_T_2M = open(file_T_2M, "rb")
+
+    # read CLON, CLAT
+    while True:
+        # Get the next message
+        recConst = codes_grib_new_from_file(fh_Const)
+        if recConst is None:
+            break
+
+        # Get the short name of the current field
+        short_name = codes_get(recConst, "shortName")
+
+        # Extract longitude and latitude of the ICON grid
+        if short_name == "CLON":
+            CLON = codes_get_array(recConst, "values")
+        if short_name == "CLAT":
+            CLAT = codes_get_array(recConst, "values")
+
+        # Delete the message
+        codes_release(recConst)
+
+    # Close the GRIB file
+    fh_Const.close()
+
+    # read POV and extract available fields
+    specs = ["ALNU", "BETU", "POAC", "CORY"]
+    fields= ["tthrs", "tthre", "saisn", "ctsum"]
+    pol_fields = [x + y for x in specs for y in fields]
+    pol_fields[9] = "POACsaisl"
+
+    calFields = {}
+    while True:
+        # Get the next message
+        recPOV = codes_grib_new_from_file(fh_POV)
+        if recPOV is None:
+            break
+
+        # Get the short name of the current field
+        short_name = codes_get(recPOV, "shortName")
+
+        # Extract and alter fields if they present
+        for pol_field in pol_fields:
+            if short_name == pol_field:
+                calFields[pol_field] = codes_get_array(recPOV, "values")
+
+        # Delete the message
+        codes_release(recPOV)
+
+    # Close the GRIB files
+    fh_POV.close()
+
+
+    # Read T_2M
+    while True:
+        # Get the next message
+        recX = codes_grib_new_from_file(fh_T_2M)
+        if recX is None:
+            break
+
+        # Get the short name of the current field
+        short_name = codes_get(recX, "shortName")
+
+        if short_name == "T_2M":
+            calFields["T_2M"] = codes_get_array(recX, "values")
+
+            # timestamp is needed. Take it from the T_2M field
+            dataDate = str(codes_get(recX, "dataDate"))
+            dataTime = str(str(codes_get(recX, "dataTime")).zfill(2))
+            dataDateTime = dataDate + dataTime
+
+            # Convert the string to a datetime object
+            date_obj = datetime.strptime(dataDateTime, '%Y%m%d%H')
+            date_obj_fmt = date_obj.strftime('%Y-%m-%dT%H:00:00.000000000')
+            time_values = np.datetime64(date_obj_fmt)
+
+        codes_release(recX)
+
+    # Close the GRIB files
+    fh_T_2M.close()
+
+    # Dictionary to hold DataArrays for each variable
+    calFields_arrays = {}
+
+    # Loop through variables to create DataArrays
+    for var_name, data in calFields.items():
+        data_array = xr.DataArray(data, coords={'index': np.arange(len(data))}, dims=['index'])
+        data_array.coords['latitude'] = (('index'), CLAT)
+        data_array.coords['longitude'] = (('index'), CLON)
+        data_array.coords['time'] = time_values
+        calFields_arrays[var_name] = data_array
+
+    # Create an xarray Dataset with the DataArrays
+    ds = xr.Dataset(calFields_arrays)
+
     ptype_present = utils.get_pollen_type(ds)
     if verbose:
         print(f"Detected pollen types in the DataSet provided: {ptype_present}")
@@ -55,4 +166,4 @@ def update_phenology_realtime(file_obs, file_in, file_out, verbose=False):
                     obs_mod_data.coord_stns,
                     method="sum",
                 )
-    utils.to_grib(file_in, file_out, dict_fields)
+    utils.to_grib(file_POV, file_out, dict_fields)

src/realtime_pollen_calibration/utils.py

@@ -30,7 +30,7 @@
 thr_con_24 = {"ALNU": 240, "BETU": 240, "POAC": 72, "CORY": 240}
 thr_con_120 = {"ALNU": 720, "BETU": 720, "POAC": 216, "CORY": 720}
 failsafe = {"ALNU": 1000, "BETU": 2500, "POAC": 6000, "CORY": 2500}
-jul_days_excl = {"ALNU": 14, "BETU": 40, "POAC": 3, "CORY": 46}
+jul_days_excl = {"ALNU": 14, "BETU": 40, "POAC": 46, "CORY": 3}
 
 
 def count_to_log_level(count: int) -> int:
@@ -414,15 +414,15 @@ def get_change_phenol(  # pylint: disable=R0912,R0914,R0915
                 f"and last 120H {sum_obs}",
             )
             print(
-                f"Cumulative temperature sum {ctsum_stns.values[0][0]} ",
-                f"and threshold (start): {tthrs_stns.values[0][0]}",
-                f" and saisn: {saisn_stns.values[0][0]}",
+                f"Cumulative temperature sum {ctsum_stns.values[0]} ",
+                f"and threshold (start): {tthrs_stns.values[0]}",
+                f" and saisn: {saisn_stns.values[0]}",
             )
             if pollen_type != "POAC":
-                print(f"Cumsum temp threshold end: {tthre_stns.values[0][0]}")
+                print(f"Cumsum temp threshold end: {tthre_stns.values[0]}")
             else:
-                print(f"Saisl: {saisl_stns.values[0][0]}")
-            print(f"Temperature at station {t_2m_stns.values[0][0]}, " f"date: {date}")
+                print(f"Saisl: {saisl_stns.values[0]}")
+            print(f"Temperature at station {t_2m_stns.values[0]}, " f"date: {date}")
             print("-----------------------------------------")
         # ADJUSTMENT OF SEASON START AND END AT THE BEGINNING OF THE SEASON
         if (
@@ -545,15 +545,20 @@ def to_grib(inp: str, outp: str, dict_fields: dict) -> None:
                 break
             # clone record
             clone_id = eccodes.codes_clone(gid)
-            # get short_name
 
+            # get short_name
             short_name = eccodes.codes_get_string(clone_id, "shortName")
+
             # read values
             values = eccodes.codes_get_values(clone_id)
             eccodes.codes_set(
                 clone_id, "dataTime", eccodes.codes_get(clone_id, "dataTime") + 100
             )
             if short_name in dict_fields:
+
+                #set values in dict_fields[short_name] to zero where values is zero (edge values)
+                #This is because COSMO-1E was slightly smaller than ICON-CH1
+                dict_fields[short_name][values == 0] = 0
                 eccodes.codes_set_values(clone_id, dict_fields[short_name].flatten())
             else:
                 eccodes.codes_set_values(clone_id, values)