diff --git a/README.md b/README.md index bb0d400..f8aa5f2 100644 --- a/README.md +++ b/README.md @@ -5,7 +5,7 @@ Checks ISMIP7 NetCDF simulation datasets for compliance with the [ISMIP7 data re 1. **Naming** — variable name, region field, ISM member id (`mNNN`), ESM name (CMIP6/CMIP7 registry), forcing member id (`fNNN`), set counter (`[C|E|P]NNN`), and year range (`YYYY-YYYY` matching the actual time axis). 2. **Numerical** — units match the data request; all values lie within the allowed min/max range for the relevant region; array is not entirely fill values. 3. **Spatial** *(xyt variables only)* — grid corners lie within the expected AIS or GrIS extents; resolution is one of the allowed values; x and y cell size are equal. -4. **Time** — time dimension is present, unlimited, and monotonically increasing; annual cadence; experiment start/end dates and duration match `experiments_ismip7.csv`, accounting for the variable type (ST or FL, see [Time encoding](#time-encoding)). +4. **Time** — time dimension is present, unlimited, and monotonically increasing; annual cadence for `x,y,t` and `t` variables; sparse snapshot timestamps for `x,y,z,t` variables (see [Time encoding](#time-encoding)); experiment start/end dates and duration match `experiments_ismip7.csv` for `x,y,t` and `t` variables. 5. **Attributes** — required global and coordinate attributes are present and have correct values; `standard_name` matches data request; `_FillValue` equals the NetCDF4 default for the variable's dtype; variable and time are float32; `scale_factor` and `add_offset` are not allowed. Compliance criteria are defined in `conventions/ISMIP7_variable_request.csv` (variable metadata) and `experiments_ismip7.csv` (valid experiment year ranges and durations). @@ -27,6 +27,17 @@ For example, a 286-year `ctrl` simulation covering nominal years 2015–2300: The filename year range (`YYYY-YYYY`) always refers to the **nominal simulation years** (2015–2300 in the example above), regardless of variable type. The checker accounts for this when validating the filename against the time axis. +### Snapshot variables (`x,y,z,t`, e.g. `litemp`) + +`x,y,z,t` variables carry a sparse set of ST snapshots rather than a full annual time series. The required snapshot nominal years depend on the experiment type: + +| Experiment | Required snapshots | +|---|---| +| `historical` | first year of run, 1900 (if in range), 2000 (if in range), last year of run | +| projection (e.g. `ssp585`, `ctrl`) | 2100, 2200, 2300 (each if within the experiment's year range) | + +Together, a `historical` run ending at 2014 and a projection starting at 2015 provide snapshots at 100-year intervals (1900, 2000, 2100, 2200, 2300) as well as at the first and last years of the historical run. The filename year range for `litemp` reflects the full simulation year range (e.g. `2015-2300`), not the first/last snapshot year. Duration and start/end timestamp checks are not applied to snapshot variables. + Reference lookup tables are available in the companion repository [`ismip7-time-encoding`](https://github.com/ismip/ismip7-time-encoding). --- @@ -48,18 +59,18 @@ The script must be run from the repository root. It writes `compliance_checker_l ```bash # Check x,y,t (3D spatial) variables -python compliance_checker.py --source-path ./Models/GrIS/ISMIP7/SYNTH1/CORE --variable-list ismip7_xyt +python compliance_checker.py --source-path ./Models/GrIS/ISMIP7/SYNTH1/CORE/C001 --variable-list ismip7_xyt # Check scalar (time-only) variables -python compliance_checker.py --source-path ./Models/AIS/ISMIP7/SYNTH1/CORE --variable-list ismip7_scalars +python compliance_checker.py --source-path ./Models/AIS/ISMIP7/SYNTH1/CORE/C001 --variable-list ismip7_scalars # Check both -python compliance_checker.py --source-path ./Models/GrIS/ISMIP7/SYNTH1/CORE --variable-list ismip7 +python compliance_checker.py --source-path ./Models/GrIS/ISMIP7/SYNTH1/CORE/C001 --variable-list ismip7 ``` | Option | Default | Description | |--------|---------|-------------| -| `--source-path` | `./Models/GrIS/ISMIP7/SYNTH1/CORE` | Directory containing `.nc` files to check | +| `--source-path` | `./Models/GrIS/ISMIP7/SYNTH1/CORE/C001` | Set-counter subdirectory containing `.nc` files to check | | `--variable-list` | `ismip7_scalars` | `ismip7_xyt`, `ismip7_scalars`, or `ismip7` (both) | `experiments_ismip7.csv` defines the allowed nominal year ranges and durations for each experiment. The checker derives the expected FL and ST timestamps from these year values at runtime (see [Time encoding](#time-encoding)). diff --git a/compliance_checker.py b/compliance_checker.py index 85fad6b..44a149c 100755 --- a/compliance_checker.py +++ b/compliance_checker.py @@ -50,7 +50,7 @@ from tqdm import tqdm -DEFAULT_SOURCE_PATH = "./Models/GrIS/ISMIP7/SYNTH1/CORE" +DEFAULT_SOURCE_PATH = "./Models/GrIS/ISMIP7/SYNTH1/CORE/C001" DEFAULT_VARIABLE_LIST = "ismip7_scalars" VARIABLE_LIST_CHOICES = ("ismip7_scalars", "ismip7_xyt", "ismip7") @@ -219,7 +219,7 @@ def _load_criteria(workdir: str, variable_list: str): df = df.dropna(subset=["Variable Name"]) if variable_list == "ismip7_xyt": - df = df[df["Dim"] == "x,y,t"] + df = df[df["Dim"].isin(["x,y,t", "x,y,z,t", "x,y"])] elif variable_list == "ismip7_scalars": df = df[df["Dim"] == "t"] @@ -739,9 +739,11 @@ def _check_naming( ) errors += 1 + is_static = not ({"time", "t"} & dim) year_range_field = parts[ISMIP7_FILENAME_YEAR_RANGE_IDX].removesuffix(".nc") - year_range_match = re.fullmatch(r"(\d{4})-(\d{4})", year_range_field) - if not year_range_match: + if is_static: + log_file.write(f" - Filename year range: N/A (static spatial variable)\n") + elif not (year_range_match := re.fullmatch(r"(\d{4})-(\d{4})", year_range_field)): log_file.write( f" - ERROR: year range '{year_range_field}' (field {ISMIP7_FILENAME_YEAR_RANGE_IDX}) does not match expected format YYYY-YYYY (e.g. 2015-2300).\n" ) @@ -760,17 +762,23 @@ def _check_naming( _yr_offset = 1 if var_type == "ST" else 0 actual_start = min(_decoded_time).item().year - _yr_offset actual_end = max(_decoded_time).item().year - _yr_offset - if fn_start_year != actual_start: - log_file.write( - f" - ERROR: filename start year {fn_start_year} does not match first time step year {actual_start}.\n" - ) - errors += 1 + is_snapshot = "z" in dim + yr_range_ok = True + if not is_snapshot: + # For regular time-series variables the filename start year must match the first time step. + if fn_start_year != actual_start: + log_file.write( + f" - ERROR: filename start year {fn_start_year} does not match first time step year {actual_start}.\n" + ) + errors += 1 + yr_range_ok = False if fn_end_year != actual_end: log_file.write( f" - ERROR: filename end year {fn_end_year} does not match last time step year {actual_end}.\n" ) errors += 1 - if fn_start_year == actual_start and fn_end_year == actual_end: + yr_range_ok = False + if yr_range_ok: log_file.write( f" - Filename year range {fn_start_year}-{fn_end_year} matches time axis: OK\n" ) @@ -915,6 +923,10 @@ def _check_time( log_file.write("TIME Tests \n") if not ({"t"}.issubset(dim) or {"time"}.issubset(dim)): + if {"x", "y"}.issubset(dim): + # Static spatial variable (x,y) — no time axis is expected. + log_file.write(" - Time axis: N/A (static spatial variable)\n") + return errors log_file.write( " - ERROR: The time dimension is missing. Time Tests have been ignored.\n" ) @@ -960,7 +972,34 @@ def _check_time( ) return errors + 1 - if len(ds["time"].values) > 1: + is_snapshot_var = "z" in dim + if is_snapshot_var: + # x,y,z,t variable (e.g. litemp): sparse snapshot time axis. + # Each time value must be an ST-encoded timestamp for one of the known snapshot nominal years. + _LITEMP_SNAPSHOT_NOMINAL_YEARS = {1900, 2000, 2100, 2200, 2300} + _yr_offset = 1 if var_type == "ST" else 0 + nominal_years = [t.year - _yr_offset for t in ds["time"].values] + exp_for_snaps = experiments[index_exp] + valid_snapshot_years = {nominal_years[-1]} | { + y for y in _LITEMP_SNAPSHOT_NOMINAL_YEARS + if exp_for_snaps["start_year_min"] <= y <= exp_for_snaps["end_year"] + } + if exp_for_snaps["experiment"] == "historical": + valid_snapshot_years.add(nominal_years[0]) + snap_errors = 0 + for ny in nominal_years: + if ny not in valid_snapshot_years: + log_file.write( + f" - ERROR: time snapshot nominal year {ny} is not in the valid" + f" snapshot set {sorted(valid_snapshot_years)}.\n" + ) + snap_errors += 1 + if snap_errors == 0: + log_file.write( + f" - Snapshot time axis nominal years {nominal_years} are all valid: OK\n" + ) + errors += snap_errors + elif len(ds["time"].values) > 1: if isinstance(ds["time"].values[1] - ds["time"].values[0], datetime.timedelta): time_step = (ds["time"].values[1] - ds["time"].values[0]).days elif isinstance(ds["time"].values[1] - ds["time"].values[0], np.timedelta64): @@ -991,6 +1030,12 @@ def _check_time( end_exp.item().year, end_exp.item().month, end_exp.item().day ) + if is_snapshot_var: + log_file.write( + " - Duration / start / end timestamp checks: N/A (snapshot variable — time axis holds sparse snapshots only).\n" + ) + return errors + # Compute expected timestamp windows from nominal years + var_type (±1 day tolerance). _tol = datetime.timedelta(days=1) startinf = _nominal_to_timestamp(exp["start_year_min"], var_type) - _tol @@ -1012,14 +1057,15 @@ def _check_time( ) errors += 1 + _yr_off = 1 if var_type == "ST" else 0 dateformat_start_exp = datetime.datetime( start_exp.item().year, start_exp.item().month, start_exp.item().day ) if startinf <= dateformat_start_exp <= startsup: log_file.write( - " - Experiment starts correctly on " + " - First time stamp correctly set to " + start_exp.item().strftime("%Y-%m-%d") - + ".\n" + + f" (nominal year {start_exp.item().year - _yr_off}).\n" ) else: log_file.write( @@ -1035,9 +1081,9 @@ def _check_time( if endinf <= dateformat_end_exp <= endsup: log_file.write( - " - Experiment ends correctly on " + " - Last time stamp correctly set to " + end_exp.item().strftime("%Y-%m-%d") - + ".\n" + + f" (nominal year {end_exp.item().year - _yr_off}).\n" ) else: log_file.write( @@ -1053,6 +1099,7 @@ def _check_time( else: if duration_years == exp["duration"]: log_file.write(" - Experiment lasts " + str(duration_years) + " years.\n") + _yr_off = 1 if var_type == "ST" else 0 dateformat_start_exp = datetime.datetime( start_exp.item().year, start_exp.item().month, @@ -1060,9 +1107,9 @@ def _check_time( ) if startinf <= dateformat_start_exp <= startsup: log_file.write( - " - Experiment starts correctly on " + " - First time stamp correctly set to " + start_exp.item().strftime("%Y-%m-%d") - + ".\n" + + f" (nominal year {start_exp.item().year - _yr_off}).\n" ) else: log_file.write( @@ -1078,9 +1125,9 @@ def _check_time( if endinf <= dateformat_end_exp <= endsup: log_file.write( - " - Experiment ends correctly on " + " - Last time stamp correctly set to " + end_exp.item().strftime("%Y-%m-%d") - + ".\n" + + f" (nominal year {end_exp.item().year - _yr_off}).\n" ) else: log_file.write( @@ -1157,9 +1204,12 @@ def _check_attributes( if name in ds.coords: time_coord = name break - if time_coord is None: + is_static_spatial = time_coord is None and {"x", "y"}.issubset(set(ds.coords)) + if time_coord is None and not is_static_spatial: log_file.write(" - ERROR (attributes): coordinate 'time' not found.\n") coord_errors += 1 + elif time_coord is None and is_static_spatial: + log_file.write(" - Time coordinate: N/A (static spatial variable)\n") else: # xarray decodes 'units' and 'calendar' into .encoding; 'bounds' stays in .attrs time_var = ds[time_coord] @@ -1184,7 +1234,8 @@ def _check_attributes( ) coord_errors += 1 if not isscalar: - for coord in ("x", "y"): + spatial_coords = ("x", "y", "z") if "z" in ds.coords else ("x", "y") + for coord in spatial_coords: if coord in ds.coords: if "units" not in ds[coord].attrs: log_file.write( diff --git a/conventions/ISMIP7_variable_request.csv b/conventions/ISMIP7_variable_request.csv index 22e8cd4..8d4db8f 100644 --- a/conventions/ISMIP7_variable_request.csv +++ b/conventions/ISMIP7_variable_request.csv @@ -10,31 +10,34 @@ Basal mass balance flux beneath floating ice,"x,y,t",FL,libmassbffl,land_ice_bas Ice thickness imbalance,"x,y,t",FL,dlithkdt,tendency_of_land_ice_thickness,m s-1,yes,dHdt,-0.0001,0.0001,-0.0001,0.0001 Surface velocity in x,"x,y,t",ST,xvelsurf,land_ice_surface_x_velocity,m s-1,no,u-velocity at land ice surface,-0.0004,0.0004,-0.0008,0.0008 Surface velocity in y,"x,y,t",ST,yvelsurf,land_ice_surface_y_velocity,m s-1,no,v-velocity at land ice surface,-0.0004,0.0004,-0.0008,0.0008 -Surface velocity in z,"x,y,t",ST,zvelsurf,land_ice_surface_upward_velocity,m s-1,no,w-velocity at land ice surface,-4.00E-06,4e-06,-8.00E-06,8e-06 +Surface velocity in z,"x,y,t",ST,zvelsurf,land_ice_surface_upward_velocity,m s-1,no,w-velocity at land ice surface,-4.00E-06,4.00E-06,-8.00E-06,8.00E-06 Basal velocity in x,"x,y,t",ST,xvelbase,land_ice_basal_x_velocity,m s-1,no,u-velocity at land ice base,-0.0004,0.0004,-0.0008,0.0008 Basal velocity in y,"x,y,t",ST,yvelbase,land_ice_basal_y_velocity,m s-1,no,v-velocity at land ice base,-0.0004,0.0004,-0.0008,0.0008 -Basal velocity in z,"x,y,t",ST,zvelbase,land_ice_basal_upward_velocity,m s-1,no,w-velocity at land ice base,-5.00E-06,5e-06,-8.00E-06,8e-06 +Basal velocity in z,"x,y,t",ST,zvelbase,land_ice_basal_upward_velocity,m s-1,no,w-velocity at land ice base,-5.00E-06,5.00E-06,-8.00E-06,8.00E-06 Mean velocity in x,"x,y,t",ST,xvelmean,land_ice_vertical_mean_x_velocity,m s-1,yes,Vertical mean land ice velocity ,-0.0004,0.0004,-0.0008,0.0008 Mean velocity in y,"x,y,t",ST,yvelmean,land_ice_vertical_mean_y_velocity,m s-1,yes,The vertical mean land ice velocity is the average from the bedrock to the surface of the ice,-0.0004,0.0004,-0.0008,0.0008 Surface temperature,"x,y,t",ST,litemptop,temperature_at_top_of_ice_sheet_model,K,no,Ice temperature at surface,183,290,183,290 Depth average temperature,"x,y,t",ST,litempavg,,K,no,Should get standard name land_ice_temperature,183,290,183,290 -Vertical Basal temperature gradient beneath grounded ice sheet,"x,y,t",ST,litempgradgr,,K m-1,no,"Vertical Basal temperature gradient beneath grounded ice sheet, should get a standard name",-10,10,-10,10 -Vertical Basal temperature gradient beneath floating ice shelf,"x,y,t",ST,litempgradfl,,K m-1,no,"Vertical Basal temperature gradient beneath floating ice shelf, should get a standard name",-10,10,-10,10 Basal temperature beneath grounded ice sheet,"x,y,t",ST,litempbotgr,temperature_at_base_of_ice_sheet_model,K,no,Ice temperature at base of grounded ice sheet,183,290,183,290 Basal temperature beneath floating ice shelf,"x,y,t",ST,litempbotfl,temperature_at_base_of_ice_sheet_model,K,no,Ice temperature at base of floating ice shelf,183,290,183,290 Basal drag,"x,y,t",ST,strbasemag,land_ice_basal_drag,Pa,yes,Basal drag,0,300000,0,300000 -Calving flux,"x,y,t",FL,licalvf,land_ice_specific_mass_flux_due_to_calving,kg m-2 s-1,yes,Ice mass change resulting from iceberg calving. Only for grid cells in contact with ocean,-1E+11,0,-1E+11,0 -Ice front melt flux,"x,y,t",FL,lifmassbf,,kg m-2 s-1,yes,Ice mass change resulting from ice front melting. Only for grid cells in contact with ocean. Should get standard name: land_ice_specific_mass_flux_due_to_ice_front_melting,-1E+11,0,-1E+11,0 +Calving flux,"x,y,t",FL,licalvf,land_ice_specific_mass_flux_due_to_calving,kg m-2 s-1,yes,Ice mass change resulting from iceberg calving. Only for grid cells in contact with ocean,-1.00E+11,0,-1.00E+11,0 +Grounding line flux,"x,y,t",FL,ligroundf,,kg m-2 s-1,yes,Flux of ice mass across the grounding line. Only for grounding line grid cells. Should get standard name: land_ice_specific_grounding_line_flux,0,1.00E+11,0,1.00E+11 +Ice front melt flux,"x,y,t",FL,lifmassbf,,kg m-2 s-1,yes,Ice mass change resulting from ice front melting. Only for grid cells in contact with ocean. Should get standard name: land_ice_specific_mass_flux_due_to_ice_front_melting,-1.00E+11,0,-1.00E+11,0 Land ice area fraction,"x,y,t",ST,sftgif,land_ice_area_fraction,1,yes,"Fraction of grid cell covered by land ice (ice sheet, ice shelf, ice cap, glacier)",0,1,0,1 Grounded ice sheet area fraction,"x,y,t",ST,sftgrf,grounded_ice_sheet_area_fraction,1,yes,"Fraction of grid cell covered by grounded ice sheet, where grounded indicates that the quantity correspond to the ice sheet that flows over bedrock",0,1,0,1 Floating ice sheet area fraction,"x,y,t",ST,sftflf,floating_ice_shelf_area_fraction,1,yes,Fraction of grid cell covered by ice sheet flowing over seawater,0,1,0,1 -Total ice mass,t,ST,lim,land_ice_mass,kg,yes,"spatial integration, volume times density",0,1e+25,0,1e+25 -Mass above floatation,t,ST,limnsw,land_ice_mass_not_displacing_sea_water,kg,yes,"spatial integration, volume times density",0,1e+25,0,1e+25 -Grounded ice area,t,ST,iareagr,grounded_ice_sheet_area,m^2,yes,spatial integration,0,1e+25,0,1e+25 -Floating ice area,t,ST,iareafl,floating_ice_shelf_area,m^2,yes,spatial integration,0,1e+25,0,1e+25 -Total SMB flux,t,FL,tendacabf,tendency_of_land_ice_mass_due_to_surface_mass_balance,kg s-1,yes,spatial integration,0,1e+25,0,1e+25 -Total BMB flux beneath grounded ice,t,FL,tendlibmassbfgr,tendency_of_land_ice_mass_due_to_basal_mass_balance,kg s-1,yes,spatial integration,0,1e+25,0,1e+25 -Total BMB flux beneath floating ice,t,FL,tendlibmassbffl,tendency_of_land_ice_mass_due_to_basal_mass_balance,kg s-1,yes,spatial integration (computed beneath floating ice only),0,1e+25,0,1e+25 -Total calving flux,t,FL,tendlicalvf,tendency_of_land_ice_mass_due_to_calving,kg s-1,yes,spatial integration,0,1e+25,0,1e+25 -Total ice front melting flux,t,FL,tendlifmassbf,,kg s-1,yes,"spatial integration, should get standard name: tendency_of_land_ice_mass_due_to_ice_front_melting",0,1e+25,0,1e+25 -Total grounding line flux,t,FL,tendligroundf,,kg s-1,yes,"spatial integration, should get a standard name",0,1e+25,0,1e+25 \ No newline at end of file +Thermal forcing at the ice base under floating ice shelves,"x,y,t",ST,tfbase,,K,no,Thermal forcing interpolated to the ice draft under floating ice shelves (fill value for purely grounded ice or not ice),-1,30,-1,30 +Anomaly in geopotential height from reference geoid,"x,y,t",ST,deltag,geopotential_height_anomaly,m,no,The change in geoid height should be relative to the reference geoid,-4000,4000,-4000,4000 +Reference geoid,"x,y",ST,refgeoid,geoid_height_above_reference_ellipsoid,m,no,This field is calculated with respect to the WGS84 reference ellipsoid,-500,500,-500,500 +Ice temperature,"x,y,z,t",ST,litemp,land_ice_temperature,K,no,"3D temperature snapshots at (1900 if in historical), initial state, 2014, 2100, 2200, and 2300. Output on the model native vertical coordinate. Only for models that have 3D temperature calculation",183,274,183,274 +Total ice mass,t,ST,lim,land_ice_mass,kg,yes,"spatial integration, volume times density",0,1.00E+25,0,1.00E+25 +Mass above floatation,t,ST,limnsw,land_ice_mass_not_displacing_sea_water,kg,yes,"spatial integration, volume times density",0,1.00E+25,0,1.00E+25 +Grounded ice area,t,ST,iareagr,grounded_ice_sheet_area,m^2,yes,spatial integration,0,1.00E+25,0,1.00E+25 +Floating ice area,t,ST,iareafl,floating_ice_shelf_area,m^2,yes,spatial integration,0,1.00E+25,0,1.00E+25 +Total SMB flux,t,FL,tendacabf,tendency_of_land_ice_mass_due_to_surface_mass_balance,kg s-1,yes,spatial integration,0,1.00E+25,0,1.00E+25 +Total BMB flux beneath grounded ice,t,FL,tendlibmassbfgr,tendency_of_land_ice_mass_due_to_basal_mass_balance,kg s-1,yes,spatial integration,0,1.00E+25,0,1.00E+25 +Total BMB flux beneath floating ice,t,FL,tendlibmassbffl,tendency_of_land_ice_mass_due_to_basal_mass_balance,kg s-1,yes,spatial integration (computed beneath floating ice only),0,1.00E+25,0,1.00E+25 +Total calving flux,t,FL,tendlicalvf,tendency_of_land_ice_mass_due_to_calving,kg s-1,yes,spatial integration,0,1.00E+25,0,1.00E+25 +Total ice front melting flux,t,FL,tendlifmassbf,,kg s-1,yes,"spatial integration, should get standard name: tendency_of_land_ice_mass_due_to_ice_front_melting",0,1.00E+25,0,1.00E+25 +Total grounding line flux,t,FL,tendligroundf,,kg s-1,yes,"spatial integration, should get a standard name",0,1.00E+25,0,1.00E+25 \ No newline at end of file diff --git a/generate/README.md b/generate/README.md index ed482c8..8557bd3 100644 --- a/generate/README.md +++ b/generate/README.md @@ -2,7 +2,7 @@ `generate/generate_test_files.py` creates ISMIP7-style NetCDF test files with synthetic data, one file per variable, following the naming convention and grid definitions used by the compliance checker. -Files are written to `Models/{GrIS|AIS}/ISMIP7/SYNTH1/CORE/`. +Files are written to `Models/{GrIS|AIS}/ISMIP7/SYNTH1/CORE/{set_counter}/` (default `C001`). ## Usage @@ -49,14 +49,25 @@ python generate/generate_test_files.py --grid GrIS_16000m --scenario ctrl \ # Generate both 3D and scalar variables, including non-mandatory ones python generate/generate_test_files.py --grid GrIS_16000m --scenario ctrl \ --xyt --scalars --include-non-mandatory --nyears 286 --start-year 2015 + +# Generate testdata for ismip7-scalar-processing +python generate/generate_test_files.py --grid AIS_16000m --scenario historical \ + --set-counter C001 --xyt --include-non-mandatory --nyears 1 --start-year 2014 +python generate/generate_test_files.py --grid AIS_16000m --scenario ssp585 \ + --set-counter C007 --xyt --include-non-mandatory --nyears 286 --start-year 2015 +python generate/generate_test_files.py --grid GrIS_16000m --scenario historical \ + --set-counter C001 --xyt --include-non-mandatory --nyears 55 --start-year 1960 +python generate/generate_test_files.py --grid GrIS_16000m --scenario ctrl \ + --set-counter C009 --xyt --include-non-mandatory --nyears 286 --start-year 2015 ``` ## Implemented conventions - CF-1.7 as baseline. - Time encoding: `days since 1850-01-01`, `calendar='standard'`. - - State (ST) variables: snapshot at year end (Dec 31). - - Flux (FL) variables: mid-year (Jul 1) with bounds from Jan 1 to Jan 1 next year. + - State (ST) variables: timestamp is Jan 1 of year N+1 (end-of-year snapshot). No `time_bounds`. + - Flux (FL) variables: timestamp is Jul 1 of year N (mid-year), with `time_bounds` = [Jan 1 of N, Jan 1 of N+1]. + - `x,y,z,t` variables (e.g. `litemp`): ST snapshots at a sparse set of nominal years. For `historical`: first year of run, 1900 (if in range), 2000 (if in range), last year of run. For projection scenarios: 2100, 2200, 2300 (each if within the simulation year range). The filename year range reflects the full simulation period, not the first/last snapshot year. - Single precision (`float32`) for all variables and time. - `_FillValue` and `missing_value` set to NetCDF4 default `f4` fill value. - `time` is an unlimited (record) dimension. diff --git a/generate/generate_test_files.py b/generate/generate_test_files.py index 828cf57..f05f1f8 100644 --- a/generate/generate_test_files.py +++ b/generate/generate_test_files.py @@ -122,6 +122,10 @@ def read_variable_criteria(csv_file, include_non_mandatory=False): dimensions = ['x', 'y', 't'] elif dim_str == 't': dimensions = ['t'] + elif dim_str == 'x,y,z,t': + dimensions = ['x', 'y', 'z', 't'] + elif dim_str == 'x,y': + dimensions = ['x', 'y'] else: dimensions = ['x', 'y', 't'] # Default @@ -185,7 +189,7 @@ def generate_synthetic_data(shape, min_val, max_val, dtype=np.float32, eps=1e-6) def create_netcdf_file(output_file, grid_name='GrIS_16000m', scenario='ctrl', start_year=2015, nyears=5, conventions_dir=None, include_non_mandatory=False, include_scalars=False, - include_xyt=False, output_root=None, + include_xyt=False, output_root=None, nz=5, ism_member_id='m001', esm_id='CESM2-WACCM', forcing_member_id='f001', set_counter='C001'): """ Create NetCDF files with ISMIP7 variables (one file per variable). @@ -241,9 +245,9 @@ def create_netcdf_file(output_file, grid_name='GrIS_16000m', scenario='ctrl', st # Determine output directory models_dir = Path(output_root) if output_root is not None else Path(__file__).parent.parent / 'Models' if grid_type == 'AIS': - output_dir = models_dir / 'AIS' / group / model / 'CORE' + output_dir = models_dir / 'AIS' / group / model / 'CORE' / set_counter else: # GrIS - output_dir = models_dir / 'GrIS' / group / model / 'CORE' + output_dir = models_dir / 'GrIS' / group / model / 'CORE' / set_counter output_dir.mkdir(parents=True, exist_ok=True) @@ -290,6 +294,14 @@ def create_netcdf_file(output_file, grid_name='GrIS_16000m', scenario='ctrl', st xyt_vars = {var: info for var, info in variables.items() if info['dimensions'] == ['x', 'y', 't']} + # 4D variables with vertical z axis (e.g. litemp) + xyzt_vars = {var: info for var, info in variables.items() + if info['dimensions'] == ['x', 'y', 'z', 't']} + + # Static 2D spatial variables (e.g. ref_geoid) + static_vars = {var: info for var, info in variables.items() + if info['dimensions'] == ['x', 'y']} + # Create scalar time-series variables (t dimension only) scalar_vars = {var: info for var, info in variables.items() if info['dimensions'] == ['t']} @@ -451,6 +463,151 @@ def create_netcdf_file(output_file, grid_name='GrIS_16000m', scenario='ctrl', st else: print(f"Skipping x,y,t variables (include_xyt=False); {len(xyt_vars)} variables not written") + # Process x,y,z,t variables (4D snapshots, e.g. litemp) + if include_xyt: + _SNAPSHOT_NOMINAL_YEARS = {1900, 2000, 2100, 2200, 2300} + end_year = start_year + nyears - 1 + snap_set = {end_year} | {y for y in _SNAPSHOT_NOMINAL_YEARS if start_year <= y <= end_year} + if scenario == 'historical': + snap_set.add(start_year) + snapshot_years = sorted(snap_set) + origin = datetime(1850, 1, 1).date() + + for var_name, var_info in xyzt_vars.items(): + # ST only for x,y,z,t variables per ISMIP7 spec + time_days = [float((datetime(y + 1, 1, 1).date() - origin).days) for y in snapshot_years] + time_coord = np.array(time_days, dtype=np.float32) + n_snapshots = len(snapshot_years) + + min_val = var_info.get(val_key_min) + max_val = var_info.get(val_key_max) + if min_val is None and max_val is None: + min_val, max_val = -1e6, 1e6 + else: + if min_val is None: + min_val = max_val - 1.0 + if max_val is None: + max_val = min_val + 1.0 + if min_val >= max_val: + max_val = min_val + 1.0 + + data = generate_synthetic_data((n_snapshots, nz, ny, nx), min_val, max_val) + + data_vars = { + var_name: ( + ('time', 'z', 'y', 'x'), + data.astype(np.float32), + { + 'long_name': var_info['description'], + 'units': var_info['units'], + 'standard_name': var_info['standard_name'], + } + ) + } + + z = np.arange(nz, dtype=np.float32) + coords = { + 'x': ('x', x, {'long_name': 'x-coordinate', 'units': 'm'}), + 'y': ('y', y, {'long_name': 'y-coordinate', 'units': 'm'}), + 'z': ('z', z, {'long_name': 'z-coordinate', 'units': '1'}), + 'time': ('time', time_coord, {'long_name': 'time', 'units': 'days since 1850-01-01', 'calendar': 'standard'}), + } + + ds = xr.Dataset(data_vars, coords=coords) + ds.attrs.update({ + 'title': f'ISMIP7 synthetic data - {var_name}', + 'history': f'Generated on {datetime.now().isoformat()}', + 'Conventions': 'CF-1.7', + 'grid_type': grid_type, + 'grid_resolution': resolution, + 'group': group, + 'model': model, + 'contact_name': contact_names, + 'contact_email': contact_emails, + 'crs': domain_crs, + }) + + snap_time_range = f"{start_year}-{end_year}" + snap_filename_template = ( + f"{domain_id}_{source_id}_{ism_id}_{ism_member_id}_{esm_id}_{forcing_member_id}_" + f"{experiment_id}_{set_counter}_{snap_time_range}.nc" + ) + filename = f"{var_name}_{snap_filename_template}" + output_path = output_dir / filename + + encoding = { + var_name: {'dtype': 'f4', '_FillValue': fillval}, + 'time': {'dtype': 'f4', '_FillValue': fillval}, + } + + ds.to_netcdf(output_path, unlimited_dims=('time',), encoding=encoding) + with netCDF4.Dataset(output_path, 'a') as nc: + nc[var_name].missing_value = fillval + created_files.append(str(output_path)) + + # Process static x,y variables (e.g. ref_geoid) + if include_xyt: + for var_name, var_info in static_vars.items(): + min_val = var_info.get(val_key_min) + max_val = var_info.get(val_key_max) + if min_val is None and max_val is None: + min_val, max_val = -1e6, 1e6 + else: + if min_val is None: + min_val = max_val - 1.0 + if max_val is None: + max_val = min_val + 1.0 + if min_val >= max_val: + max_val = min_val + 1.0 + + data = generate_synthetic_data((ny, nx), min_val, max_val) + + data_vars = { + var_name: ( + ('y', 'x'), + data.astype(np.float32), + { + 'long_name': var_info['description'], + 'units': var_info['units'], + 'standard_name': var_info['standard_name'], + } + ) + } + + coords = { + 'x': ('x', x, {'long_name': 'x-coordinate', 'units': 'm'}), + 'y': ('y', y, {'long_name': 'y-coordinate', 'units': 'm'}), + } + + ds = xr.Dataset(data_vars, coords=coords) + ds.attrs.update({ + 'title': f'ISMIP7 synthetic data - {var_name}', + 'history': f'Generated on {datetime.now().isoformat()}', + 'Conventions': 'CF-1.7', + 'grid_type': grid_type, + 'grid_resolution': resolution, + 'group': group, + 'model': model, + 'contact_name': contact_names, + 'contact_email': contact_emails, + 'crs': domain_crs, + }) + + # Static field: use 0000-0000 as year-range placeholder (checker skips this check) + static_filename_template = ( + f"{domain_id}_{source_id}_{ism_id}_{ism_member_id}_{esm_id}_{forcing_member_id}_" + f"{experiment_id}_{set_counter}_0000-0000.nc" + ) + filename = f"{var_name}_{static_filename_template}" + output_path = output_dir / filename + + encoding = {var_name: {'dtype': 'f4', '_FillValue': fillval}} + + ds.to_netcdf(output_path, unlimited_dims=(), encoding=encoding) + with netCDF4.Dataset(output_path, 'a') as nc: + nc[var_name].missing_value = fillval + created_files.append(str(output_path)) + # Process scalar variables if include_scalars: for var_name, var_info in scalar_vars.items(): @@ -683,7 +840,7 @@ def create_multiple_files(output_dir=None, n_files=3, conventions_dir=None, parser.add_argument( '--xyt', action='store_true', - help='Write x,y,t variables' + help='Write spatial variables: x,y,t (3D), x,y,z,t (4D snapshot), and x,y (static)' ) parser.add_argument( '--include-non-mandatory', diff --git a/tests/test_compliance_checker.py b/tests/test_compliance_checker.py index 2609866..886af24 100644 --- a/tests/test_compliance_checker.py +++ b/tests/test_compliance_checker.py @@ -44,7 +44,7 @@ def baseline_core_dir(tmp_path_factory): ) assert created_files, "Synthetic baseline generation did not create any files." - core_dir = baseline_root / "GrIS" / "ISMIP7" / "SYNTH1" / "CORE" + core_dir = baseline_root / "GrIS" / "ISMIP7" / "SYNTH1" / "CORE" / "C001" baseline_summary = compliance_checker.run_checker( source_path=str(core_dir), variable_list="ismip7_scalars", @@ -60,7 +60,7 @@ def baseline_core_dir(tmp_path_factory): @pytest.fixture def case_dir(tmp_path, baseline_core_dir): - case_root = tmp_path / "CORE" + case_root = tmp_path / "CORE" / "C001" shutil.copytree(baseline_core_dir, case_root) return case_root