Preprocess and assimilate the GHCN snow data

.gitignore

@@ -175,6 +175,7 @@ ush/global_cycle_driver.sh
 ush/gsi_satbias2ioda_all.sh
 ush/jediinc2fv3.py
 ush/imsfv3_scf2ioda.py
+ush/ghcn_snod2ioda.py
 ush/bufr_snocvr_snomad.py
 ush/atparse.bash
 ush/bufr2ioda_insitu*

dev/parm/config/gfs/config.snowanl.j2

@@ -14,9 +14,13 @@ export APPLY_INCR_NML_TMPL="${PARMgfs}/gdas/snow/apply_incr_nml.j2"
 
 export TASK_CONFIG_YAML="${PARMgfs}/gdas/snow/snow_det_config.yaml.j2"
 export OBS_LIST_YAML="${PARMgfs}/gdas/snow/snow_obs_list.yaml.j2"
+export ims_scf_obs_suffix="asc" # asc-ascii; nc-netcdf
+export fail_on_missing=False    # False: just warn; True: fail & exit
 
 export PREP_SNOCVR_SNOMAD_YAML="${PARMgfs}/gdas/snow/prep/prep_snocvr_snomad.yaml.j2"
 export OBSBUILDER="${USHgfs}/bufr_snocvr_snomad.py"
+export PREP_GHCN_YAML="${PARMgfs}/gdas/snow/prep/prep_ghcn.yaml.j2"
+export GHCN2IODACONV="${USHgfs}/ghcn_snod2ioda.py"
 
 export io_layout_x="{{ IO_LAYOUT_X }}"
 export io_layout_y="{{ IO_LAYOUT_Y }}"

dev/scripts/exglobal_snow_analysis.py

@@ -31,6 +31,10 @@
     if snow_anl.task_config.DO_IMS_SCF:
         snow_anl.execute('scf_to_ioda')
 
+    # Process GHCN (if applicable)
+    if snow_anl.task_config.DO_GHCN:
+        snow_anl.prepare_GHCN()
+
     # Execute JEDI snow analysis
     snow_anl.execute('snowanlvar')
 

sorc/link_workflow.sh

@@ -276,6 +276,7 @@ if [[ -d "${HOMEgfs}/sorc/gdas.cd/build" ]]; then
     ${LINK_OR_COPY} "${HOMEgfs}/sorc/gdas.cd/ush/gsi_satbias2ioda_all.sh" .
     ${LINK_OR_COPY} "${HOMEgfs}/sorc/gdas.cd/ush/snow/bufr_snocvr_snomad.py" .
     ${LINK_OR_COPY} "${HOMEgfs}/sorc/gdas.cd/build/bin/imsfv3_scf2ioda.py" .
+    ${LINK_OR_COPY} "${HOMEgfs}/sorc/gdas.cd/ush/snow/ghcn_snod2ioda.py" .
 fi
 
 #------------------------------
@@ -439,7 +440,6 @@ fi
 # GDASApp executables
 if [[ -d "${HOMEgfs}/sorc/gdas.cd/install" ]]; then
     cp -f "${HOMEgfs}/sorc/gdas.cd/install/bin"/gdas* ./
-    cp -f "${HOMEgfs}/sorc/gdas.cd/install/bin/calcfIMS.exe" ./gdas_calcfIMS.x
     cp -f "${HOMEgfs}/sorc/gdas.cd/install/bin/apply_incr.exe" ./gdas_apply_incr.x
 fi
 

ush/python/pygfs/task/snow_analysis.py

@@ -48,14 +48,46 @@ def __init__(self, config: Dict[str, Any]):
         super().__init__(config)
 
         _res = int(self.task_config['CASE'][1:])
+        _fail_on_missing = str(self.task_config.fail_on_missing[0]).lower() == "true" \
+            if isinstance(self.task_config.fail_on_missing, list) \
+            else bool(self.task_config.fail_on_missing)
 
         # if 00z, do SCF preprocessing
-        _ims_file = os.path.join(self.task_config.COMIN_OBS, f'{self.task_config.OPREFIX}imssnow96.asc')
+        _ims_file = os.path.join(
+            self.task_config.COMIN_OBS,
+            f'{self.task_config.OPREFIX}imssnow96.{self.task_config.ims_scf_obs_suffix}'
+        )
         logger.info(f"Checking for IMS file: {_ims_file}")
