diff --git a/landice/output_processing_li/ismip7_postprocessing/grid_and_mapping.py b/landice/output_processing_li/ismip7_postprocessing/grid_and_mapping.py new file mode 100644 index 000000000..a47ec862d --- /dev/null +++ b/landice/output_processing_li/ismip7_postprocessing/grid_and_mapping.py @@ -0,0 +1,213 @@ +from mpas_tools.scrip.from_mpas import scrip_from_mpas +from subprocess import check_call +import os +import netCDF4 + + +VALID_EXPERIMENTS = [f"C{i:03d}" for i in range(1, 12)] +VALID_RESOLUTIONS = [1, 2, 4, 8, 16] + + +def check_res(res): + """ + Validate the ISMIP7 grid resolution. + + Parameters + ---------- + res : str or int + Resolution in kilometres to validate. + + Raises + ------ + ValueError + If the resolution is not in the list of valid values. + """ + try: + res_int = int(res) + except (ValueError, TypeError): + raise ValueError( + f"Resolution '{res}' is not a valid integer. " + f"Valid resolutions (km) are: {VALID_RESOLUTIONS}") + if res_int not in VALID_RESOLUTIONS: + raise ValueError( + f"Invalid resolution '{res_int}' km. " + f"Valid resolutions (km) are: {VALID_RESOLUTIONS}") + print(f"Resolution {res_int} km is valid.") + + +def check_exp_name(exp): + """ + Validate the ISMIP7 experiment name. + + Parameters + ---------- + exp : str + Experiment name to validate (e.g. 'C001'). + + Raises + ------ + ValueError + If the experiment name is not in the list of valid experiments. + """ + if exp not in VALID_EXPERIMENTS: + raise ValueError( + f"Invalid experiment name '{exp}'. " + f"Valid experiments are: {', '.join(VALID_EXPERIMENTS)}") + print(f"Experiment name '{exp}' is valid.") + + +def check_ismip7_grid_file(ismip7_grid_file_path, res_ismip7_grid): + """ + Ensure the ISMIP7 grid file has 'x' and 'y' coordinate variables and that + its corners match the ISMIP7-required extents. + + Parameters + ---------- + ismip7_grid_file_path : str + Path to the ISMIP7 grid file supplied by the user. + res_ismip7_grid : str + Resolution of the ISMIP7 grid in kilometres (e.g. '8'). + + Returns + ------- + None + """ + print("\n---Checking the coordinate variables of the ismip7 grid file---") + data_ismip7 = netCDF4.Dataset(ismip7_grid_file_path, "r") + + if 'x' in data_ismip7.variables and 'y' in data_ismip7.variables: + ismip7_grid_file = ismip7_grid_file_path + print("'x' and 'y' coordinates exist in the file.") + else: + data_ismip7.close() + raise ValueError( + f"'x' and/or 'y' coordinate variables are missing from the " + f"ISMIP7 grid file '{ismip7_grid_file_path}'. Please provide a " + f"grid file that includes 'x' and 'y' coordinate variables.") + + data_ismip7.close() + + print("Checking the grid corners...") + check_ds = netCDF4.Dataset(ismip7_grid_file, "r") + x = check_ds.variables["x"] + y = check_ds.variables["y"] + if not x[0] == -3040000 or not y[0] == -3040000: + raise ValueError( + f"The lower left corner values must be at " + f"-3040000m and -3040000m. But the values are at " + f"{x[0]}m and {y[0]}m. Check the value you " + f"provided for '--res' matches with the resolution of " + f"the MALI output files. ") + elif not x[-1] == 3040000 or not y[-1] == 3040000: + raise ValueError( + f"The upper right corner values must be at " + f"3040000m and 3040000m. But the values are at " + f"{x[-1]}m and {y[-1]}m. Check the value you " + f"provided for '--res' matches with the resolution of " + f"the MALI output files. ") + else: + print(f"Grid corners are as ismip7-required: " + f"lower left corner values at {x[0]}m and {y[0]}m, and " + f"upper right corner values at {x[-1]}m and {y[-1]}m") + check_ds.close() + + +def build_mapping_file(mali_mesh_file, + mapping_file, res_ismip7_grid, + ismip7_grid_file=None, + method_remap=None): + """ + Build a mapping file if it does not exist. + Mapping file is then used to remap the MALI source file to the + ISMIP7 ppolarstero grid + + Parameters + ---------- + + mali_mesh_file : str + mali file + + mapping_file : str + weights for interpolation from mali_mesh_file to ismip7_grid_file + + res_ismip7_grid: str + resolution of the ismip7 grid in kilometers + + ismip7_grid_file : str, optional + The ISMIP7 file if mapping file does not exist + + method_remap : str, optional + Remapping method used in building a mapping file + """ + + if os.path.exists(mapping_file): + print("Mapping file exists. Not building a new one.") + return + + if ismip7_grid_file is None: + raise ValueError("Mapping file does not exist. To build one, ISMIP7 " + "grid file with '-f' should be provided. " + "Type --help for info") + + if method_remap is None: + method_remap = "conserve" + + ismip7_scripfile = f"temp_ismip7_{res_ismip7_grid}km_scrip.nc" + mali_scripfile = "temp_mali_scrip.nc" + ismip7_projection = "ais-bedmap2" + + # create the ismip7 scripfile if mapping file does not exist + # this is the projection of ismip7 data for Antarctica + print("Mapping file does not exist. Building one based on " + "the input/ouptut meshes") + print("Creating temporary scripfiles " + "for ismip7 grid and mali mesh...") + + # create a scripfile for ismip7 grid + args = ["create_scrip_file_from_planar_rectangular_grid", + "--input", ismip7_grid_file, + "--scrip", ismip7_scripfile, + "--proj", ismip7_projection, + "--rank", "2"] + + check_call(args) + + # create a MALI mesh scripfile + scrip_from_mpas(mali_mesh_file, mali_scripfile) + + # create a mapping file using ESMF weight gen + print(f"Creating a mapping file... " + f"Mapping method used: {method_remap}") + + # generate a mapping file using the scrip files + + # On compute node, ESMG regridder needs to be called with srun + # On head nodes or local machines it does not. + # Here, assuming a compute node has hostname starting with nid, + # which is the case on Cori. Modify as needed for other machines. + # Also assuming we can use multiple cores on a compute node. + hostname = os.uname()[1] + if hostname.startswith('nid'): + args = (["srun", "-n", "12", "ESMF_RegridWeightGen", + "-s", mali_scripfile, + "-d", ismip7_scripfile, + "-w", mapping_file, + "-m", method_remap, + "-i", "-64bit_offset", + "--dst_regional", "--src_regional"]) + else: + args = (["ESMF_RegridWeightGen", + "-s", mali_scripfile, + "-d", ismip7_scripfile, + "-w", mapping_file, + "-m", method_remap, + "-i", "-64bit_offset", + "--dst_regional", "--src_regional"]) + + print(f"Running remapping command: {' '.join(args)}") + check_call(args) + + # remove the temporary scripfiles once the mapping file is generated + print("Removing the temporary mesh and scripfiles...") + os.remove(ismip7_scripfile) + os.remove(mali_scripfile) diff --git a/landice/output_processing_li/ismip7_postprocessing/post_process_mali_to_ismip7.py b/landice/output_processing_li/ismip7_postprocessing/post_process_mali_to_ismip7.py new file mode 100644 index 000000000..5fb0bd02e --- /dev/null +++ b/landice/output_processing_li/ismip7_postprocessing/post_process_mali_to_ismip7.py @@ -0,0 +1,292 @@ +#!/usr/bin/env python + +""" +This script processes MALI simulation outputs into the required format +for submission to ISMIP7. +There are 3 flavors of variables that can be processed: +1D, 2D state, and 2D flux variables. +Any combination can be specified. +Processing assumes simulations were run with the output streams defined +by the ismip7_run compass test case. +Multifile output file sets should be specified with a glob pattern. +Note: wildcard paths should be quoted to avoid shell expansion. + +An ISMIP submimssion grid resolution and grid file need to be specified. +A mapping file will be created unless an existing one is specified. + +More info at: +https://github.com/ismip/ISM_SimulationChecker/blob/main/conventions/ISMIP7_variable_request.csv +https://www.ismip.org/research/ismip7 + +Example usage from testing on initial ISMIP7 submission data: +python post_process_mali_to_ismip7.py \ + -o /pscratch/sd/h/hoffman2/ISMIP7-postprocessing-June30-deadline/AIS/test \ + -g "/pscratch/sd/t/trhille/ISMIP7/AIS_runs/no_slm_20260616/landice/ismip7_run/ismip7_ais/historical_CESM2-WACCM/output/globalStats_*.nc" \ + --icesheet AIS \ + -e C001 \ + --res 04 \ + -m /pscratch/sd/t/trhille/ISMIP7/AIS_runs/no_slm_20260616/landice/ismip7_run/ismip7_ais/historical_CESM2-WACCM/relaxed_10yrs_4km.nc \ + -s "/pscratch/sd/t/trhille/ISMIP7/AIS_runs/no_slm_20260616/landice/ismip7_run/ismip7_ais/historical_CESM2-WACCM/output/output_state_*.nc" \ + --ismip7_grid_file /pscratch/sd/h/hoffman2/ISMIP7-postprocessing-June30-deadline/AIS/misc/af2_AIS_04000m_v1.nc \ + -f "/pscratch/sd/t/trhille/ISMIP7/AIS_runs/no_slm_20260616/landice/ismip7_run/ismip7_ais/historical_CESM2-WACCM/output/output_state_*.nc" + +And here is an example of only processing 2d state and flux after the mapping file has been created: +python post_process_mali_to_ismip7.py \ + -o /pscratch/sd/h/hoffman2/ISMIP7-postprocessing-June30-deadline/AIS/test \ + --icesheet AIS \ + -e C001 \ + --res 04 \ + -m /pscratch/sd/t/trhille/ISMIP7/AIS_runs/no_slm_20260616/landice/ismip7_run/ismip7_ais/historical_CESM2-WACCM/relaxed_10yrs_4km.nc \ + -s "/pscratch/sd/t/trhille/ISMIP7/AIS_runs/no_slm_20260616/landice/ismip7_run/ismip7_ais/historical_CESM2-WACCM/output/output_state_*.nc" \ + -f "/pscratch/sd/t/trhille/ISMIP7/AIS_runs/no_slm_20260616/landice/ismip7_run/ismip7_ais/historical_CESM2-WACCM/output/output_state_*.nc" \ + --ismip7_grid_file /pscratch/sd/h/hoffman2/ISMIP7-postprocessing-June30-deadline/AIS/misc/af2_AIS_04000m_v1.nc \ + --reuse_mapping_file mapping_mali_to_ismip7.conserve.20260626T130744.nc + +""" + +import argparse +import glob +import os +from datetime import datetime + +from grid_and_mapping import ( + build_mapping_file, + check_exp_name, + check_ismip7_grid_file, + check_res, +) +from process_1d_variables_ismip7 import ( + check_global_stats_files, + generate_output_1d_vars, +) +from process_flux_variables_ismip7 import process_flux_pipeline +from process_state_variables_ismip7 import process_state_pipeline + +DEFAULT_AUTHORS = 'Matthew Hoffman, Trevor Hillebrand, Holly Kyeore Han' +DEFAULT_GROUP = 'Los Alamos National Laboratory, U.S. Department of Energy' +DEFAULT_MODEL = 'MALI (MPAS-Albany Land Ice model)' +DEFAULT_GROUP_NICKNAME = 'DOE' + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument( + "-e", + "--exp_name", + dest="exp", + required=True, + help="ISMIP7 experiment name (e.g., exp05)", + ) + parser.add_argument( + "-s", + "--input_state_pattern", + dest="input_state_pattern", + required=False, + help=( + "glob pattern matching one or more MALI state output files. " + "Note: wildcard paths should be quoted to avoid shell expansion." + ), + ) + parser.add_argument( + "-f", + "--input_flux_pattern", + dest="input_flux_pattern", + required=False, + help=( + "glob pattern matching one or more MALI flux output files. " + "Note: wildcard paths should be quoted to avoid shell expansion." + ), + ) + parser.add_argument( + "-m", + "--input_mesh", + dest="input_file_mesh", + required=False, + help="MALI file with mesh information", + ) + parser.add_argument( + "-g", + "--global_stats_pattern", + dest="global_stats_pattern", + required=False, + help=( + "glob pattern matching one or more globalStats.nc files " + "(e.g. 'globalStats_*.nc'). " + "Note: wildcard paths should be quoted to avoid shell expansion." + ), + ) + parser.add_argument( + "-o", + "--output_path", + dest="output_path", + required=True, + help="path to which the final output files will be saved", + ) + parser.add_argument( + "--reuse_mapping_file", + dest="reuse_mapping_file", + required=False, + help="existing mapping file name to reuse", + ) + parser.add_argument( + "--ismip7_grid_file", + dest="ismip7_grid_file", + required=True, + help="Input ismip7 mesh file.", + ) + parser.add_argument( + "--method", + dest="method_remap", + default="conserve", + required=False, + help="mapping method. Default='conserve'", + ) + parser.add_argument( + "--res", + dest="res_ismip7_grid", + required=True, + help="resolution of the ismip7 grid, in kilometers: 16, 8, 4, 2, 1", + ) + parser.add_argument( + "--icesheet", + dest="icesheet", + required=True, + choices=['AIS', 'GIS'], + help="ice sheet domain: 'AIS' (Antarctica) or 'GIS' (Greenland)", + ) + parser.add_argument( + "--authors", + dest="authors", + required=False, + default=DEFAULT_AUTHORS, + help=( + "author string for output file metadata " + f"(default: '{DEFAULT_AUTHORS}')" + ), + ) + parser.add_argument( + "--group", + dest="group", + required=False, + default=DEFAULT_GROUP, + help=( + "group/institution string for output file metadata " + f"(default: '{DEFAULT_GROUP}')" + ), + ) + parser.add_argument( + "--group_nickname", + dest="group_nickname", + required=False, + default=DEFAULT_GROUP_NICKNAME, + help=( + "short group nickname used in output filenames " + f"(default: '{DEFAULT_GROUP_NICKNAME}')" + ), + ) + args = parser.parse_args() + + check_exp_name(args.exp) + check_res(args.res_ismip7_grid) + + metadata = { + 'exp': args.exp, + 'icesheet': args.icesheet, + 'authors': args.authors, + 'group': args.group, + 'group_nickname': args.group_nickname, + 'model': DEFAULT_MODEL, + 'date': datetime.now().strftime("%d-%b-%Y"), + } + + print("\n---Processing remapping file---") + # Only do remapping steps if we have 2d files to process + if ( + args.input_state_pattern is not None or + args.input_flux_pattern is not None): + + check_ismip7_grid_file(args.ismip7_grid_file, args.res_ismip7_grid) + + method_remap = args.method_remap + if args.reuse_mapping_file is not None: + if not os.path.exists(args.reuse_mapping_file): + raise FileNotFoundError( + "Mapping file to reuse not found: " + f"{args.reuse_mapping_file}" + ) + print(f"Reusing existing mapping file: {args.reuse_mapping_file}") + mapping_file = args.reuse_mapping_file + else: + if args.input_file_mesh is None: + raise ValueError( + "--input_mesh is required when creating " + "a new mapping file." + ) + + created_at = datetime.now().strftime("%Y%m%dT%H%M%S") + mapping_file = ( + f"mapping_mali_to_ismip7.{method_remap}.{created_at}.nc" + ) + + print("Creating new mapping file. " + f"Mapping method used: {method_remap}") + build_mapping_file( + args.input_file_mesh, + mapping_file, + args.res_ismip7_grid, + args.ismip7_grid_file, + method_remap, + ) + + print("---Processing remapping file complete---\n") + + output_path = args.output_path + print(f"Using output path: {output_path}") + if not os.path.isdir(output_path): + os.makedirs(output_path) + + if args.global_stats_pattern is None: + print("--- No global stats pattern provided; " + "skipping 1D variable processing.") + else: + print("\n---Processing global stats file(s)---") + global_stats_files = sorted(glob.glob(args.global_stats_pattern)) + check_global_stats_files(global_stats_files) + generate_output_1d_vars(global_stats_files, output_path, metadata) + print("---Processing global stats file(s) complete---\n") + + if args.input_state_pattern is None: + print("--- MALI state pattern is not provided, " + "thus it will not be processed.") + else: + print("\n---Processing state file(s)---") + state_files = sorted(glob.glob(args.input_state_pattern)) + process_state_pipeline( + state_files, + mapping_file, + args.ismip7_grid_file, + output_path, + metadata, + ) + print("---Processing state file(s) complete---\n") + + if args.input_flux_pattern is None: + print("--- MALI flux pattern is not provided, " + "thus it will not be processed.") + else: + print("\n---Processing flux file(s)---") + flux_files = sorted(glob.glob(args.input_flux_pattern)) + process_flux_pipeline( + flux_files, + mapping_file, + args.ismip7_grid_file, + output_path, + metadata, + ) + print("---Processing flux file(s) complete---\n") + + print("---All processing complete---") + + +if __name__ == "__main__": + main() diff --git a/landice/output_processing_li/ismip7_postprocessing/process_1d_variables_ismip7.py b/landice/output_processing_li/ismip7_postprocessing/process_1d_variables_ismip7.py new file mode 100644 index 000000000..97934fb81 --- /dev/null +++ b/landice/output_processing_li/ismip7_postprocessing/process_1d_variables_ismip7.py @@ -0,0 +1,372 @@ +""" +Functions for processing and writing 1D (scalar time-series) output variables +for ISMIP7 submissions from MALI globalStats output. +""" + +from netCDF4 import Dataset +import numpy as np +import os +import sys +import xarray as xr +from validate import validate_mali_files + + +EXPECTED_VARIABLES = [ + 'daysSinceStart', 'deltat', 'simulationStartTime', + 'totalIceVolume', 'volumeAboveFloatation', + 'groundedIceArea', 'floatingIceArea', + 'totalSfcMassBal', 'totalGroundedBasalMassBal', + 'totalFloatingBasalMassBal', 'totalCalvingFlux', + 'totalFaceMeltingFlux', 'groundingLineFlux', +] + + +def check_global_stats_files(files): + """ + Validate a list of globalStats files before processing. + + See validate.validate_mali_files for full details of checks performed. + """ + validate_mali_files(files, EXPECTED_VARIABLES, label='globalStats') + + +def _build_output_filename(output_path, varname, metadata): + """Build a standard ISMIP7 output filename for a given variable.""" + return ( + f'{output_path}/{varname}_{metadata["icesheet"]}_' + f'{metadata["group_nickname"]}_MALI_{metadata["exp"]}.nc' + ) + + +def _write_state_var( + varname, + data_values, + time_days, + standard_name, + units, + variable_desc, + output_path, + simulationStartDate, + metadata): + """Write one 1D state (snapshot) variable to NETCDF4_CLASSIC.""" + nt = len(data_values) + filename = _build_output_filename(output_path, varname, metadata) + ds_out = Dataset(filename, 'w', format='NETCDF4_CLASSIC') + ds_out.createDimension('time', nt) + var_out = ds_out.createVariable(varname, 'd', ('time',)) + time_out = ds_out.createVariable('time', 'd', ('time',)) + var_out[:] = data_values + time_out[:] = time_days + time_out.units = f'days since {simulationStartDate}' + time_out.calendar = 'noleap' + time_out.standard_name = 'time' + time_out.long_name = 'time' + var_out.standard_name = standard_name + var_out.units = units + ds_out.AUTHORS = metadata['authors'] + ds_out.MODEL = metadata['model'] + ds_out.GROUP = metadata['group'] + ds_out.VARIABLE = variable_desc + ds_out.DATE = metadata['date'] + ds_out.close() + + +def _write_flux_var(varname, data_values, days_min, days_max, standard_name, + units, variable_desc, output_path, simulationStartDate, + metadata): + """Write one 1D flux (time-averaged) variable to NETCDF4_CLASSIC.""" + nt = len(data_values) + filename = _build_output_filename(output_path, varname, metadata) + ds_out = Dataset(filename, 'w', format='NETCDF4_CLASSIC') + ds_out.createDimension('time', nt) + ds_out.createDimension('bnds', 2) + var_out = ds_out.createVariable(varname, 'd', ('time',)) + time_out = ds_out.createVariable('time', 'd', ('time',)) + bnds_out = ds_out.createVariable('time_bnds', 'd', ('time', 'bnds')) + var_out[:] = data_values + time_out[:] = (days_min + days_max) / 2.0 + bnds_out[:, 0] = days_min + bnds_out[:, 1] = days_max + time_out.units = f'days since {simulationStartDate}' + time_out.calendar = 'noleap' + time_out.standard_name = 'time' + time_out.long_name = 'time' + var_out.standard_name = standard_name + var_out.units = units + ds_out.AUTHORS = metadata['authors'] + ds_out.MODEL = metadata['model'] + ds_out.GROUP = metadata['group'] + ds_out.VARIABLE = variable_desc + ds_out.DATE = metadata['date'] + ds_out.close() + + +def generate_output_1d_vars(files, output_path, metadata): + """ + Process and write 1D (scalar time-series) state and flux variables. + + Parameters + ---------- + files : list of str + Sorted list of globalStats.nc file paths to process. + output_path : str + Directory for output files. + metadata : dict + Submission metadata with keys: + exp, icesheet, authors, group, model, date. + """ + + if output_path is None or not os.path.exists(output_path): + output_path = os.getcwd() + + ds = xr.open_mfdataset(files, combine='nested', concat_dim='Time', + decode_cf=False, data_vars='minimal', + coords='minimal', compat='override') + with xr.open_dataset(files[0], decode_cf=False) as ds_first: + simulationStartTime = ( + ds_first['simulationStartTime'] + .values.tobytes().decode('utf-8').strip().strip('\x00') + ) + daysSinceStart = ds['daysSinceStart'].values + dt = ds['deltat'].values + simulationStartDate = simulationStartTime.split("_")[0] + if simulationStartDate[5:10] != '01-01': + ds.close() + sys.exit( + "Error: simulationStartTime for globalStats file " + "is not on Jan. 1." + ) + refYear = int(simulationStartDate[0:4]) + decYears = refYear + daysSinceStart / 365.0 + endYr = decYears[-1] + if endYr != np.round(endYr): + ds.close() + sys.exit("Error: end year not an even year in globalStats file.") + + # Determine processed time levels for state and flux fields + # The historical state fields should include the initial time (Jan. 1). + # Projection state fields should not include the initial time (Jan. 1) + # of the projection because it's a restart from the historical. + # Flux fields should never use the Jan. 1 time level at the start of the + # year as part of the averaging. + # For year conventions here, for state fields, the year is the snapshot at + # the start of the year, e.g., state year 2000 means the snapshot at + # Jan. 1, 2000. + # For flux fields, the years is the calendar year being averaged over, + # e.g., flux year 2000 is the average between Jan. 1, 2000, + # and Jan. 1, 2001. + # Note this year convention differs from the first column in table in + # A2.3.2 at + # https://www.climate-cryosphere.org/wiki/index.php?title=ISMIP7-Projections2300-Antarctica#A2.3.3_Table_A1:_Variable_request_for_ISMIP6 + # but that year indexing convention ultimately doesn't matter because the + # time coordinates in these files uses units of days since a + # reference date, + # and it does not use a year indexing convention at all. + if decYears[0] == np.round(decYears[0]): + # The initial time level will only be on an even year (Jan. 1) + # for the hist run. In that case, we want to include that initial + # even year in the state processing. We also want the state snapshot + # at the final (even) year in the output. + # The flux processing should start with the first year, which covers a + # full 12 months. We exclude the final year, which is just a Jan. 1 + # posting. + years_state = np.arange(decYears[0], endYr + 1) + years_flux = np.arange(decYears[0], endYr) + else: + # For projection runs, the first state snapshot we want is the + # first Jan. 1, + # which we be the first even year after the initial time in the file. + # For flux files, the first full year we want to process is the + # year of the + # first time level in the file. As with hist, we exclude the + # final year, + # which is just a Jan. 1 posting. + years_state = np.arange(np.ceil(decYears[0]), endYr + 1) + years_flux = np.arange(np.floor(decYears[0]), endYr) + nt_state = len(years_state) + nt_flux = len(years_flux) + print( + "For state processing, using start " + f"year={years_state[0]} and end year={years_state[-1]}." + ) + print( + "For flux processing, using start " + f"year={years_flux[0]} and end year={years_flux[-1]}." + ) + + # read in state variables + vol = ds['totalIceVolume'].values + vaf = ds['volumeAboveFloatation'].values + gia = ds['groundedIceArea'].values + fia = ds['floatingIceArea'].values + + # read in flux variables over which yearly average will be taken + smb = ds['totalSfcMassBal'].values + bmbGr = ds['totalGroundedBasalMassBal'].values.copy() + # clean out some garbage values we can't account for + ind = np.nonzero(bmbGr > 1.0e18)[0] + if len(ind) > 0: + print( + f"WARNING: Found { + len(ind)} values of totalGroundedBasalMassBal>1.0e18") + bmbGr[ind] = np.nan + ind = np.nonzero(bmbGr < -1.0e18)[0] + if len(ind) > 0: + print( + f"WARNING: Found { + len(ind)} values of totalGroundedBasalMassBal<-1.0e18") + bmbGr[ind] = np.nan + bmbFlt = ds['totalFloatingBasalMassBal'].values.copy() + # clean out some garbage values we can't account for + ind = np.nonzero(bmbFlt > 1.0e18)[0] + if len(ind) > 0: + print( + f"WARNING: Found { + len(ind)} values of totalFloatingBasalMassBal>1.0e18") + bmbFlt[ind] = np.nan + ind = np.nonzero(bmbFlt < -1.0e18)[0] + if len(ind) > 0: + print( + f"WARNING: Found { + len(ind)} values of totalFloatingBasalMassBal<-1.0e18") + bmbFlt[ind] = np.nan + cfx = ds['totalCalvingFlux'].values + fmfx = ds['totalFaceMeltingFlux'].values + gfx = ds['groundingLineFlux'].values + + # initialize 1D variables that will store data value on the + # January 1st of each year + vol_snapshot = np.zeros(nt_state) * np.nan + vaf_snapshot = np.zeros(nt_state) * np.nan + gia_snapshot = np.zeros(nt_state) * np.nan + fia_snapshot = np.zeros(nt_state) * np.nan + days_snapshot = np.zeros(nt_state) * np.nan + smb_avg = np.zeros(nt_flux) * np.nan + bmbGr_avg = np.zeros(nt_flux) * np.nan + bmbFlt_avg = np.zeros(nt_flux) * np.nan + cfx_avg = np.zeros(nt_flux) * np.nan + fmfx_avg = np.zeros(nt_flux) * np.nan + gfx_avg = np.zeros(nt_flux) * np.nan + days_min = np.zeros(nt_flux) * np.nan + days_max = np.zeros(nt_flux) * np.nan + + # this is for the state variables + for i in range(nt_state): + # Use isclose to avoid floating-point equality issues. + ind_snap = np.where(np.isclose(decYears, years_state[i]))[0] + if len(ind_snap) == 0: + raise ValueError( + f"No state snapshot found for year { + years_state[i]}.") + if len(ind_snap) > 1: + print(f"WARNING: Found {len(ind_snap)} snapshots for year " + f"{years_state[i]}; using the first one.") + idx_snap = ind_snap[0] + + vol_snapshot[i] = vol[idx_snap] + vaf_snapshot[i] = vaf[idx_snap] + gia_snapshot[i] = gia[idx_snap] + fia_snapshot[i] = fia[idx_snap] + days_snapshot[i] = daysSinceStart[idx_snap] + + if decYears[idx_snap] == endYr: + break + + # this is for the flux variables + for i in range(nt_flux): + ind_avg = np.where( + np.logical_and( + decYears > years_flux[i], + decYears <= ( + years_flux[i] + + 1.0)))[0] + if len(ind_avg) == 0: + raise ValueError(f"No flux averaging samples found for year " + f"{years_flux[i]}.") + smbi = smb[ind_avg] + bmbGri = bmbGr[ind_avg] + bmbFlti = bmbFlt[ind_avg] + cfxi = cfx[ind_avg] + fmfxi = fmfx[ind_avg] + gfxi = gfx[ind_avg] + dti = dt[ind_avg] + + # take the average of the flux variables + smb_avg[i] = np.nansum(smbi * dti) / np.nansum(dti) + bmbGr_avg[i] = np.nansum(bmbGri * dti) / np.nansum(dti) + bmbFlt_avg[i] = np.nansum(bmbFlti * dti) / np.nansum(dti) + cfx_avg[i] = np.nansum(cfxi * dti) / np.nansum(dti) + fmfx_avg[i] = np.nansum(fmfxi * dti) / np.nansum(dti) + gfx_avg[i] = np.nansum(gfxi * dti) / np.nansum(dti) + days_min[i] = (years_flux[i] - refYear) * 365.0 + days_max[i] = (years_flux[i] + 1.0 - refYear) * 365.0 + + if decYears[ind_avg][-1] == endYr: + break + + # shared keyword arguments for both writer helpers + common = dict( + output_path=output_path, + simulationStartDate=simulationStartDate, + metadata=metadata, + ) + + # --- state (snapshot) variables --- + _write_state_var('lim', vol_snapshot * 910, days_snapshot, + 'land_ice_mass', 'kg', 'Total ice mass', **common) + _write_state_var('limnsw', vaf_snapshot * 910, days_snapshot, + 'land_ice_mass_not_displacing_sea_water', 'kg', + 'Mass above floatation', **common) + _write_state_var( + 'iareagr', + gia_snapshot, + days_snapshot, + 'grounded_ice_sheet_area', + 'm2', + 'Grounded ice area', + **common) + _write_state_var( + 'iareafl', + fia_snapshot, + days_snapshot, + 'floating_ice_shelf_area', + 'm2', + 'Floating ice area', + **common) + + # --- flux (time-averaged) variables --- + _write_flux_var('tendacabf', smb_avg / 31536000.0, days_min, days_max, + 'tendency_of_land_ice_mass_due_to_surface_mass_balance', + 'kg s-1', 'Total SMB flux', **common) + _write_flux_var( + 'tendlibmassbfgr', + bmbGr_avg / 31536000.0, + days_min, + days_max, + 'tendency_of_land_ice_mass_due_to_basal_mass_balance', + 'kg s-1', + 'Total BMB flux beneath grounded ice', + **common) + _write_flux_var( + 'tendlibmassbffl', + bmbFlt_avg / 31536000.0, + days_min, + days_max, + 'tendency_of_land_ice_mass_due_to_basal_mass_balance', + 'kg s-1', + 'Total BMB flux beneath floating ice', + **common) + # tendlicalvf: sign convention — calving removes mass, so negate + _write_flux_var('tendlicalvf', -cfx_avg / 31536000.0, days_min, days_max, + 'tendency_of_land_ice_mass_due_to_calving', + 'kg s-1', 'Total calving flux', **common) + # tendlifmassbf: in ISMIP7 this is ice-front melting only (not calving) + _write_flux_var('tendlifmassbf', -fmfx_avg / 31536000.0, days_min, + days_max, + 'tendency_of_land_ice_mass_due_to_ice_front_melting', + 'kg s-1', 'Total ice front melting flux', **common) + _write_flux_var('tendligroundf', gfx_avg / 31536000.0, days_min, days_max, + 'tendency_of_grounded_ice_mass', + 'kg s-1', 'Total grounding line flux', **common) + + ds.close() diff --git a/landice/output_processing_li/ismip7_postprocessing/process_flux_variables_ismip7.py b/landice/output_processing_li/ismip7_postprocessing/process_flux_variables_ismip7.py new file mode 100644 index 000000000..0b9b75c80 --- /dev/null +++ b/landice/output_processing_li/ismip7_postprocessing/process_flux_variables_ismip7.py @@ -0,0 +1,307 @@ +""" +This script has functions that are needed to post-process and write flux +output variables from ISMIP7 simulations. +""" + +from netCDF4 import Dataset +import xarray as xr +import numpy as np +from subprocess import check_call +import os +from validate import validate_mali_files + + +EXPECTED_FLUX_VARIABLES = [ + 'daysSinceStart', 'simulationStartTime', 'cellMask', + 'avgSMBFlux', 'avgFloatingBMBFlux', 'avgGroundedBMBFlux', + 'avgDhdt', 'avgCalvingFlux', 'avgFaceMeltFlux', + 'avgGroundingLineFlux', +] + +REQUIRED_FLUX_REMAPPING_VARIABLES = [ + 'daysSinceStart', 'simulationStartTime', + 'avgSMBFlux', 'avgFloatingBMBFlux', 'avgGroundedBMBFlux', + 'avgDhdt', 'avgCalvingFlux', 'avgFaceMeltFlux', + 'avgGroundingLineFlux', + # masks not currently used - see note below in process_flux_vars + #'ice_mask', 'floating_mask', 'dynamic_mask', 'grounded_mask', + #'gl_mask' +] + + +def check_flux_files(files): + """ + Validate a list of MALI flux output files before processing. + + See validate.validate_mali_files for full details of checks performed. + """ + validate_mali_files(files, EXPECTED_FLUX_VARIABLES, label='flux') + + +def process_flux_vars(files, tmp_file): + """ + Prepare flux files into a temporary file for remapping. + This is very simple, because flux variables are already + time-averaged by MALI. + + Parameters + ---------- + files : list of str + Sorted list of MALI flux output file paths. + tmp_file : str + Temporary output file name. + """ + ds_flux = xr.open_mfdataset(files, combine='nested', + concat_dim='Time', + engine='netcdf4', + decode_cf=False, + data_vars='minimal', + coords='minimal', + compat='override') + if 'units' in ds_flux['daysSinceStart'].attrs: + del ds_flux['daysSinceStart'].attrs['units'] + + # variables that need to be modified prior to remapping + # these masks are not currently used, but would be needed if the + # flux vars are meant to be masked. + # However there is a potential issue that fluxes could have occurred + # in locations that fall outside of the final mask, because + # the fluxes are time-averaged over a year, and the mask is only + # valid at the end of the year. + #ds_flux['ice_mask'] = (ds_flux['cellMask'][:, :] & 2) / 2 # grounded: dynamic ice + #ds_flux['floating_mask'] = (ds_flux['cellMask'][:, :] & 4) / 4 + #ds_flux['dynamic_mask'] = (ds_flux['cellMask'][:, :] & 2) / 2 + #ds_flux['grounded_mask'] = (ds_flux['cellMask'][:, :] * 0 + 1) - ds_flux['floating_mask'] + #ds_flux['gl_mask'] = (ds_flux['cellMask'][:, :] & 256) / 256 + + missing = [ + var for var in REQUIRED_FLUX_REMAPPING_VARIABLES + if var not in ds_flux + ] + if missing: + raise ValueError( + "Processed flux dataset is missing required output variables: " + f"{missing}" + ) + + # subset to only required variables to keep file small for remapping + ds_flux_out = ds_flux[REQUIRED_FLUX_REMAPPING_VARIABLES] + # remove the first time step (which is always 0) + time_mask = ~np.isclose(ds_flux_out['daysSinceStart'].values, 0.0) + if not np.any(time_mask): + raise ValueError("No flux time records remain after dropping daysSinceStart==0.") + ds_flux_out = ds_flux_out.isel(Time=time_mask) + + ds_flux_out.to_netcdf(tmp_file) + + ds_flux_out.close() + ds_flux.close() + + +def write_netcdf_2d_flux_vars(mali_var_name, ismip7_var_name, var_std_name, + var_units, var_varname, remapped_mali_flux_file, + ismip7_grid_file, output_path, metadata): + """ + mali_var_name: variable name on MALI side + ismip7_var_name: variable name required by ISMIP7 + var_std_name: standard variable name + var_units: variable units + var_varname: variable variable name + remapped_mali_flux_file: mali flux file remapped on the ISMIP7 grid + ismip7_grid_file: original ISMIP7 file + exp: experiment name + output_path: output path to which the output files will be saved + """ + + data_ismip7 = Dataset(ismip7_grid_file, 'r') + var_x = data_ismip7.variables['x'][:] + var_y = data_ismip7.variables['y'][:] + + data = Dataset(remapped_mali_flux_file, 'r') + data.set_auto_mask(False) + simulationStartTime = data.variables['simulationStartTime'][:].tobytes( + ).decode('utf-8').strip().strip('\x00') + simulationStartDate = simulationStartTime.split("_")[0] + daysSinceStart = data.variables['daysSinceStart'][:] + refYear = int(simulationStartDate[0:4]) + decYears = refYear + daysSinceStart / 365.0 + years_flux = np.floor(decYears) + + # Flux outputs are annual means over each calendar year. + # Bounds are Jan 1 of the previous year through Jan 1 of the year indexed. + timeBndsMin = (years_flux - refYear - 1.0) * 365.0 + timeBndsMax = (years_flux - refYear) * 365.0 + if mali_var_name not in data.variables: + print(f"WARNING: {mali_var_name} not present. Skipping.") + data.close() + return + var_mali = data.variables[mali_var_name][:, :, :] + var_mali[np.where(abs(var_mali + 1e34) < 1e33)] = np.nan + timeSteps, latN, lonN = np.shape(var_mali) + + dataOut = Dataset( + f'{output_path}/{ismip7_var_name}_{ + metadata["icesheet"]}_{ + metadata["group_nickname"]}_MALI_{ + metadata["exp"]}.nc', + 'w', + format='NETCDF4_CLASSIC') + dataOut.createDimension('time', timeSteps) + dataOut.createDimension('bnds', 2) + timebndsValues = dataOut.createVariable('time_bnds', 'd', ('time', 'bnds')) + dataOut.createDimension('x', lonN) + dataOut.createDimension('y', latN) + dataValues = dataOut.createVariable(ismip7_var_name, 'd', + ('time', 'y', 'x'), fill_value=np.nan) + xValues = dataOut.createVariable('x', 'd', ('x')) + yValues = dataOut.createVariable('y', 'd', ('y')) + timeValues = dataOut.createVariable('time', 'd', ('time')) + + for i in range(timeSteps): + # mask not currently used - see note above in process_flux_vars + #mask = data.variables['ice_mask'][i, :, :] + tmp = var_mali[i, :, :] + #tmp[mask == 0] = np.nan + dataValues[i, :, :] = tmp + timeValues[i] = (timeBndsMin[i] + timeBndsMax[i]) / 2.0 + timebndsValues[i, 0] = timeBndsMin[i] + timebndsValues[i, 1] = timeBndsMax[i] + + for i in range(latN): + xValues[i] = var_x[i] + + for i in range(lonN): + yValues[i] = var_y[i] + + dataValues.standard_name = var_std_name + dataValues.units = var_units + timeValues.bounds = 'time_bnds' + timeValues.units = f'days since {simulationStartDate}' + timeValues.calendar = 'noleap' + timeValues.standard_name = 'time' + timeValues.long_name = 'time' + timebndsValues.units = f'days since {simulationStartDate}' + timebndsValues.calendar = 'noleap' + xValues.units = 'm' + xValues.standard_name = 'x' + xValues.long_name = 'x' + yValues.units = 'm' + yValues.standard_name = 'y' + yValues.long_name = 'y' + dataOut.AUTHORS = metadata['authors'] + dataOut.MODEL = metadata['model'] + dataOut.GROUP = metadata['group'] + dataOut.VARIABLE = var_varname + dataOut.DATE = metadata['date'] + dataOut.close() + data.close() + + +def generate_output_2d_flux_vars(file_remapped_mali_flux, + ismip7_grid_file, output_path, metadata): + """ + file_remapped_mali_flux: flux output file on mali mesh remapped + onto the ismip7 grid + ismip7 grid + ismip7_grid_file: ismip7 original file + exp: ISMIP7 experiment name + output_path: path to which the final output files are saved + """ + + # ----------- acabf ------------------ + write_netcdf_2d_flux_vars('avgSMBFlux', 'acabf', + 'land_ice_surface_specific_mass_balance_flux', + 'kg m-2 s-1', 'Surface mass balance flux', + file_remapped_mali_flux, + ismip7_grid_file, output_path, metadata) + + # ----------- libmassbffl ------------------ + write_netcdf_2d_flux_vars('avgFloatingBMBFlux', 'libmassbffl', + 'land_ice_basal_specific_mass_balance_flux', + 'kg m-2 s-1', + 'Basal mass balance flux beneath floating ice', + file_remapped_mali_flux, + ismip7_grid_file, output_path, metadata) + + # ----------- libmassbfgr ------------------ + write_netcdf_2d_flux_vars('avgGroundedBMBFlux', 'libmassbfgr', + 'land_ice_basal_specific_mass_balance_flux', + 'kg m-2 s-1', + 'Basal mass balance flux beneath grounded ice', + file_remapped_mali_flux, + ismip7_grid_file, output_path, metadata) + + # ----------- dlithkdt ------------------ + write_netcdf_2d_flux_vars('avgDhdt', 'dlithkdt', + 'tendency_of_land_ice_thickness', + 'm s-1', + 'Ice thickness imbalance', + file_remapped_mali_flux, + ismip7_grid_file, output_path, metadata) + + # ----------- licalvf ------------------ + write_netcdf_2d_flux_vars('avgCalvingFlux', 'licalvf', + 'land_ice_specific_mass_flux_due_to_calving', + 'kg m-2 s-1', + 'Calving flux', + file_remapped_mali_flux, + ismip7_grid_file, output_path, metadata) + + # ----------- lifmassbf ------------------ + # Note: facemelting and calving flux are combined above + write_netcdf_2d_flux_vars('avgFaceMeltFlux', 'lifmassbf', + 'TBD by ISMIP7', + 'kg m-2 s-1', + 'Ice front melt flux', + file_remapped_mali_flux, + ismip7_grid_file, output_path, metadata) + + # ----------- ligroundf ------------------ + write_netcdf_2d_flux_vars('avgGroundingLineFlux', 'ligroundf', + 'TBD by ISMIP7', + 'kg m-2 s-1', + 'Grounding line flux', + file_remapped_mali_flux, + ismip7_grid_file, output_path, metadata) + + +def process_flux_pipeline(flux_files, mapping_file, ismip7_grid_file, + output_path, metadata): + """ + Full flux-variable processing pipeline: + validate, concatenate, remap, write. + + Parameters + ---------- + flux_files : list of str + Sorted list of MALI flux output file paths. + mapping_file : str + Path to the ESMF mapping/weights file. + ismip7_grid_file : str + Path to the ISMIP7 grid file. + output_path : str + Directory for output files. + metadata : dict + Submission metadata dict. + """ + check_flux_files(flux_files) + + print("Preparing concatenated flux file for remapping.") + tmp_flux_file = 'tmp_flux.nc' + process_flux_vars(flux_files, tmp_flux_file) + + print("Remapping flux file.") + remapped_file_flux = 'remapped_flux.nc' + check_call(["ncremap", + "-i", tmp_flux_file, + "-o", remapped_file_flux, + "-m", mapping_file, + "-P", "mpas"]) + + print("Writing processed and remapped flux fields to ISMIP7 file format.") + generate_output_2d_flux_vars(remapped_file_flux, + ismip7_grid_file, + output_path, metadata) + + os.remove(tmp_flux_file) + os.remove(remapped_file_flux) diff --git a/landice/output_processing_li/ismip7_postprocessing/process_state_variables_ismip7.py b/landice/output_processing_li/ismip7_postprocessing/process_state_variables_ismip7.py new file mode 100644 index 000000000..ef6166e94 --- /dev/null +++ b/landice/output_processing_li/ismip7_postprocessing/process_state_variables_ismip7.py @@ -0,0 +1,398 @@ +""" +This script has functions that are needed to post-process and write state +output variables from ISMIP7 simulations. +""" + +from netCDF4 import Dataset +import xarray as xr +import numpy as np +import os +from subprocess import check_call +from validate import validate_mali_files + + +EXPECTED_STATE_VARIABLES = [ + 'daysSinceStart', 'simulationStartTime', + 'cellMask', 'basalTemperature', 'betaSolve', + 'uReconstructX', 'uReconstructY', 'upperSurface', +] + +# Keep only variables needed downstream by generate_output_2d_state_vars +# and write_netcdf_2d_state_vars. +REQUIRED_STATE_REMAPPING_VARIABLES = [ + 'daysSinceStart', 'simulationStartTime', + 'thickness', 'upperSurface', 'lowerSurface', 'bedTopography', + 'uReconstructX_sfc', 'uReconstructY_sfc', + 'uReconstructX_base', 'uReconstructY_base', + 'xvelmean', 'yvelmean', 'surfaceTemperature', + 'litempbotgr', 'litempbotfl', 'strbasemag', + 'sftgif', 'sftgrf', 'sftflf', +] + + +def check_state_files(files): + """ + Validate a list of MALI state output files before processing. + + See validate.validate_mali_files for full details of checks performed. + """ + validate_mali_files(files, EXPECTED_STATE_VARIABLES, label='state') + + +def process_state_vars(files, tmp_file): + """ + files: list of MALI state output file paths + tmp_file: temporary output file name + """ + + inputfile_state_vars = xr.open_mfdataset(files, combine='nested', + concat_dim='Time', + engine='netcdf4', + decode_cf=False, + data_vars='minimal', + coords='minimal', + compat='override') + # Delete the 'units' attr on daysSinceStart to prevent xarray from + # encoding it as a timedelta type (corrupts values after ~250 years). + if 'units' in inputfile_state_vars['daysSinceStart'].attrs: + del inputfile_state_vars['daysSinceStart'].attrs['units'] + + # get the mesh description data + nLayer = inputfile_state_vars.sizes['nVertLevels'] + nInterface = nLayer + 1 # inputfile_state_vars.sizes['nVertInterfaces'] + cellMask = inputfile_state_vars['cellMask'][:, :] + basalTemperature = inputfile_state_vars['basalTemperature'][:, :] + betaSolve = inputfile_state_vars['betaSolve'][:, :] + + inputfile_state_vars['litempbotfl'] = basalTemperature * \ + (cellMask[:, :] & 4) / 4 + inputfile_state_vars['litempbotgr'] = basalTemperature * \ + (1 - (cellMask[:, :] & 4) / 4) + + uxsurf = inputfile_state_vars['uReconstructX'][:, :, 0] + uysurf = inputfile_state_vars['uReconstructY'][:, :, 0] + uxbase = inputfile_state_vars['uReconstructX'][:, :, nInterface - 1] + uybase = inputfile_state_vars['uReconstructY'][:, :, nInterface - 1] + inputfile_state_vars['uReconstructX_sfc'] = uxsurf + inputfile_state_vars['uReconstructY_sfc'] = uysurf + inputfile_state_vars['uReconstructX_base'] = uxbase + inputfile_state_vars['uReconstructY_base'] = uybase + + inputfile_state_vars['upperSurface'] = np.maximum( + 0.0, inputfile_state_vars['upperSurface']) + + floating_mask = (cellMask[:, :] & 4) / 4 + dynamic_mask = (cellMask[:, :] & 2) / 2 + grounded_mask = (cellMask[:, :] * 0 + 1) - floating_mask + + # floating and dynamic + inputfile_state_vars['sftflf'] = floating_mask * dynamic_mask + # grounded: not-floating and dynamic + inputfile_state_vars['sftgrf'] = grounded_mask * dynamic_mask + # dynamic ice + inputfile_state_vars['sftgif'] = dynamic_mask + + speed_base = (uxbase[:, :] ** 2 + uybase[:, :] ** 2) ** 0.5 + seconds_per_year = 3600.0 * 24.0 * 365.0 + inputfile_state_vars['strbasemag'] = ( + betaSolve[:, :] * speed_base * seconds_per_year * grounded_mask * + dynamic_mask + ) + + missing = [ + var for var in REQUIRED_STATE_REMAPPING_VARIABLES + if var not in inputfile_state_vars + ] + if missing: + raise ValueError( + "Processed state dataset is missing required output variables: " + f"{missing}" + ) + + output_state_vars = inputfile_state_vars[REQUIRED_STATE_OUTPUT_VARIABLES] + output_state_vars.to_netcdf(tmp_file) + output_state_vars.close() + inputfile_state_vars.close() + + +def write_netcdf_2d_state_vars( + mali_var_name, + ismip7_var_name, + var_std_name, + var_units, + var_varname, + remapped_mali_outputfile, + ismip7_grid_file, + output_path, + metadata): + """ + mali_var_name: variable name on MALI side + ismip7_var_name: variable name required by ISMIP7 + var_std_name: standard variable name + var_units: variable units + var_varname: variable variable name + remapped_mali_outputfile: mali state file remapped on the ISMIP7 grid + ismip7_grid_file: original ISMIP7 file + exp: experiment name + output_path: output path to which the output files will be saved + """ + + data_ismip7 = Dataset(ismip7_grid_file, 'r') + var_x = data_ismip7.variables['x'][:] + var_y = data_ismip7.variables['y'][:] + + data = Dataset(remapped_mali_outputfile, 'r') + data.set_auto_mask(False) + simulationStartTime = data.variables['simulationStartTime'][:].tobytes( + ).decode('utf-8').strip().strip('\x00') + simulationStartDate = simulationStartTime.split("_")[0] + daysSinceStart = data.variables['daysSinceStart'][:] + var_sftgif = data.variables['sftgif'][:, :, :] + var_sftgrf = data.variables['sftgrf'][:, :, :] + var_sftflf = data.variables['sftflf'][:, :, :] + var_mali = data.variables[mali_var_name][:, :, :] + var_mali[np.where(abs(var_mali + 1e34) < 1e33)] = np.nan + timeSteps, latN, lonN = np.shape(var_mali) + + dataOut = Dataset( + f'{output_path}/{ismip7_var_name}_{ + metadata["icesheet"]}_{ + metadata["group_nickname"]}_MALI_{ + metadata["exp"]}.nc', + 'w', + format='NETCDF4_CLASSIC') + dataOut.createDimension('time', timeSteps) + dataOut.createDimension('x', lonN) + dataOut.createDimension('y', latN) + dataValues = dataOut.createVariable(ismip7_var_name, 'd', + ('time', 'y', 'x'), fill_value=np.nan) + xValues = dataOut.createVariable('x', 'd', ('x')) + yValues = dataOut.createVariable('y', 'd', ('y')) + timeValues = dataOut.createVariable('time', 'd', ('time')) + timeValues[:] = daysSinceStart + for i in range(timeSteps): + if ismip7_var_name == 'sftgif': + dataValues[i, :, :] = var_mali[i, :, :] + else: + if ismip7_var_name == 'litempbotgr': + mask = var_sftgrf[i, :, :] + elif ismip7_var_name == 'litempbotfl': + mask = var_sftflf[i, :, :] + elif ismip7_var_name == 'topg': + mask = np.ones(var_mali.shape[1:]) # don't mask topg + else: + mask = var_sftgif[i, :, :] + tmp = var_mali[i, :, :] + tmp[mask == 0] = np.nan + dataValues[i, :, :] = tmp + + for i in range(latN): + xValues[i] = var_x[i] + + for i in range(lonN): + yValues[i] = var_y[i] + + dataValues.standard_name = var_std_name + dataValues.units = var_units + timeValues.units = f'days since {simulationStartDate}' + timeValues.calendar = 'noleap' + timeValues.standard_name = 'time' + timeValues.long_name = 'time' + xValues.units = 'm' + xValues.standard_name = 'x' + xValues.long_name = 'x' + yValues.units = 'm' + yValues.standard_name = 'y' + yValues.long_name = 'y' + dataOut.AUTHORS = metadata['authors'] + dataOut.MODEL = metadata['model'] + dataOut.GROUP = metadata['group'] + dataOut.VARIABLE = var_varname + dataOut.DATE = metadata['date'] + dataOut.close() + + +def process_state_pipeline(state_files, mapping_file, ismip7_grid_file, + output_path, metadata): + """ + Full state-variable processing pipeline: validate, adjust, remap, write. + + Parameters + ---------- + state_files : list of str + Sorted list of MALI state output file paths. + mapping_file : str + Path to the ESMF mapping/weights file. + ismip7_grid_file : str + Path to the ISMIP7 grid file. + output_path : str + Directory for output files. + metadata : dict + Submission metadata (exp, icesheet, authors, group, model, date, ...). + """ + check_state_files(state_files) + + print("Calculating needed state file adjustments.") + tmp_file = "tmp_state.nc" + process_state_vars(state_files, tmp_file) + + processed_and_remapped_file_state = 'processed_and_remapped_state.nc' + print("Remapping state file.") + check_call(["ncremap", + "-i", tmp_file, + "-o", processed_and_remapped_file_state, + "-m", mapping_file, + "-P", "mpas"]) + + print("Writing processed and remapped state fields to ISMIP7 file format.") + generate_output_2d_state_vars(processed_and_remapped_file_state, + ismip7_grid_file, output_path, metadata) + + os.remove(tmp_file) + os.remove(processed_and_remapped_file_state) + + +def generate_output_2d_state_vars(file_remapped_mali_state, + ismip7_grid_file, output_path, metadata): + """ + file_remapped_mali_state: output files on mali mesh remapped + on the ismip7 grid + ismip7_grid_file: ismip7 original file + exp: ISMIP7 experiment name + output_path: path to which the final output files are saved + """ + + # ----------- lithk ------------------ + write_netcdf_2d_state_vars('thickness', 'lithk', 'land_ice_thickness', + 'm', 'Ice thickness', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- orog ------------------ + write_netcdf_2d_state_vars('upperSurface', 'orog', 'surface_altitude', 'm', + 'Surface elevation', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- base ------------------ + write_netcdf_2d_state_vars('lowerSurface', 'base', 'base_altitude', 'm', + 'Base elevation', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- topg ------------------ + write_netcdf_2d_state_vars( + 'bedTopography', + 'topg', + 'bedrock_altitude', + 'm', + 'Bedrock elevation', + file_remapped_mali_state, + ismip7_grid_file, + output_path, + metadata) + + # ----------- hfgeoubed------------------ + # Note: even though this is a flux variable, we are taking a snapshot of it + # as it does not change with time. + # Uncomment once basalHeatFlux is outputted in the output stream. + # write_netcdf_2d_state_vars('basalHeatFlux', 'hfgeoubed', + # 'upward_geothermal_heat_flux_in_land_ice', + # 'W m-2', 'Geothermal heat flux', + # file_remapped_mali_state, ismip7_grid_file, + # exp, output_path) + + # ----------- xvelsurf ------------------ + write_netcdf_2d_state_vars('uReconstructX_sfc', 'xvelsurf', + 'land_ice_surface_x_velocity', + 'm s-1', 'Surface velocity in x', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # -----------yxvelsurf ------------------ + write_netcdf_2d_state_vars('uReconstructY_sfc', 'yvelsurf', + 'land_ice_surface_y_velocity', + 'm s-1', 'Surface velocity in x', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- xvelbase ------------------ + write_netcdf_2d_state_vars('uReconstructX_base', 'xvelbase', + 'land_ice_basal_x_velocity', + 'm s-1', 'Basal velocity in x', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- yvelbase ------------------ + write_netcdf_2d_state_vars('uReconstructY_base', 'yvelbase', + 'land_ice_basal_y_velocity', + 'm s-1', 'Basal velocity in y', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- zvelsurf & zvelbase ------------------ + # ISMIP7 requires these variables, but MALI does not output them. + # So, we are not processing/writing these variables out. + + # ----------- xvelmean ------------------ + write_netcdf_2d_state_vars('xvelmean', 'xvelmean', + 'land_ice_vertical_mean_x_velocity', + 'm s-1', 'Mean velocity in x', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- yvelmean ------------------ + write_netcdf_2d_state_vars('yvelmean', 'yvelmean', + 'land_ice_vertical_mean_y_velocity', + 'm s-1', 'Mean velocity in y', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- litemptop ------------------ + write_netcdf_2d_state_vars('surfaceTemperature', 'litemptop', + 'temperature_at_top_of_ice_sheet_model', 'K', + 'Surface temperature', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- litempbotgr ------------------ + write_netcdf_2d_state_vars('litempbotgr', 'litempbotgr', + 'temperature_at_base_of_ice_sheet_model', 'K', + 'Basal temperature beneath grounded ice sheet', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- litempbotfl ------------------ + write_netcdf_2d_state_vars('litempbotfl', 'litempbotfl', + 'temperature_at_base_of_ice_sheet_model', 'K', + 'Basal temperature beneath floating ice shelf', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- strbasemag ------------------ + write_netcdf_2d_state_vars('strbasemag', 'strbasemag', + 'land_ice_basal_drag ', 'Pa', + 'Basal drag', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- sftgif ------------------ + write_netcdf_2d_state_vars('sftgif', 'sftgif', + 'land_ice_area_fraction', '1', + 'Land ice area fraction', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- sftgrf ------------------ + write_netcdf_2d_state_vars('sftgrf', 'sftgrf', + 'grounded_ice_sheet_area_fraction', '1', + 'Grounded ice sheet area fraction', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) + + # ----------- sftflf ------------------ + write_netcdf_2d_state_vars('sftflf', 'sftflf', + 'floating_ice_shelf_area_fraction', '1', + 'Floating ice shelf area fraction', + file_remapped_mali_state, + ismip7_grid_file, output_path, metadata) diff --git a/landice/output_processing_li/ismip7_postprocessing/validate.py b/landice/output_processing_li/ismip7_postprocessing/validate.py new file mode 100644 index 000000000..f7482852f --- /dev/null +++ b/landice/output_processing_li/ismip7_postprocessing/validate.py @@ -0,0 +1,84 @@ +""" +Generic validation utilities for sets of MALI output files. +""" + +import os +import xarray as xr + + +def validate_mali_files(files, required_vars, label=''): + """ + Validate a sorted list of MALI output files before processing. + + Checks that: + - The list is not empty + - Each file exists + - Each file contains all required variables + - simulationStartTime is consistent across all files + - No time overlaps (daysSinceStart) exist between consecutive files + - No unexpectedly large time gaps (> 366 days) exist between + consecutive files + + Parameters + ---------- + files : list of str + Sorted list of file paths. + required_vars : list of str + Variable names that must be present in every file. + label : str, optional + Human-readable label for the file type, used in messages. + + Raises + ------ + ValueError + If any validation check fails. + FileNotFoundError + If any file does not exist. + """ + tag = f' ({label})' if label else '' + + if len(files) == 0: + raise ValueError(f"No files provided{tag}.") + + for f in files: + if not os.path.exists(f): + raise FileNotFoundError(f"File not found{tag}: {f}") + + # Check required variables in each file + for f in files: + with xr.open_dataset(f, decode_cf=False) as ds: + missing = [v for v in required_vars if v not in ds] + if missing: + raise ValueError( + f"File '{f}'{tag} is missing expected variables: {missing}") + + # Check simulationStartTime consistency across all files + if len(files) > 1: + start_times = [] + for f in files: + with xr.open_dataset(f, decode_cf=False) as ds: + start_times.append( + ds['simulationStartTime'].values.tobytes() + .decode('utf-8').strip().strip('\x00')) + if len(set(start_times)) > 1: + raise ValueError( + f"Inconsistent simulationStartTime across files{tag}: " + f"{set(start_times)}") + + # Check for time overlaps or large gaps between consecutive files + for i in range(len(files) - 1): + with xr.open_dataset(files[i], decode_cf=False) as ds_a: + end_a = float(ds_a['daysSinceStart'].values[-1]) + with xr.open_dataset(files[i + 1], decode_cf=False) as ds_b: + start_b = float(ds_b['daysSinceStart'].values[0]) + if start_b <= end_a: + raise ValueError( + f"Time overlap detected between files{tag}:\n" + f" {files[i]} (ends at day {end_a})\n" + f" {files[i + 1]} (starts at day {start_b})") + gap_days = start_b - end_a + if gap_days > 366: + print(f"WARNING: Gap of {gap_days:.1f} days between files{tag}:\n" + f" {files[i]}\n {files[i + 1]}") + + print(f"Validated {len(files)} file(s){tag}.") diff --git a/mesh_tools/create_SCRIP_files/create_SCRIP_file_from_planar_rectangular_grid_latlon.py b/mesh_tools/create_SCRIP_files/create_SCRIP_file_from_planar_rectangular_grid_latlon.py new file mode 100755 index 000000000..0f4627ab9 --- /dev/null +++ b/mesh_tools/create_SCRIP_files/create_SCRIP_file_from_planar_rectangular_grid_latlon.py @@ -0,0 +1,74 @@ +#!/usr/bin/env python +# Create a SCRIP file from a planar rectanfular mesh. +# See for details: http://www.earthsystemmodeling.org/esmf_releases/public/ESMF_5_2_0rp1/ESMF_refdoc/node3.html#SECTION03024000000000000000 + +import sys +import netCDF4 +import numpy as np +from optparse import OptionParser + + +print ("== Gathering information. (Invoke with --help for more details. All arguments are optional)") +parser = OptionParser() +parser.description = "This script takes an MPAS grid file and generates a SCRIP grid file." +parser.add_option("-i", "--input", dest="inputFile", help="input grid file name used as input.", default="input.nc", metavar="FILENAME") +parser.add_option("-s", "--scrip", dest="scripFile", help="SCRIP grid file to output.", default="scrip.nc", metavar="FILENAME") + +for option in parser.option_list: + if option.default != ("NO", "DEFAULT"): + option.help += (" " if option.help else "") + "[default: %default]" +options, args = parser.parse_args() + +if not options.inputFile: + sys.exit('Error: Data input grid file is required. Specify with -c command line argument.') +if not options.scripFile: + sys.exit('Error: SCRIP output grid file is required. Specify with -s command line argument.') +print ('') # make a space in stdout before further output + +# =================================== + +fin = netCDF4.Dataset(options.inputFile, 'r') +fout = netCDF4.Dataset(options.scripFile, 'w') # This will clobber existing files + +# Get info from input file +nx = len(fin.dimensions['x']) +ny = len(fin.dimensions['y']) + +# Write to output file +# Dimensions +fout.createDimension("grid_size", nx * ny) +fout.createDimension("grid_corners", 4 ) + +print('grid rank is 2') +fout.createDimension("grid_rank", 2) + +# Variables +grid_center_lat = fout.createVariable('grid_center_lat', 'f8', ('grid_size',)) +grid_center_lat.units = 'degrees' +grid_center_lon = fout.createVariable('grid_center_lon', 'f8', ('grid_size',)) +grid_center_lon.units = 'degrees' +grid_corner_lat = fout.createVariable('grid_corner_lat', 'f8', ('grid_size', 'grid_corners')) +grid_corner_lat.units = 'degrees' +grid_corner_lon = fout.createVariable('grid_corner_lon', 'f8', ('grid_size', 'grid_corners')) +grid_corner_lon.units = 'degrees' +grid_imask = fout.createVariable('grid_imask', 'i4', ('grid_size',)) +grid_imask.units = 'unitless' +grid_dims = fout.createVariable('grid_dims', 'i4', ('grid_rank',)) + + +grid_center_lat[:] = fin.variables['lat'][:].flatten() +grid_center_lon[:] = fin.variables['lon'][:].flatten() + +lat_bnds = fin.variables['lat_bnds'][:] +grid_corner_lat[:] = lat_bnds.reshape(-1, lat_bnds.shape[-1]) +lon_bnds = fin.variables['lon_bnds'][:] +grid_corner_lon[:] = lon_bnds.reshape(-1, lon_bnds.shape[-1]) + +grid_imask[:] = 1 # For now, assume we don't want to mask anything out - but eventually may want to exclude certain cells from the input mesh during interpolation + +# set the grid dimension based on the grid rank +grid_dims[:] = [nx , ny] + +fin.close() +fout.close() +print('scrip file generation complete')