Add ismip7_forcing test group

compass/landice/__init__.py

@@ -11,6 +11,7 @@
 from compass.landice.tests.hydro_radial import HydroRadial
 from compass.landice.tests.ismip6_forcing import Ismip6Forcing
 from compass.landice.tests.ismip6_run import Ismip6Run
+from compass.landice.tests.ismip7_forcing import Ismip7Forcing
 from compass.landice.tests.isunnguata_sermia import IsunnguataSermia
 from compass.landice.tests.kangerlussuaq import Kangerlussuaq
 from compass.landice.tests.koge_bugt_s import KogeBugtS
@@ -46,6 +47,7 @@ def __init__(self):
         self.add_test_group(HydroRadial(mpas_core=self))
         self.add_test_group(Ismip6Forcing(mpas_core=self))
         self.add_test_group(Ismip6Run(mpas_core=self))
+        self.add_test_group(Ismip7Forcing(mpas_core=self))
         self.add_test_group(IsunnguataSermia(mpas_core=self))
         self.add_test_group(Kangerlussuaq(mpas_core=self))
         self.add_test_group(KogeBugtS(mpas_core=self))

compass/landice/tests/ismip7_forcing/__init__.py

@@ -0,0 +1,24 @@
+from compass.landice.tests.ismip7_forcing.atmosphere import Atmosphere
+from compass.landice.tests.ismip7_forcing.ocean_thermal import OceanThermal
+from compass.testgroup import TestGroup
+
+
+class Ismip7Forcing(TestGroup):
+    """
+    A test group for processing ISMIP7 forcing data
+    for the Antarctic Ice Sheet
+    """
+
+    def __init__(self, mpas_core):
+        """
+        Create the test group
+
+        Parameters
+        ----------
+        mpas_core : compass.landice.Landice
+            the MPAS core that this test group belongs to
+        """
+        super().__init__(mpas_core=mpas_core, name="ismip7_forcing")
+
+        self.add_test_case(Atmosphere(test_group=self))
+        self.add_test_case(OceanThermal(test_group=self))

compass/landice/tests/ismip7_forcing/atmosphere/__init__.py

@@ -0,0 +1,53 @@
+from compass.landice.tests.ismip7_forcing.atmosphere.process_runoff import (
+    ProcessRunoff,
+)
+from compass.landice.tests.ismip7_forcing.atmosphere.process_smb import (
+    ProcessSmb,
+)
+from compass.landice.tests.ismip7_forcing.atmosphere.process_smb_gradient import (  # noqa: E501
+    ProcessSmbGradient,
+)
+from compass.landice.tests.ismip7_forcing.atmosphere.process_temperature import (  # noqa: E501
+    ProcessTemperature,
+)
+from compass.landice.tests.ismip7_forcing.atmosphere.process_temperature_gradient import (  # noqa: E501
+    ProcessTemperatureGradient,
+)
+from compass.landice.tests.ismip7_forcing.configure import (
+    configure as configure_testgroup,
+)
+from compass.testcase import TestCase
+
+
+class Atmosphere(TestCase):
+    """
+    A test case for processing ISMIP7 atmosphere forcing data.
+    Remaps monthly SMB, temperature, and annual gradients (SMB and
+    temperature) from the ISMIP7 polar stereographic grid to the
+    MALI unstructured mesh. Also processes runoff (mrro).
+    """
+
+    def __init__(self, test_group):
+        """
+        Create the test case
+
+        Parameters
+        ----------
+        test_group : compass.landice.tests.ismip7_forcing.Ismip7Forcing
+            The test group that this test case belongs to
+        """
+        name = "atmosphere"
+        subdir = name
+        super().__init__(test_group=test_group, name=name, subdir=subdir)
+
+        self.add_step(ProcessSmb(test_case=self))
+        self.add_step(ProcessTemperature(test_case=self))
+        self.add_step(ProcessSmbGradient(test_case=self))
+        self.add_step(ProcessTemperatureGradient(test_case=self))
+        self.add_step(ProcessRunoff(test_case=self))
+
+    def configure(self):
+        """
+        Configures test case
+        """
+        configure_testgroup(config=self.config)

compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py