-        if self.task_config.cyc == 0 and os.path.exists(_ims_file):
-            _DO_IMS_SCF = True
+        _DO_IMS_SCF = False
+        if self.task_config.cyc == 0:
+            if os.path.exists(_ims_file):
+                _DO_IMS_SCF = True
+            else:
+                if _fail_on_missing:
+                    raise FileNotFoundError(
+                        f"IMS obs file required but not found: {_ims_file}"
+                    )
+                else:
+                    logger.warning(f"IMS obs file missing: {_ims_file}")
+        else:
+            logger.info("Not 00z cycle — Skipping IMS preprocessing.")
+
+        # if 00z, do GHCN preprocessing
+        _ghcn_file = os.path.join(self.task_config.COMIN_OBS, f'{self.task_config.OPREFIX}ghcn_snow.csv')
+        logger.info(f"Checking for GHCN csv file: {_ghcn_file}")
+        _DO_GHCN = False
+        if self.task_config.cyc == 0:
+            if os.path.exists(_ghcn_file):
+                _DO_GHCN = True
+            else:
+                if _fail_on_missing:
+                    raise FileNotFoundError(
+                        f"GHCN obs file required but not found: {_ghcn_file}"
+                    )
+                else:
+                    logger.warning(f"GHCN obs file missing: {_ghcn_file}")
         else:
-            _DO_IMS_SCF = False
+            logger.info("Not 00z cycle — Skipping GHCN preprocessing.")
 
         # Extend task_config with variables repeatedly used across this class
         self.task_config.update(AttrDict(
@@ -68,6 +100,7 @@ def __init__(self, config: Dict[str, Any]):
                 'snow_prepobs_path': os.path.join(self.task_config.DATA, 'prep'),
                 'ims_file': _ims_file,
                 'DO_IMS_SCF': _DO_IMS_SCF,  # Boolean to decide if IMS snow cover processing is done
+                'DO_GHCN': _DO_GHCN,  # Boolean to decide if GHCN processing is done
             }
         ))
 
@@ -222,6 +255,74 @@ def prepare_SNOCVR_SNOMAD(self) -> None:
             logger.info(f"Copy {output_file} successfully generated")
             FileHandler(prep_snocvr_snomad_config.netcdf).sync()
 
