From 9d93fd17d9f17a7ce7128a4973ee64b6ca18ee0d Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 3 Sep 2026 09:11:51 -0700 Subject: [PATCH 1/6] Add ISMIP7 GrIS OCX (reanalysis) forcing support Handle the OCX scenario in the ismip7_forcing test group. OCX has no distinct ESM model: it uses RACMO2.3p2-ERA (atmosphere) and EN4 (ocean) at data version v1, with a named grid resolution (ocean-1000m) in the ocean file names. get_params() now applies OCX-specific overrides (sources, versions, ocean grid token) when scenario == OCX, so the [ismip7] model option is ignored and a single config file drives both the atmosphere and ocean_thermal test cases. Adds ismip7_forcing_ocx_gis.cfg and updates the user and developer docs. --- .../atmosphere/process_runoff.py | 2 + .../ismip7_forcing/atmosphere/process_smb.py | 2 + .../atmosphere/process_smb_gradient.py | 2 + .../atmosphere/process_temperature.py | 2 + .../process_temperature_gradient.py | 2 + .../tests/ismip7_forcing/ice_sheet_params.py | 32 +++++++- .../ismip7_forcing/ismip7_forcing_ocx_gis.cfg | 77 +++++++++++++++++++ .../ocean_thermal/process_thermal_forcing.py | 5 +- .../landice/test_groups/ismip7_forcing.rst | 14 +++- .../landice/test_groups/ismip7_forcing.rst | 22 ++++++ 10 files changed, 156 insertions(+), 4 deletions(-) create mode 100644 compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py index e03e85ecd6..ef7ac82164 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -72,6 +72,8 @@ def run(self): prefix = params['prefix'] resolution = params['atm_resolution'] version = params['atm_version'] + if params['atm_model'] is not None: + model = params['atm_model'] input_path = os.path.join(base_path_ismip7, "mrro", version) file_pattern = (f"mrro_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py index 222fb6e92e..014e098bc9 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -69,6 +69,8 @@ def run(self): prefix = params['prefix'] resolution = params['atm_resolution'] version = params['atm_version'] + if params['atm_model'] is not None: + model = params['atm_model'] input_path = os.path.join(base_path_ismip7, "acabf", version) file_pattern = (f"acabf_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py index 2f3e4c2cb4..6d39bb304e 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -70,6 +70,8 @@ def run(self): prefix = params['prefix'] resolution = params['atm_resolution'] version = params['atm_version'] + if params['atm_model'] is not None: + model = params['atm_model'] input_path = os.path.join(base_path_ismip7, "dacabfdz", version) file_pattern = (f"dacabfdz_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index caa8f6be9f..13ec7c15d9 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -69,6 +69,8 @@ def run(self): prefix = params['prefix'] resolution = params['atm_resolution'] version = params['atm_version'] + if params['atm_model'] is not None: + model = params['atm_model'] input_path = os.path.join(base_path_ismip7, "ts", version) file_pattern = (f"ts_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py index 5dec7b78bc..7399eda925 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -71,6 +71,8 @@ def run(self): prefix = params['prefix'] resolution = params['atm_resolution'] version = params['atm_version'] + if params['atm_model'] is not None: + model = params['atm_model'] input_path = os.path.join(base_path_ismip7, "dtsdz", version) file_pattern = (f"dtsdz_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") diff --git a/compass/landice/tests/ismip7_forcing/ice_sheet_params.py b/compass/landice/tests/ismip7_forcing/ice_sheet_params.py index 8324236c56..0d57094e13 100644 --- a/compass/landice/tests/ismip7_forcing/ice_sheet_params.py +++ b/compass/landice/tests/ismip7_forcing/ice_sheet_params.py @@ -10,8 +10,11 @@ 'atm_resolution': '2000m', 'atm_version': 'v2', 'ocean_version': 'v3', + 'ocean_grid': 'ocean', 'ocean_3d': True, 'ocean_temporal': 'decade', + 'atm_model': None, + 'ocean_model': None, }, 'gis': { 'projection': 'gis-bamber', @@ -19,8 +22,25 @@ 'atm_resolution': '1000m', 'atm_version': 'v2', 'ocean_version': 'v2', + 'ocean_grid': 'ocean', 'ocean_3d': False, 'ocean_temporal': 'yearly', + 'atm_model': None, + 'ocean_model': None, + }, +} + +# Overrides applied for the OCX (reanalysis) scenario. OCX has no distinct +# ESM model: it uses fixed reanalysis products (RACMO for the atmosphere and +# EN4 for the ocean) at data version v1. When scenario is 'OCX' the [ismip7] +# model option is ignored and these sources are used instead. +_OCX_OVERRIDES = { + 'gis': { + 'atm_version': 'v1', + 'ocean_version': 'v1', + 'ocean_grid': 'ocean-1000m', + 'atm_model': 'RACMO2.3p2-ERA', + 'ocean_model': 'EN4', }, } @@ -44,4 +64,14 @@ def get_params(config): raise ValueError( f"Unknown ice_sheet '{ice_sheet}'. " f"Must be one of: {list(_PARAMS.keys())}") - return _PARAMS[ice_sheet] + params = dict(_PARAMS[ice_sheet]) + + scenario = config.get("ismip7", "scenario") + if scenario == "OCX": + if ice_sheet not in _OCX_OVERRIDES: + raise ValueError( + f"The OCX scenario is not yet supported for ice_sheet " + f"'{ice_sheet}'.") + params.update(_OCX_OVERRIDES[ice_sheet]) + + return params diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg new file mode 100644 index 0000000000..d3658027db --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg @@ -0,0 +1,77 @@ +# Example config for processing the GrIS ISMIP7 OCX (reanalysis) forcing. +# OCX has no distinct ESM model: it uses RACMO2.3p2-ERA for the atmosphere and +# EN4 for the ocean, selected automatically when scenario = OCX. The [ismip7] +# model option is ignored for OCX (set it to None). The same config file drives +# both test cases; set them up individually, e.g.: +# compass setup -t landice/ismip7_forcing/atmosphere -w WORKDIR -f ismip7_forcing_ocx_gis.cfg +# compass setup -t landice/ismip7_forcing/ocean_thermal -w WORKDIR -f ismip7_forcing_ocx_gis.cfg + +# config options for ismip7 forcing data +[ismip7] + +# Ice sheet: ais (Antarctic) or gis (Greenland) +ice_sheet = gis + +# Base path to the input ISMIP7 forcing files. User has to supply. +base_path_ismip7 = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/forcing/GIS/OCX + +# Base path to the MALI mesh. User has to supply. +base_path_mali = /global/cfs/cdirs/fanssie/MALI_input_files/GIS_1to10km_r02/ + +# Base path to which output forcing files are saved. +output_base_path = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/test_processing/GIS + +# OCX has no distinct model (RACMO for atmosphere, EN4 for ocean are selected +# automatically), so this option is ignored. +model = None + +# Scenario for forcing data. +scenario = OCX + +# Name of the MALI mesh. Used to name mapping and output files. +mali_mesh_name = GIS_1to10km_r02_20230202 + +# MALI mesh file. User has to supply. +mali_mesh_file = GIS_1to10km_r02_20230202.nc + +# Number of MPI tasks for ESMF_RegridWeightGen +esmf_ntasks = 512 + +# Whether to process time-varying ocean thermal forcing (ESM scenario data) +process_ocean_thermal = true + +# Whether to process observational ocean thermal forcing climatology +process_ocean_climatology = false + +# config options for ismip7 atmosphere forcing +[ismip7_atmosphere] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = conserve + +# Start year for processing +start_year = 1990 + +# End year for processing +end_year = 2025 + +# config options for ismip7 ocean thermal forcing +[ismip7_ocean_thermal] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = bilinear + +# Start year for processing +start_year = 1990 + +# End year for processing +end_year = 2025 + +# config options for ismip7 ocean thermal forcing climatology +[ismip7_ocean_climatology] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = bilinear + +# Base path to observational climatology data +base_path_climatology = None diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py index a1469c8f2d..3e9aea9c87 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -91,11 +91,14 @@ def _run_scenario(self): # Discover input files prefix = params['prefix'] ocean_version = params['ocean_version'] + ocean_grid = params['ocean_grid'] ocean_3d = params['ocean_3d'] + if params['ocean_model'] is not None: + model = params['ocean_model'] input_path = os.path.join(base_path_ismip7, "ocean", "tf", ocean_version) file_pattern = (f"tf_{prefix}_{model}_{scenario}_" - f"ocean_{ocean_version}_*.nc") + f"{ocean_grid}_{ocean_version}_*.nc") all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) if not all_files: diff --git a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst index b535bf081a..0a37029290 100644 --- a/docs/developers_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/developers_guide/landice/test_groups/ismip7_forcing.rst @@ -24,11 +24,19 @@ ice_sheet_params The module :py:mod:`compass.landice.tests.ismip7_forcing.ice_sheet_params` defines a dictionary of ice-sheet-specific parameters (projection, file -naming prefix, grid resolution, data version, ocean dimensionality) and -provides the function +naming prefix, grid resolution, data version, ocean dimensionality, and the +atmosphere/ocean source names) and provides the function :py:func:`compass.landice.tests.ismip7_forcing.ice_sheet_params.get_params` to retrieve them based on the ``ice_sheet`` config option. +When ``scenario = OCX``, ``get_params`` applies a set of OCX overrides on top +of the ice-sheet defaults: data version ``v1``, the ocean file-name grid token +(e.g. ``ocean-1000m``), and the fixed reanalysis sources (``atm_model`` = +``RACMO2.3p2-ERA`` and ``ocean_model`` = ``EN4``). The processing steps use +``atm_model`` / ``ocean_model`` in place of the ``[ismip7] model`` option when +they are set, so the OCX ``model`` option is ignored. This keeps OCX handling +centralized and lets a single config file drive both test cases. + configure ~~~~~~~~~ @@ -39,6 +47,8 @@ the user (i.e., are not ``NotAvailable``). Repository-local example user configs are available at ``compass/landice/tests/ismip7_forcing/ismip7_forcing_test.cfg`` (AIS) and ``compass/landice/tests/ismip7_forcing/ismip7_forcing_test_gis.cfg`` (GrIS). +The GrIS OCX scenario has a dedicated example +``ismip7_forcing_ocx_gis.cfg`` in the same directory. These are intended for development/testing and include environment-specific paths. diff --git a/docs/users_guide/landice/test_groups/ismip7_forcing.rst b/docs/users_guide/landice/test_groups/ismip7_forcing.rst index fdec1eee14..0f9ec60eba 100644 --- a/docs/users_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/users_guide/landice/test_groups/ismip7_forcing.rst @@ -69,6 +69,14 @@ The AIS example enables both ocean processing modes and uses a 2015-2300 processing window for both atmosphere and ocean thermal scenario forcing. The GrIS example enables scenario ocean processing only and uses 1980-2015. +For the GrIS OCX (reanalysis) scenario, use the dedicated example config +``compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_gis.cfg``. OCX has +no distinct ESM model: it uses ``RACMO2.3p2-ERA`` for the atmosphere and +``EN4`` for the ocean, both selected automatically when ``scenario = OCX``. +The ``[ismip7] model`` option is ignored for OCX (set it to ``None``), so a +single config file processes both the ``atmosphere`` and ``ocean_thermal`` +test cases, just like the ESM scenarios. + .. _landice_ismip7_forcing_input_data: Input Data @@ -121,6 +129,20 @@ For GrIS ocean thermal (same 1km grid, 2D, yearly files): ocean/tf/v2/tf_GrIS_{model}_{scenario}_ocean_v2_{year}.nc +The OCX (reanalysis) scenario follows the same directory layout as the ESM +scenarios, but uses fixed reanalysis sources, data version ``v1``, and a +named grid resolution in the ocean file names. For GrIS OCX the atmosphere +source is ``RACMO2.3p2-ERA`` and the ocean source is ``EN4``: + +.. code-block:: none + + acabf/v1/acabf_GrIS_RACMO2.3p2-ERA_OCX_SDBN1-1000m_v1_{year}.nc + ocean/tf/v1/tf_GrIS_EN4_OCX_ocean-1000m_v1_{year}.nc + +Set ``base_path_ismip7`` to the ``OCX`` directory and ``scenario = OCX``. +The sources, version, and ocean grid token are selected automatically for +OCX, and the ``[ismip7] model`` option is ignored (set it to ``None``). + .. _landice_ismip7_forcing_config: config options From 4c7116c215006206450d0a1521f4215caf9b8871 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 3 Sep 2026 10:39:19 -0700 Subject: [PATCH 2/6] Fix ISMIP7 forcing atmosphere/ocean processing bugs Two fixes that block processing of forcing whose remapped files or source directories differ from the CESM test data: - Drop the ncremap-mangled 'time_bnds' variable when combining remapped files. ncremap collides the bounds 'nv' dimension with the destination mesh corner dimension and pads with a fill value that overflows time decoding. The steps build their own xtime and never use time_bnds. - Skip input files whose trailing token is not a year (e.g. climatology averages like *_1978-2007_avg.nc) instead of crashing on int('avg'). --- .../tests/ismip7_forcing/atmosphere/process_runoff.py | 9 +++++++-- .../tests/ismip7_forcing/atmosphere/process_smb.py | 11 ++++++++--- .../ismip7_forcing/atmosphere/process_smb_gradient.py | 9 +++++++-- .../ismip7_forcing/atmosphere/process_temperature.py | 9 +++++++-- .../atmosphere/process_temperature_gradient.py | 9 +++++++-- .../ocean_thermal/process_thermal_forcing.py | 6 ++++-- 6 files changed, 40 insertions(+), 13 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py index ef7ac82164..06852c6690 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -87,7 +87,11 @@ def run(self): # Filter to requested year range input_files = [] for f in all_files: - year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + # skip non-yearly files such as climatology averages (*_avg.nc) + token = os.path.basename(f).split("_")[-1].replace(".nc", "") + if not token.isdigit(): + continue + year = int(token) if start_year <= year <= end_year: input_files.append(f) @@ -178,7 +182,8 @@ def _combine_and_rename(self, remapped_files, output_file): Output file path """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Rename dimensions to MALI conventions rename_dims = {} diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py index 014e098bc9..57444cbb09 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -84,8 +84,12 @@ def run(self): # Filter to requested year range input_files = [] for f in all_files: - # Extract year from filename (last part before .nc) - year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + # Extract year from filename (last part before .nc); skip + # non-yearly files such as climatology averages (e.g. *_avg.nc) + token = os.path.basename(f).split("_")[-1].replace(".nc", "") + if not token.isdigit(): + continue + year = int(token) if start_year <= year <= end_year: input_files.append(f) @@ -166,7 +170,8 @@ def _combine_and_rename(self, remapped_files, output_file): Output file path """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Rename dimensions to MALI conventions rename_dims = {} diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py index 6d39bb304e..e08d11713a 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -85,7 +85,11 @@ def run(self): # Filter to requested year range input_files = [] for f in all_files: - year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + # skip non-yearly files such as climatology averages (*_avg.nc) + token = os.path.basename(f).split("_")[-1].replace(".nc", "") + if not token.isdigit(): + continue + year = int(token) if start_year <= year <= end_year: input_files.append(f) @@ -167,7 +171,8 @@ def _combine_and_rename(self, remapped_files, output_file): Output file path """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Rename dimensions to MALI conventions rename_dims = {} diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index 13ec7c15d9..9e8f0490ba 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -84,7 +84,11 @@ def run(self): # Filter to requested year range input_files = [] for f in all_files: - year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + # skip non-yearly files such as climatology averages (*_avg.nc) + token = os.path.basename(f).split("_")[-1].replace(".nc", "") + if not token.isdigit(): + continue + year = int(token) if start_year <= year <= end_year: input_files.append(f) @@ -166,7 +170,8 @@ def _combine_and_rename(self, remapped_files, output_file): Output file path """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Rename dimensions to MALI conventions rename_dims = {} diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py index 7399eda925..9968101810 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -86,7 +86,11 @@ def run(self): # Filter to requested year range input_files = [] for f in all_files: - year = int(os.path.basename(f).split("_")[-1].replace(".nc", "")) + # skip non-yearly files such as climatology averages (*_avg.nc) + token = os.path.basename(f).split("_")[-1].replace(".nc", "") + if not token.isdigit(): + continue + year = int(token) if start_year <= year <= end_year: input_files.append(f) @@ -168,7 +172,8 @@ def _combine_and_rename(self, remapped_files, output_file): Output file path """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Rename dimensions to MALI conventions rename_dims = {} diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py index 3e9aea9c87..9a67071bc9 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -308,7 +308,8 @@ def _combine_and_rename_3d(self, remapped_files, output_file, Last year to include in output """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Subset to requested year range years = ds.time.dt.year @@ -422,7 +423,8 @@ def _combine_and_rename_2d(self, remapped_files, output_file, Last year to include in output """ ds = xr.open_mfdataset(remapped_files, concat_dim="time", - combine="nested", engine="netcdf4") + combine="nested", engine="netcdf4", + drop_variables="time_bnds") # Subset to requested year range years = ds.time.dt.year From 198473b6a6acea98830982c7d030d4618679079a Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 3 Sep 2026 10:43:02 -0700 Subject: [PATCH 3/6] Write ISMIP7 forcing output in the layout ismip7_run ingests Restructure processed-forcing output from {output_base_path}/atmosphere_forcing/{model}_{scenario}/ (and the analogous ocean_thermal_forcing path) to {output_base_path}/{group}/atmosphere/ and {output_base_path}/{group}/ocean_thermal_forcing/, matching the directory layout the ismip7_run test group globs for. 'group' is {model}_{scenario} for the ESM scenarios and {scenario} (i.e. OCX) for OCX, whose atmosphere and ocean sources differ but must share one directory. Documents the layout and how to point ismip7_run's forcing_basepath / ocx_forcing_path at it. --- .../atmosphere/process_runoff.py | 7 ++++-- .../ismip7_forcing/atmosphere/process_smb.py | 7 ++++-- .../atmosphere/process_smb_gradient.py | 7 ++++-- .../atmosphere/process_temperature.py | 7 ++++-- .../process_temperature_gradient.py | 7 ++++-- .../ocean_thermal/process_thermal_forcing.py | 7 ++++-- .../landice/test_groups/ismip7_forcing.rst | 25 +++++++++++++++++++ 7 files changed, 55 insertions(+), 12 deletions(-) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py index 06852c6690..3f0e7b9e5f 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -73,7 +73,10 @@ def run(self): resolution = params['atm_resolution'] version = params['atm_version'] if params['atm_model'] is not None: + forcing_group = scenario model = params['atm_model'] + else: + forcing_group = f"{model}_{scenario}" input_path = os.path.join(base_path_ismip7, "mrro", version) file_pattern = (f"mrro_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") @@ -158,8 +161,8 @@ def run(self): os.remove(f) # Place output in appropriate directory - output_path = os.path.join(output_base_path, "atmosphere_forcing", - f"{model}_{scenario}") + output_path = os.path.join(output_base_path, forcing_group, + "atmosphere") if not os.path.exists(output_path): os.makedirs(output_path) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py index 57444cbb09..04de34180f 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -70,7 +70,10 @@ def run(self): resolution = params['atm_resolution'] version = params['atm_version'] if params['atm_model'] is not None: + forcing_group = scenario model = params['atm_model'] + else: + forcing_group = f"{model}_{scenario}" input_path = os.path.join(base_path_ismip7, "acabf", version) file_pattern = (f"acabf_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") @@ -146,8 +149,8 @@ def run(self): os.remove(f) # Place output in appropriate directory - output_path = os.path.join(output_base_path, "atmosphere_forcing", - f"{model}_{scenario}") + output_path = os.path.join(output_base_path, forcing_group, + "atmosphere") if not os.path.exists(output_path): os.makedirs(output_path) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py index e08d11713a..d6c291e4d5 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -71,7 +71,10 @@ def run(self): resolution = params['atm_resolution'] version = params['atm_version'] if params['atm_model'] is not None: + forcing_group = scenario model = params['atm_model'] + else: + forcing_group = f"{model}_{scenario}" input_path = os.path.join(base_path_ismip7, "dacabfdz", version) file_pattern = (f"dacabfdz_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") @@ -147,8 +150,8 @@ def run(self): os.remove(f) # Place output in appropriate directory - output_path = os.path.join(output_base_path, "atmosphere_forcing", - f"{model}_{scenario}") + output_path = os.path.join(output_base_path, forcing_group, + "atmosphere") if not os.path.exists(output_path): os.makedirs(output_path) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index 9e8f0490ba..3168c69f28 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -70,7 +70,10 @@ def run(self): resolution = params['atm_resolution'] version = params['atm_version'] if params['atm_model'] is not None: + forcing_group = scenario model = params['atm_model'] + else: + forcing_group = f"{model}_{scenario}" input_path = os.path.join(base_path_ismip7, "ts", version) file_pattern = (f"ts_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") @@ -146,8 +149,8 @@ def run(self): os.remove(f) # Place output in appropriate directory - output_path = os.path.join(output_base_path, "atmosphere_forcing", - f"{model}_{scenario}") + output_path = os.path.join(output_base_path, forcing_group, + "atmosphere") if not os.path.exists(output_path): os.makedirs(output_path) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py index 9968101810..16181f02e5 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -72,7 +72,10 @@ def run(self): resolution = params['atm_resolution'] version = params['atm_version'] if params['atm_model'] is not None: + forcing_group = scenario model = params['atm_model'] + else: + forcing_group = f"{model}_{scenario}" input_path = os.path.join(base_path_ismip7, "dtsdz", version) file_pattern = (f"dtsdz_{prefix}_{model}_{scenario}_" f"SDBN1-{resolution}_{version}_*.nc") @@ -148,8 +151,8 @@ def run(self): os.remove(f) # Place output in appropriate directory - output_path = os.path.join(output_base_path, "atmosphere_forcing", - f"{model}_{scenario}") + output_path = os.path.join(output_base_path, forcing_group, + "atmosphere") if not os.path.exists(output_path): os.makedirs(output_path) diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py index 9a67071bc9..fba0053b1d 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -94,7 +94,10 @@ def _run_scenario(self): ocean_grid = params['ocean_grid'] ocean_3d = params['ocean_3d'] if params['ocean_model'] is not None: + forcing_group = scenario model = params['ocean_model'] + else: + forcing_group = f"{model}_{scenario}" input_path = os.path.join(base_path_ismip7, "ocean", "tf", ocean_version) file_pattern = (f"tf_{prefix}_{model}_{scenario}_" @@ -187,8 +190,8 @@ def _run_scenario(self): os.remove(f) # Place output in appropriate directory - output_path = os.path.join(output_base_path, "ocean_thermal_forcing", - f"{model}_{scenario}") + output_path = os.path.join(output_base_path, forcing_group, + "ocean_thermal_forcing") if not os.path.exists(output_path): os.makedirs(output_path) diff --git a/docs/users_guide/landice/test_groups/ismip7_forcing.rst b/docs/users_guide/landice/test_groups/ismip7_forcing.rst index 0f9ec60eba..6afb6b7d92 100644 --- a/docs/users_guide/landice/test_groups/ismip7_forcing.rst +++ b/docs/users_guide/landice/test_groups/ismip7_forcing.rst @@ -77,6 +77,31 @@ The ``[ismip7] model`` option is ignored for OCX (set it to ``None``), so a single config file processes both the ``atmosphere`` and ``ocean_thermal`` test cases, just like the ESM scenarios. +.. _landice_ismip7_forcing_output: + +Output Layout +------------- + +Processed forcing is written under ``output_base_path`` in a layout that the +:ref:`landice_ismip7_run` test group can ingest directly: + +.. code-block:: none + + {output_base_path}/{group}/atmosphere/{mesh}_SMB_{source}_{scenario}_{years}.nc + {output_base_path}/{group}/atmosphere/{mesh}_temperature_{source}_{scenario}_{years}.nc + {output_base_path}/{group}/atmosphere/{mesh}_runoff_... (and the two gradients) + {output_base_path}/{group}/ocean_thermal_forcing/{mesh}_thermal_forcing_{source}_{scenario}_{years}.nc + +The ``group`` directory is ``{model}_{scenario}`` for the ESM scenarios and +``{scenario}`` (i.e. ``OCX``) for OCX, whose atmosphere and ocean use +different sources (``RACMO2.3p2-ERA`` and ``EN4``) but must share one +directory. To feed :ref:`landice_ismip7_run`, point its ``forcing_basepath`` +at ``output_base_path`` for ESM scenarios, or its ``ocx_forcing_path`` at +``{output_base_path}/OCX`` for OCX. + +Process only one year range into a given ``group`` directory: the run setup +expects exactly one file per forcing field there. + .. _landice_ismip7_forcing_input_data: Input Data From 828699bf0d987c1738ef1286103516df9b7de1a4 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Thu, 3 Sep 2026 11:55:35 -0700 Subject: [PATCH 4/6] Fix missing-value leakage in ISMIP7 atmosphere forcing; share extrapolation helper process_smb.py, process_temperature.py, process_smb_gradient.py, and process_temperature_gradient.py remapped the raw source file directly with ncremap, unlike process_runoff.py, ocean_thermal, and the fracture pathway, which all extrapolate fill/missing values on the source grid first. Without that step, no-data cells near the ice-sheet margin (e.g. ocean cells in the RACMO OCX source) leak the ~9.97e36 netCDF fill sentinel into valid neighboring cells during conservative/bilinear remapping. Add the same pre-remap extrapolation to the four affected steps. Also promote the previously-duplicated extrapolate_source() (from fracture/remap_utils.py, formerly re-implemented again in ocean_thermal/process_thermal_forcing.py and process_runoff.py) to a single shared compass.landice.tests.ismip7_forcing.remap_utils module, and update all five call sites to use it instead of per-file private copies. --- .../atmosphere/process_runoff.py | 59 +---------------- .../ismip7_forcing/atmosphere/process_smb.py | 12 +++- .../atmosphere/process_smb_gradient.py | 13 +++- .../atmosphere/process_temperature.py | 12 +++- .../process_temperature_gradient.py | 12 +++- .../fracture/process_excess_melt.py | 2 +- .../fracture/process_lake_properties.py | 2 +- .../ismip7_forcing/fracture/remap_utils.py | 59 ----------------- .../ocean_thermal/process_thermal_forcing.py | 65 +------------------ .../tests/ismip7_forcing/remap_utils.py | 64 ++++++++++++++++++ 10 files changed, 116 insertions(+), 184 deletions(-) create mode 100644 compass/landice/tests/ismip7_forcing/remap_utils.py diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py index 3f0e7b9e5f..801844edf0 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_runoff.py @@ -2,16 +2,15 @@ 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.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step @@ -132,8 +131,7 @@ def run(self): # 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) + extrapolate_source(input_file, extrap_file, "mrro", logger) logger.info(f" Remapping: {basename}") args = ["ncremap", @@ -235,56 +233,3 @@ def _combine_and_rename(self, remapped_files, output_file): ds = ds.drop_vars("Time") 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) diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py index 04de34180f..3bb864710e 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb.py @@ -10,6 +10,7 @@ build_mapping_file, ) from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step @@ -126,15 +127,24 @@ def run(self): 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): + extrapolate_source(input_file, extrap_file, "acabf", logger) + logger.info(f" Remapping: {basename}") args = ["ncremap", - "-i", input_file, + "-i", extrap_file, "-o", remapped_file, "-m", mapping_file, "-v", "acabf"] 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}_SMB_{model}_{scenario}_" diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py index d6c291e4d5..d44192d14c 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_smb_gradient.py @@ -10,6 +10,7 @@ build_mapping_file, ) from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step @@ -127,15 +128,25 @@ def run(self): 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): + extrapolate_source(input_file, extrap_file, "dacabfdz", + logger) + logger.info(f" Remapping: {basename}") args = ["ncremap", - "-i", input_file, + "-i", extrap_file, "-o", remapped_file, "-m", mapping_file, "-v", "dacabfdz"] 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}_SMB_gradient_{model}_{scenario}_" diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py index 3168c69f28..24c9529342 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature.py @@ -10,6 +10,7 @@ build_mapping_file, ) from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step @@ -126,15 +127,24 @@ def run(self): 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): + extrapolate_source(input_file, extrap_file, "ts", logger) + logger.info(f" Remapping: {basename}") args = ["ncremap", - "-i", input_file, + "-i", extrap_file, "-o", remapped_file, "-m", mapping_file, "-v", "ts"] 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}_temperature_{model}_{scenario}_" diff --git a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py index 16181f02e5..45d39b9a1d 100644 --- a/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py +++ b/compass/landice/tests/ismip7_forcing/atmosphere/process_temperature_gradient.py @@ -10,6 +10,7 @@ build_mapping_file, ) from compass.landice.tests.ismip7_forcing.ice_sheet_params import get_params +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step @@ -128,15 +129,24 @@ def run(self): 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): + extrapolate_source(input_file, extrap_file, "dtsdz", logger) + logger.info(f" Remapping: {basename}") args = ["ncremap", - "-i", input_file, + "-i", extrap_file, "-o", remapped_file, "-m", mapping_file, "-v", "dtsdz"] 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}_temperature_gradient_{model}_" diff --git a/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py b/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py index fbf3f75b83..8a9275c962 100644 --- a/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py +++ b/compass/landice/tests/ismip7_forcing/fracture/process_excess_melt.py @@ -12,9 +12,9 @@ ) from compass.landice.tests.ismip7_forcing.fracture.remap_utils import ( add_xtime_and_write, - extrapolate_source, open_rename_and_trim, ) +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step diff --git a/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py b/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py index 10df254c69..6867f20735 100644 --- a/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py +++ b/compass/landice/tests/ismip7_forcing/fracture/process_lake_properties.py @@ -9,9 +9,9 @@ ) from compass.landice.tests.ismip7_forcing.fracture.remap_utils import ( add_xtime_and_write, - extrapolate_source, open_rename_and_trim, ) +from compass.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step diff --git a/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py b/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py index 8da5de4ae4..14518c7c79 100644 --- a/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py +++ b/compass/landice/tests/ismip7_forcing/fracture/remap_utils.py @@ -1,67 +1,8 @@ """ Shared helpers for remapping ISMIP7 fracture forcing data to the MALI mesh. """ -import os - -import numpy as np import xarray as xr from mpas_tools.io import write_netcdf -from scipy.ndimage import distance_transform_edt - - -def extrapolate_source(input_file, output_file, varnames, 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 - - varnames : str or list of str - Name(s) of the variable(s) to extrapolate - - logger : logging.Logger - Logger for status messages - """ - if isinstance(varnames, str): - varnames = [varnames] - - logger.info(f" Extrapolating fill values on source grid: " - f"{os.path.basename(input_file)}") - - ds = xr.open_dataset(input_file, decode_times=False) - - for varname in varnames: - data = ds[varname] - values = data.values.copy() - non_spatial_shape = values.shape[:-2] - - for idx in np.ndindex(non_spatial_shape): - slab = values[idx] - 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 - if "_FillValue" in ds[varname].encoding: - del ds[varname].encoding["_FillValue"] - - write_netcdf(ds, output_file) - ds.close() def open_rename_and_trim(remapped_file, rename_vars, start_year, end_year): diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py index fba0053b1d..c0eed44c97 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -2,16 +2,15 @@ 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.landice.tests.ismip7_forcing.remap_utils import extrapolate_source from compass.step import Step @@ -156,8 +155,7 @@ def _run_scenario(self): # 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, "tf", - logger) + extrapolate_source(input_file, extrap_file, "tf", logger) logger.info(f" Remapping: {basename}") args = ["ncremap", @@ -252,8 +250,7 @@ def _run_climatology(self): if not os.path.exists(remapped_file): extrap_file = f"extrap_{basename}" if not os.path.exists(extrap_file): - self._extrapolate_source(input_file, extrap_file, "tf", - logger) + extrapolate_source(input_file, extrap_file, "tf", logger) logger.info(f" Remapping: {basename}") args = ["ncremap", @@ -571,59 +568,3 @@ def _rename_climatology_3d(self, remapped_file, output_file): ds = ds.drop_vars("nISMIP6OceanLayers") 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 interpolation from valid cells. 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., "tf") - - 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 (and z level if 3D) - # Source files have dims like (time, z, y, x) or (time, y, x) - values = data.values.copy() - non_spatial_shape = values.shape[:-2] # (time,) or (time, z) - - # Use distance_transform_edt with return_indices to find the - # nearest valid cell index for each invalid cell. This is O(n) - # on the grid and much faster than KD-tree approaches. - 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) diff --git a/compass/landice/tests/ismip7_forcing/remap_utils.py b/compass/landice/tests/ismip7_forcing/remap_utils.py new file mode 100644 index 0000000000..e54f9e9fbf --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/remap_utils.py @@ -0,0 +1,64 @@ +""" +Shared helpers for remapping ISMIP7 forcing data to the MALI mesh. +""" +import os + +import numpy as np +import xarray as xr +from mpas_tools.io import write_netcdf +from scipy.ndimage import distance_transform_edt + + +def extrapolate_source(input_file, output_file, varnames, 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 + + varnames : str or list of str + Name(s) of the variable(s) to extrapolate + + logger : logging.Logger + Logger for status messages + """ + if isinstance(varnames, str): + varnames = [varnames] + + logger.info(f" Extrapolating fill values on source grid: " + f"{os.path.basename(input_file)}") + + ds = xr.open_dataset(input_file, decode_times=False) + + for varname in varnames: + data = ds[varname] + values = data.values.copy() + non_spatial_shape = values.shape[:-2] + + for idx in np.ndindex(non_spatial_shape): + slab = values[idx] + 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 + if "_FillValue" in ds[varname].encoding: + del ds[varname].encoding["_FillValue"] + + write_netcdf(ds, output_file) + ds.close() From 7047302c97b272fed33039a3aa7d32973f4380f7 Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Fri, 4 Sep 2026 07:54:58 -0700 Subject: [PATCH 5/6] Fix stale autosummary reference after remap_utils refactor extrapolate_source moved from fracture.remap_utils to the shared compass.landice.tests.ismip7_forcing.remap_utils module (828699bf0). Update the api.rst autosummary entry accordingly; this was breaking the strict Sphinx build (autosummary failed to import the moved symbol). --- docs/developers_guide/landice/api.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/developers_guide/landice/api.rst b/docs/developers_guide/landice/api.rst index d2cd995145..75bb6bab8c 100644 --- a/docs/developers_guide/landice/api.rst +++ b/docs/developers_guide/landice/api.rst @@ -385,6 +385,7 @@ ismip7_forcing configure.configure ice_sheet_params.get_params create_mapfile.build_mapping_file + remap_utils.extrapolate_source atmosphere.Atmosphere atmosphere.Atmosphere.configure @@ -412,7 +413,6 @@ ismip7_forcing fracture.Fracture fracture.Fracture.configure - fracture.remap_utils.extrapolate_source fracture.remap_utils.open_rename_and_trim fracture.remap_utils.add_xtime_and_write fracture.process_excess_melt.ProcessExcessMelt From ec4fdaf6b0e4e30fe3e74a83380360db3e20c05b Mon Sep 17 00:00:00 2001 From: Trevor Hillebrand Date: Wed, 9 Sep 2026 13:09:53 -0700 Subject: [PATCH 6/6] Add AIS support for OCX scenario in ismip7_forcing Enable the OCX (reanalysis) scenario for the Antarctic ice sheet: - ice_sheet_params: add an 'ais' OCX override (RACMO atmosphere at v1, reanalysis ocean at v1 with no model token) and an ocean_choice_layout flag for the idiosyncratic AIS OCX ocean layout. - process_thermal_forcing: assemble a list of ocean forcing sources so AIS OCX can process one or more ocean choices (main/cold/warm/vary), each in its own OCX_ output directory; other cases are unchanged. - ismip7_forcing.cfg: add ocean_choice option (AIS OCX only). - remap_utils: write the extrapolated source file as CDF-5 so large 3D ocean thermal forcing fields aren't limited by the NetCDF3 per-variable size cap or HDF5 chunk limits that break ncremap. - add ismip7_forcing_ocx_ais.cfg example config. --- .../tests/ismip7_forcing/ice_sheet_params.py | 9 ++ .../tests/ismip7_forcing/ismip7_forcing.cfg | 4 + .../ismip7_forcing/ismip7_forcing_ocx_ais.cfg | 84 +++++++++++ .../ocean_thermal/process_thermal_forcing.py | 136 +++++++++++++++--- .../tests/ismip7_forcing/remap_utils.py | 19 ++- 5 files changed, 232 insertions(+), 20 deletions(-) create mode 100644 compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_ais.cfg diff --git a/compass/landice/tests/ismip7_forcing/ice_sheet_params.py b/compass/landice/tests/ismip7_forcing/ice_sheet_params.py index 0d57094e13..d991eae45c 100644 --- a/compass/landice/tests/ismip7_forcing/ice_sheet_params.py +++ b/compass/landice/tests/ismip7_forcing/ice_sheet_params.py @@ -42,6 +42,15 @@ 'atm_model': 'RACMO2.3p2-ERA', 'ocean_model': 'EN4', }, + 'ais': { + 'atm_version': 'v1', + 'ocean_version': 'v1', + 'atm_model': 'RACMO2.3p2-ERA', + 'ocean_model': None, + # AIS OCX ocean files have no model token and live in per-choice + # subdirectories (main/cold/warm/vary); see process_thermal_forcing. + 'ocean_choice_layout': True, + }, } diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg index f02b79067a..be7d195aa4 100644 --- a/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing.cfg @@ -54,6 +54,10 @@ end_year = 2014 # Remapping method. Options: bilinear, neareststod, conserve method_remap = bilinear +# Ocean forcing choice(s), only used for the AIS OCX scenario. Comma-separated +# list of any of: main, cold, warm, vary. Use 'all' to process every choice. +ocean_choice = main + # Start year for processing start_year = 1850 diff --git a/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_ais.cfg b/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_ais.cfg new file mode 100644 index 0000000000..975cb89b85 --- /dev/null +++ b/compass/landice/tests/ismip7_forcing/ismip7_forcing_ocx_ais.cfg @@ -0,0 +1,84 @@ +# Example config for processing the AIS ISMIP7 OCX (reanalysis) forcing. +# OCX has no distinct ESM model: it uses RACMO2.3p2-ERA for the atmosphere and +# reanalysis ocean products, selected automatically when scenario = OCX. The +# [ismip7] model option is ignored for OCX (set it to None). Unlike GrIS, the +# AIS OCX ocean provides several forcing choices (main/cold/warm/vary) selected +# with the ocean_choice option below. The same config file drives both test +# cases; set them up individually, e.g.: +# compass setup -t landice/ismip7_forcing/atmosphere -w WORKDIR -f ismip7_forcing_ocx_ais.cfg +# compass setup -t landice/ismip7_forcing/ocean_thermal -w WORKDIR -f ismip7_forcing_ocx_ais.cfg + +# config options for ismip7 forcing data +[ismip7] + +# Ice sheet: ais (Antarctic) or gis (Greenland) +ice_sheet = ais + +# Base path to the input ISMIP7 forcing files. User has to supply. +base_path_ismip7 = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/forcing/AIS/OCX + +# Base path to the MALI mesh. User has to supply. +base_path_mali = /global/cfs/cdirs/fanssie/MALI_input_files/AIS_4to20km_r01 + +# Base path to which output forcing files are saved. +output_base_path = /global/cfs/cdirs/m4288/users/trhille/ISMIP7/test_processing/AIS + +# OCX has no distinct model (RACMO for atmosphere, reanalysis for ocean are +# selected automatically), so this option is ignored. +model = None + +# Scenario for forcing data. +scenario = OCX + +# Name of the MALI mesh. Used to name mapping and output files. +mali_mesh_name = AIS_4to20km_r01_20220907 + +# MALI mesh file. User has to supply. +mali_mesh_file = AIS_4to20km_r01_20220907.nc + +# Number of MPI tasks for ESMF_RegridWeightGen +esmf_ntasks = 512 + +# Whether to process time-varying ocean thermal forcing (ESM scenario data) +process_ocean_thermal = true + +# Whether to process observational ocean thermal forcing climatology +process_ocean_climatology = false + +# config options for ismip7 atmosphere forcing +[ismip7_atmosphere] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = conserve + +# Start year for processing +start_year = 1979 + +# End year for processing +end_year = 2025 + +# config options for ismip7 ocean thermal forcing +[ismip7_ocean_thermal] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = bilinear + +# Ocean forcing choice(s) for the AIS OCX scenario. Comma-separated list of any +# of: main, cold, warm, vary. Use 'all' to process every choice. Each choice is +# written to its own OCX_ output directory. +ocean_choice = main + +# Start year for processing +start_year = 1950 + +# End year for processing +end_year = 2025 + +# config options for ismip7 ocean thermal forcing climatology +[ismip7_ocean_climatology] + +# Remapping method. Options: bilinear, neareststod, conserve +method_remap = bilinear + +# Base path to observational climatology data +base_path_climatology = None diff --git a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py index c0eed44c97..cc71341cbb 100644 --- a/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py +++ b/compass/landice/tests/ismip7_forcing/ocean_thermal/process_thermal_forcing.py @@ -67,9 +67,9 @@ def run(self): def _run_scenario(self): """ Process time-varying ocean thermal forcing from an ESM - (e.g., CESM2-WACCM historical or ssp585). + (e.g., CESM2-WACCM historical or ssp585) or, for AIS OCX, from one + or more reanalysis ocean "choices" (main/cold/warm/vary). """ - logger = self.logger config = self.config params = get_params(config) @@ -87,20 +87,123 @@ def _run_scenario(self): start_year = section.getint("start_year") end_year = section.getint("end_year") - # Discover input files prefix = params['prefix'] ocean_version = params['ocean_version'] ocean_grid = params['ocean_grid'] ocean_3d = params['ocean_3d'] - if params['ocean_model'] is not None: - forcing_group = scenario - model = params['ocean_model'] + + # Assemble the list of forcing sources to process. AIS OCX has several + # ocean "choices" (main/cold/warm/vary) in per-choice subdirectories + # with no model token in the filename; all other cases (ESM scenarios + # and GrIS OCX) have a single source. + jobs = [] + if params.get('ocean_choice_layout', False): + for choice in self._get_ocean_choices(config): + input_path = os.path.join(base_path_ismip7, "ocean", choice, + ocean_version) + file_pattern = (f"tf_{prefix}_{scenario}_ocean_{choice}_" + f"{ocean_version}_*.nc") + jobs.append({ + 'input_path': input_path, + 'file_pattern': file_pattern, + 'forcing_group': f"{scenario}_{choice}", + 'label': f"{scenario}_{choice}", + }) else: - forcing_group = f"{model}_{scenario}" - input_path = os.path.join(base_path_ismip7, "ocean", "tf", - ocean_version) - file_pattern = (f"tf_{prefix}_{model}_{scenario}_" - f"{ocean_grid}_{ocean_version}_*.nc") + if params['ocean_model'] is not None: + forcing_group = scenario + model = params['ocean_model'] + else: + forcing_group = f"{model}_{scenario}" + input_path = os.path.join(base_path_ismip7, "ocean", "tf", + ocean_version) + file_pattern = (f"tf_{prefix}_{model}_{scenario}_" + f"{ocean_grid}_{ocean_version}_*.nc") + jobs.append({ + 'input_path': input_path, + 'file_pattern': file_pattern, + 'forcing_group': forcing_group, + 'label': f"{model}_{scenario}", + }) + + # Mapping file is shared across sources (identical ocean source grid). + mapping_file = (f"map_ismip7_{ice_sheet}_ocean_to_" + f"{mali_mesh_name}_{method_remap}.nc") + + for job in jobs: + self._process_ocean_forcing( + job, mapping_file, ocean_3d, method_remap, + start_year, end_year, mali_mesh_name, mali_mesh_file, + output_base_path) + + def _get_ocean_choices(self, config): + """ + Parse the comma-separated ``ocean_choice`` option (AIS OCX only) into a + list of ocean choices, expanding ``all`` to every available choice. + + Parameters + ---------- + config : compass.config.CompassConfigParser + Configuration options for the test case + + Returns + ------- + choices : list of str + Ordered, de-duplicated list of ocean choices to process + """ + valid = ["main", "cold", "warm", "vary"] + raw = config.get("ismip7_ocean_thermal", "ocean_choice") + tokens = [t.strip() for t in raw.split(",") if t.strip()] + if any(t.lower() == "all" for t in tokens): + return valid + + choices = [] + for token in tokens: + if token not in valid: + raise ValueError( + f"Invalid ocean_choice '{token}'. Must be a " + f"comma-separated list of {valid} or 'all'.") + if token not in choices: + choices.append(token) + if not choices: + raise ValueError( + "No ocean_choice specified for AIS OCX ocean forcing.") + return choices + + def _process_ocean_forcing(self, job, mapping_file, ocean_3d, + method_remap, start_year, end_year, + mali_mesh_name, mali_mesh_file, + output_base_path): + """ + Discover, remap, combine, and save the thermal forcing for a single + forcing source described by ``job``. + + Parameters + ---------- + job : dict + Source description with keys ``input_path``, ``file_pattern``, + ``forcing_group``, and ``label`` + mapping_file : str + Path of the shared ocean-to-MALI mapping file (built on demand) + ocean_3d : bool + Whether the thermal forcing is 3D (AIS) or 2D (GrIS) + method_remap : str + Remapping method passed to the mapping-file builder + start_year, end_year : int + Inclusive year range to process + mali_mesh_name, mali_mesh_file : str + MALI mesh name and file + output_base_path : str + Base path under which output is written + """ + logger = self.logger + config = self.config + + input_path = job['input_path'] + file_pattern = job['file_pattern'] + forcing_group = job['forcing_group'] + label = job['label'] + all_files = sorted(glob.glob(os.path.join(input_path, file_pattern))) if not all_files: @@ -109,8 +212,8 @@ def _run_scenario(self): f" {os.path.join(input_path, file_pattern)}") # Filter to files that overlap with the requested year range. - # AIS files are named with decade ranges (e.g., 1850-1859). - # GrIS files are named with single years (e.g., 2015). + # AIS files are named with decade or multi-decade ranges (e.g., + # 1850-1859, 1950-2025). GrIS files are named with single years. input_files = [] for f in all_files: # Extract year range from filename (last part before .nc) @@ -130,9 +233,6 @@ def _run_scenario(self): f"overlapping years {start_year}-{end_year}") # Build mapping file using the first input file as grid template. - mapping_file = (f"map_ismip7_{ice_sheet}_ocean_to_" - f"{mali_mesh_name}_{method_remap}.nc") - if not os.path.exists(mapping_file): logger.info("Building mapping file for ocean grid...") build_mapping_file(config, logger, @@ -140,7 +240,7 @@ def _run_scenario(self): mali_mesh_file=mali_mesh_file, method_remap=method_remap) - # Remap each decade file + # Remap each file remapped_files = [] for input_file in input_files: basename = os.path.basename(input_file) @@ -171,7 +271,7 @@ def _run_scenario(self): # Combine remapped files and rename to MALI conventions logger.info("Combining remapped files and renaming variables...") - output_file = (f"{mali_mesh_name}_thermal_forcing_{model}_{scenario}_" + output_file = (f"{mali_mesh_name}_thermal_forcing_{label}_" f"{start_year}-{end_year}.nc") if ocean_3d: diff --git a/compass/landice/tests/ismip7_forcing/remap_utils.py b/compass/landice/tests/ismip7_forcing/remap_utils.py index e54f9e9fbf..890a979ab1 100644 --- a/compass/landice/tests/ismip7_forcing/remap_utils.py +++ b/compass/landice/tests/ismip7_forcing/remap_utils.py @@ -3,9 +3,9 @@ """ import os +import netCDF4 import numpy as np import xarray as xr -from mpas_tools.io import write_netcdf from scipy.ndimage import distance_transform_edt @@ -60,5 +60,20 @@ def extrapolate_source(input_file, output_file, varnames, logger): if "_FillValue" in ds[varname].encoding: del ds[varname].encoding["_FillValue"] - write_netcdf(ds, output_file) + # Preserve a fill value for slabs that remain fully invalid after + # extrapolation so ncremap ignores them during interpolation. + encoding = {} + for varname in varnames: + dtype = ds[varname].dtype + if np.issubdtype(dtype, np.floating) and \ + bool(np.any(np.isnan(ds[varname].values))): + fill = netCDF4.default_fillvals[dtype.str[1:]] + encoding[varname] = {"_FillValue": dtype.type(fill)} + + # Write CDF-5 (NETCDF3_64BIT_DATA): a classic-model format with 64-bit + # sizes that supports very large variables (e.g. AIS 3D ocean thermal + # forcing) without the HDF5 chunk-size limits that make ncremap unable to + # open large NETCDF4 files. This is a temporary file consumed by ncremap. + ds.to_netcdf(output_file, format="NETCDF3_64BIT_DATA", engine="netcdf4", + encoding=encoding) ds.close()