@@ -0,0 +1,272 @@
+import glob
+import os
+import shutil
+
+import numpy as np
+import xarray as xr
+from mpas_tools.io import write_netcdf
+from mpas_tools.logging import check_call
+from scipy.ndimage import distance_transform_edt
+
+from compass.landice.tests.ismip7_forcing.create_mapfile import (
+    build_mapping_file,
+)
+from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params
+from compass.step import Step
+
+
+class ProcessRunoff(Step):
+    """
+    A step for processing ISMIP7 ice sheet runoff (mrro) data.
+    Remaps monthly runoff from the ISMIP7 polar stereographic
+    grid to the MALI unstructured mesh. GrIS only.
+    """
+
+    def __init__(self, test_case):
+        """
+        Create the step
+
+        Parameters
+        ----------
+        test_case : compass.landice.tests.ismip7_forcing.atmosphere.Atmosphere
+            The test case this step belongs to
+        """
+        super().__init__(test_case=test_case, name="process_runoff")
+
+    def setup(self):
+        """
+        Set up this step of the test case
+        """
+        config = self.config
+        section = config["ismip7"]
+        base_path_mali = section.get("base_path_mali")
+        mali_mesh_file = section.get("mali_mesh_file")
+
+        self.add_input_file(filename=mali_mesh_file,
+                            target=os.path.join(base_path_mali,
+                                                mali_mesh_file))
+
+    def run(self):
+        """
+        Run this step of the test case
+        """
+        logger = self.logger
+        config = self.config
+        params = get_params(config)
+
+        section = config["ismip7"]
+        base_path_ismip7 = section.get("base_path_ismip7")
+        mali_mesh_name = section.get("mali_mesh_name")
+        mali_mesh_file = section.get("mali_mesh_file")
+        model = section.get("model")
+        scenario = section.get("scenario")
+        output_base_path = section.get("output_base_path")
+        ice_sheet = section.get("ice_sheet")
+
+        section = config["ismip7_atmosphere"]
+        method_remap = section.get("method_remap")
+        start_year = section.getint("start_year")
+        end_year = section.getint("end_year")
+
+        # Discover input files
+        prefix = params['prefix']
+        resolution = params['atm_resolution']
+        version = params['atm_version']
+        input_path = os.path.join(base_path_ismip7, "mrro", version)
+        file_pattern = (f"mrro_{prefix}_{model}_{scenario}_"
+                        f"SDBN1-{resolution}_{version}_*.nc")
+        all_files = sorted(glob.glob(os.path.join(input_path, file_pattern)))
+
+        if not all_files:
+            raise FileNotFoundError(
+                f"No runoff files found matching pattern:\n"
+                f"  {os.path.join(input_path, file_pattern)}")
+
+        # Filter to requested year range
+        input_files = []
+        for f in all_files:
+            year = int(os.path.basename(f).split("_")[-1].replace(".nc", ""))
+            if start_year <= year <= end_year:
+                input_files.append(f)
+
+        if not input_files:
+            raise FileNotFoundError(
+                f"No runoff files found for year range "
+                f"{start_year}-{end_year}")
+
+        logger.info(f"Found {len(input_files)} runoff files for years "
+                    f"{start_year}-{end_year}")
+
+        # Build mapping file (reuse if already created by other atm steps)
+        mapping_file = (f"map_ismip7_{ice_sheet}_atm_to_"
+                        f"{mali_mesh_name}_{method_remap}.nc")
+
+        if not os.path.exists(mapping_file):
+            logger.info("Building mapping file...")
+            build_mapping_file(config, logger,
+                               input_files[0], mapping_file,
+                               mali_mesh_file=mali_mesh_file,
+                               method_remap=method_remap)
+
+        # Remap each year file
+        remapped_files = []
+        for input_file in input_files:
+            basename = os.path.basename(input_file)
+            remapped_file = f"remapped_{basename}"
+            remapped_files.append(remapped_file)
+
+            if os.path.exists(remapped_file):
+                logger.info(f"  Remapped file exists, skipping: {basename}")
+                continue
+
+            # Extrapolate fill values on source grid before remapping
+            # so they don't pollute neighboring cells during interpolation
+            extrap_file = f"extrap_{basename}"
+            if not os.path.exists(extrap_file):
+                self._extrapolate_source(input_file, extrap_file, "mrro",
+                                         logger)
+
+            logger.info(f"  Remapping: {basename}")
+            args = ["ncremap",
+                    "-i", extrap_file,
+                    "-o", remapped_file,
+                    "-m", mapping_file,
+                    "-v", "mrro"]
+
+            check_call(args, logger=logger)
+
+            # Clean up extrapolated source file
+            os.remove(extrap_file)
+
+        # Combine remapped files and rename to MALI conventions
+        logger.info("Combining remapped files and renaming variables...")
+        output_file = (f"{mali_mesh_name}_runoff_{model}_{scenario}_"
+                       f"{start_year}-{end_year}.nc")
+
+        self._combine_and_rename(remapped_files, output_file)
+
+        # Clean up remapped files
+        logger.info("Cleaning up temporary remapped files...")
+        for f in remapped_files:
+            if os.path.exists(f):
+                os.remove(f)
+
+        # Place output in appropriate directory
+        output_path = os.path.join(output_base_path, "atmosphere_forcing",
+                                   f"{model}_{scenario}")
+        if not os.path.exists(output_path):
+            os.makedirs(output_path)
+
+        dst = os.path.join(output_path, output_file)
+        shutil.copy(output_file, dst)
+
+        logger.info(f"Done. Output: {dst}")
+
+    def _combine_and_rename(self, remapped_files, output_file):
+        """
+        Combine yearly remapped files and rename variables/dimensions
+        to MALI conventions.
+
+        Parameters
+        ----------
+        remapped_files : list of str
+            List of remapped NetCDF file paths
+
+        output_file : str
+            Output file path
+        """
+        ds = xr.open_mfdataset(remapped_files, concat_dim="time",
+                               combine="nested", engine="netcdf4")
+
+        # Rename dimensions to MALI conventions
+        rename_dims = {}
+        if "time" in ds.dims:
+            rename_dims["time"] = "Time"
+        if "ncol" in ds.dims:
+            rename_dims["ncol"] = "nCells"
+        if rename_dims:
+            ds = ds.rename(rename_dims)
+
+        # Rename variable
+        if "mrro" in ds:
+            ds = ds.rename({"mrro": "ismip6Runoff"})
+
+        # Add xtime variable with monthly timestamps
+        xtime = []
+        for t_index in range(ds.sizes["Time"]):
+            date = ds.Time[t_index]
+            yr = int(date.dt.year.values)
+            mo = int(date.dt.month.values)
+            date_str = f"{yr:04d}-{mo:02d}-15_00:00:00".ljust(64)
+            xtime.append(date_str)
+
+        ds["xtime"] = ("Time", xtime)
+        ds["xtime"] = ds.xtime.astype("S")
+
+        # Set attributes
+        ds["ismip6Runoff"].attrs = {
+            "long_name": "ice sheet runoff",
+            "units": "kg m-2 s-1",
+        }
+
+        # Drop auxiliary variables from remapping
+        vars_to_drop = [v for v in ["lon", "lon_vertices", "lat",
+                                    "lat_vertices", "area"]
+                        if v in ds]
+        if vars_to_drop:
+            ds = ds.drop_vars(vars_to_drop)
+
+        write_netcdf(ds, output_file)
+
+    def _extrapolate_source(self, input_file, output_file, varname, logger):
+        """
+        Extrapolate fill/missing values on the source polar stereographic
+        grid using nearest-neighbor via distance_transform_edt. This must
+        be done before remapping so that fill values don't contaminate the
+        interpolation stencil.
+
+        Parameters
+        ----------
+        input_file : str
+            Path to the input NetCDF file on the source grid
+
+        output_file : str
+            Path to write the extrapolated file
+
+        varname : str
+            Name of the variable to extrapolate (e.g., "mrro")
+
+        logger : logging.Logger
+            Logger for status messages
+        """
+        logger.info(f"    Extrapolating fill values on source grid: "
+                    f"{os.path.basename(input_file)}")
+
+        ds = xr.open_dataset(input_file, engine="netcdf4")
+        data = ds[varname]
+
+        # Process each time step
+        # Source files have dims like (time, y, x)
+        values = data.values.copy()
+        non_spatial_shape = values.shape[:-2]  # (time,)
+
+        for idx in np.ndindex(non_spatial_shape):
+            slab = values[idx]  # shape (ny, nx)
+            valid_mask = np.isfinite(slab)
+            if valid_mask.all() or not valid_mask.any():
+                continue
+            nearest_inds = distance_transform_edt(
+                ~valid_mask, return_distances=False, return_indices=True)
+            invalid = ~valid_mask
+            values[idx][invalid] = slab[
+                nearest_inds[0, invalid],
+                nearest_inds[1, invalid]]
+
+        ds[varname] = (data.dims, values)
+        ds[varname].attrs = data.attrs
+
+        # Remove _FillValue encoding so output has no masked values
+        if "_FillValue" in ds[varname].encoding:
+            del ds[varname].encoding["_FillValue"]
+
+        write_netcdf(ds, output_file)