+    @logit(logger)
+    def prepare_GHCN(self) -> None:
+        """Prepare the GHCN data for a global snow analysis
+
+        This method will prepare GHCN data for a global snow analysis using JEDI.
+        This includes:
+        - creating GHCN snowdepth data in IODA format.
+
+        Parameters
+        ----------
+        Analysis: parent class for GDAS task
+
+        Returns
+        ----------
+        None
+        """
+
+        # Read and render the prep_ghcn.yaml.j2
+        logger.info(f"Reading {self.task_config.PREP_GHCN_YAML}")
+        prep_ghcn_config = parse_j2yaml(self.task_config.PREP_GHCN_YAML, self.task_config)
+        logger.debug(f"{self.task_config.PREP_GHCN_YAML}:\n{pformat(prep_ghcn_config)}")
+
+        # Define these locations in gdas/snow/prep/prep_ghcn.yaml.j2
+        logger.info("Copying GHCN obs to DATA")
+        FileHandler(prep_ghcn_config.stage).sync()
+
+        # Execute ioda converter to create the GHCN obs data in IODA format
+        logger.info("Create GHCN obs data in IODA format")
+
+        csv_file = f'{self.task_config.OPREFIX}ghcn_snow.csv'
+        station_file = f'ghcnd-stations.txt'
+        output_file = f'{self.task_config.OPREFIX}ghcn_snow.nc'
+        if os.path.exists(f"{os.path.join(self.task_config.DATA, output_file)}"):
+            rm_p(output_file)
+        if not os.path.isfile(csv_file):
+            logger.warning(f"WARNING: GHCN obs file not found.")
+            return
+
+        logger.info("Link GHCN2IODACONV into DATA/")
+        exe_src = self.task_config.GHCN2IODACONV
+        exe_dest = os.path.join(self.task_config.DATA, os.path.basename(exe_src))
+        if os.path.exists(exe_dest):
+            rm_p(exe_dest)
+        os.symlink(exe_src, exe_dest)
+
+        exe = Executable(exe_dest)
+        exe.add_default_arg(["-i", f"{os.path.join(self.task_config.DATA, csv_file)}"])
+        exe.add_default_arg(["-o", f"{os.path.join(self.task_config.DATA, output_file)}"])
+        exe.add_default_arg(["-f", f"{os.path.join(self.task_config.DATA, station_file)}"])
+        exe.add_default_arg(["-d", f"{to_YMDH(self.task_config.current_cycle)}"])
+        try:
+            logger.debug(f"Executing {exe}")
+            exe()
+        except OSError:
+            logger.exception(f"Failed to execute {exe}")
+            raise
+        except Exception as err:
+            logger.exception(f"An error occured during execution of {exe}")
+            raise WorkflowException(f"An error occured during execution of {exe}") from err
+
+        # Ensure the IODA snow depth GHCN file is produced by the IODA converter
+        # If so, copy to DATA/prep/
+        if not os.path.isfile(f"{os.path.join(self.task_config.DATA, output_file)}"):
+            logger.warning(f"{output_file} not produced - continuing without it.")
+        else:
+            logger.info(f"Copy {output_file} successfully generated")
+            FileHandler(prep_ghcn_config.ghcn2ioda).sync()
+
     @logit(logger)
     def add_increments(self) -> None:
         """Executes the program "apply_incr.exe" to create analysis "sfc_data" files by adding increments to backgrounds

ush/python/pygfs/task/snowens_analysis.py

@@ -51,14 +51,46 @@ def __init__(self, config: Dict[str, Any]):
         super().__init__(config)
 
         _res = int(self.task_config['CASE_ENS'][1:])
+        _fail_on_missing = str(self.task_config.fail_on_missing[0]).lower() == "true" \
+            if isinstance(self.task_config.fail_on_missing, list) \
+            else bool(self.task_config.fail_on_missing)
 
         # if 00z, do SCF preprocessing
-        _ims_file = os.path.join(self.task_config.COMIN_OBS, f'{self.task_config.OPREFIX}imssnow96.asc')
+        _ims_file = os.path.join(
+            self.task_config.COMIN_OBS,
+            f'{self.task_config.OPREFIX}imssnow96.{self.task_config.ims_scf_obs_suffix}'
+        )
         logger.info(f"Checking for IMS file: {_ims_file}")
-        if self.task_config.cyc == 0 and os.path.exists(_ims_file):
-            _DO_IMS_SCF = True
+        _DO_IMS_SCF = False
+        if self.task_config.cyc == 0:
+            if os.path.exists(_ims_file):
+                _DO_IMS_SCF = True
+            else:
+                if _fail_on_missing:
+                    raise FileNotFoundError(
+                        f"IMS obs file required but not found: {_ims_file}"
+                    )
+                else:
+                    logger.warning(f"IMS obs file missing: {_ims_file}")
+        else:
+            logger.info("Not 00z cycle — Skipping IMS preprocessing.")
+
+        # if 00z, do GHCN preprocessing
+        _ghcn_file = os.path.join(self.task_config.COMIN_OBS, f'{self.task_config.OPREFIX}ghcn_snow.csv')
+        logger.info(f"Checking for GHCN csv file: {_ghcn_file}")
+        _DO_GHCN = False
+        if self.task_config.cyc == 0:
+            if os.path.exists(_ghcn_file):
+                _DO_GHCN = True
+            else:
+                if _fail_on_missing:
+                    raise FileNotFoundError(
+                        f"GHCN obs file required but not found: {_ghcn_file}"
+                    )
+                else:
+                    logger.warning(f"GHCN obs file missing: {_ghcn_file}")
         else:
-            _DO_IMS_SCF = False
+            logger.info("Not 00z cycle — Skipping GHCN preprocessing.")
 
         # Extend task_config with variables repeatedly used across this class
         self.task_config.update(AttrDict(
@@ -69,8 +101,10 @@ def __init__(self, config: Dict[str, Any]):
                 'npz': self.task_config.LEVS - 1,
                 'CASE': self.task_config.CASE_ENS,
                 'snow_bkg_path': os.path.join('.', 'bkg', 'ensmean/'),
+                'snow_prepobs_path': os.path.join(self.task_config.DATA, 'prep'),
                 'ims_file': _ims_file,
                 'DO_IMS_SCF': _DO_IMS_SCF,  # Boolean to decide if IMS snow cover processing is done
+                'DO_GHCN': _DO_GHCN,  # Boolean to decide if GHCN processing is done
             }
         ))
 
@@ -226,6 +260,71 @@ def prepare_SNOCVR_SNOMAD(self) -> None:
             logger.info(f"Copy {output_file} successfully generated")
             FileHandler(prep_snocvr_snomad_config.netcdf).sync()
 
+    @logit(logger)
+    def prepare_GHCN(self) -> None:
+        """Prepare the GHCN data for a global snow analysis
+        This method will prepare GHCN data for a global snow analysis using JEDI.
+        This includes:
+        - creating GHCN snowdepth data in IODA format.
+        Parameters
+        ----------
+        Analysis: parent class for GDAS task
+        Returns
+        ----------
+        None
+        """
+
+        # Read and render the prep_ghcn.yaml.j2
+        logger.info(f"Reading {self.task_config.PREP_GHCN_YAML}")
+        prep_ghcn_config = parse_j2yaml(self.task_config.PREP_GHCN_YAML, self.task_config)
+        logger.debug(f"{self.task_config.PREP_GHCN_YAML}:\n{pformat(prep_ghcn_config)}")
+
+        # Define these locations in gdas/snow/prep/prep_ghcn.yaml.j2
+        logger.info("Copying GHCN obs to DATA")
+        FileHandler(prep_ghcn_config.stage).sync()
+
+        # Execute ioda converter to create the GHCN obs data in IODA format
+        logger.info("Create GHCN obs data in IODA format")
+
+        csv_file = f'{self.task_config.OPREFIX}ghcn_snow.csv'
+        station_file = f'ghcnd-stations.txt'
+        output_file = f'{self.task_config.OPREFIX}ghcn_snow.nc'
+        if os.path.exists(f"{os.path.join(self.task_config.DATA, output_file)}"):
+            rm_p(output_file)
+        if not os.path.isfile(csv_file):
+            logger.warning(f"WARNING: GHCN obs file not found.")
+            return
+
+        logger.info("Link GHCN2IODACONV into DATA/")
+        exe_src = self.task_config.GHCN2IODACONV
+        exe_dest = os.path.join(self.task_config.DATA, os.path.basename(exe_src))
+        if os.path.exists(exe_dest):
+            rm_p(exe_dest)
+        os.symlink(exe_src, exe_dest)
+
+        exe = Executable(exe_dest)
+        exe.add_default_arg(["-i", f"{os.path.join(self.task_config.DATA, csv_file)}"])
+        exe.add_default_arg(["-o", f"{os.path.join(self.task_config.DATA, output_file)}"])
+        exe.add_default_arg(["-f", f"{os.path.join(self.task_config.DATA, station_file)}"])
+        exe.add_default_arg(["-d", f"{to_YMDH(self.task_config.current_cycle)}"])
+        try:
+            logger.debug(f"Executing {exe}")
+            exe()
+        except OSError:
+            logger.exception(f"Failed to execute {exe}")
+            raise
+        except Exception as err:
+            logger.exception(f"An error occured during execution of {exe}")
+            raise WorkflowException(f"An error occured during execution of {exe}") from err
+
+        # Ensure the IODA snow depth GHCN file is produced by the IODA converter
+        # If so, copy to DATA/prep/
+        if not os.path.isfile(f"{os.path.join(self.task_config.DATA, output_file)}"):
+            logger.warning(f"{output_file} not produced - continuing without it.")
+        else:
+            logger.info(f"Copy {output_file} successfully generated")
+            FileHandler(prep_ghcn_config.ghcn2ioda).sync()
+
     @logit(logger)
     def add_increments(self) -> None:
         """Executes the program "apply_incr.exe" to create analysis "sfc_data" files by adding increments to backgrounds