diff --git a/src/CSET/cset_workflow/includes/metplus_point_stat.cylc b/src/CSET/cset_workflow/includes/metplus_point_stat.cylc index 894c04cf0..9343880c1 100644 --- a/src/CSET/cset_workflow/includes/metplus_point_stat.cylc +++ b/src/CSET/cset_workflow/includes/metplus_point_stat.cylc @@ -103,7 +103,7 @@ # Runs VerPy metloader utility to create the VerPy databases inherit = VERPY_METLOADER [[[environment]]] - VERPY_DIR = {{VERPY_DIR}} + VERPY_DIR = {{VERPY_DIR|default("")}} VER_METHOD = area STAT_TYPE = cnt STREAM = point_stat diff --git a/src/CSET/cset_workflow/meta/verification/rose-meta.conf b/src/CSET/cset_workflow/meta/verification/rose-meta.conf index bdd307a73..11c2b930a 100644 --- a/src/CSET/cset_workflow/meta/verification/rose-meta.conf +++ b/src/CSET/cset_workflow/meta/verification/rose-meta.conf @@ -54,36 +54,71 @@ sort-key=met5 ################################### # Scores -[Verification/Scores] -ns=Verification/Scores +[Verification/Scores spatial] +ns=Verification/Scores spatial sort-key=scores -title=Scores +title=Scores spatial [template variables=SCORES_ALL] -ns=Verification/Scores +ns=Verification/Scores spatial description=If selected, all scores shown below will be calculated. NOTE scores currently only available for model-model intercomparisons. type=python_boolean trigger=template variables=SCORES_RMSE_SPATIAL: False; template variables=SCORES_RMSE_TIMESERIES: False; compulsory=true -sort-key=scores1 +sort-key=scoresspatial1 [template variables=SCORES_RMSE_SPATIAL] -ns=Verification/Scores +ns=Verification/Scores spatial description=Create spatially mapped RMSE plots for the specified surface fields. The analysis method(s) set in SPATIAL_SURFACE_FIELD_METHOD will be used. NOTE currently only available for model-model intercomparisons. type=python_boolean compulsory=true -sort-key=scores2 +sort-key=scoresspatial2 + +[Verification/Scores timeseries] +ns=Verification/Scores timeseries +sort-key=scorestimeseries +title=Scores timeseries +[Verification/Scores timeseries] [template variables=SCORES_RMSE_TIMESERIES] -ns=Verification/Scores +ns=Verification/Scores timeseries description=Create a timeseries of the RMSE for the specified surface fields. NOTE currently only available for model-model intercomparisons. type=python_boolean compulsory=true +sort-key=scorestimeseries1 + +[Verification/Scores vertical profiles] +ns=Verification/Scores vertical profiles +sort-key=scoresvertical +title=Scores vertical profiles + +[Verification/Scores vertical profiles] +[template variables=SCORES_RMSE_VERTICAL_PROFILES] +ns=Verification/Scores vertical profiles +description=Create vertical profile RMSE plots for the specified level fields. + RMSE profiles are calculated at each grid point for all timesteps within a case study at each pressure level. +type=python_boolean +compulsory=true +sort-key=scoresprofile1 + +[template variables=SCORES_RMSE_VERTICAL_PROFILES_SEQUENCE] +ns=Verification/Scores vertical profiles +description=Also produce one RMSE vertical profile per time step. +help=When True, the RMSE is calculated at each grid point and level for every + time step, collapsing only the spatial (horizontal) dimensions. This produces + a sequence of vertical RMSE profile plots, one per time step, in addition to + the profile that is aggregated across the whole case study. + + When False (default), only the aggregated vertical RMSE profile is produced. +type=python_boolean +compulsory=true +#trigger=template variables=SCORES_RMSE_VERTICAL_PROFILES: True; +sort-key=scoresprofile2 sort-key=scores3 [template variables=SCORES_CRPS_FOR_ENSEMBLE] diff --git a/src/CSET/cset_workflow/rose-suite.conf.example b/src/CSET/cset_workflow/rose-suite.conf.example index 0d892e9f4..a23a78bc7 100644 --- a/src/CSET/cset_workflow/rose-suite.conf.example +++ b/src/CSET/cset_workflow/rose-suite.conf.example @@ -133,6 +133,10 @@ SATURATION_FRACTION=False SCORES_ALL=False SCORES_RMSE_SPATIAL=False SCORES_RMSE_TIMESERIES=False +SCORES_RMSE_VERTICAL_PROFILES=False +SCORES_RMSE_VERTICAL_PROFILES_SEQUENCE=False +SCORES_CRPS_FOR_ENSEMBLE=False +SATURATION_FRACTION=False SCORES_SPATIAL_AB=False SCORES_SPATIAL_MAE=False SCORES_SPATIAL_RMSE=False diff --git a/src/CSET/cset_workflow/site/nci-gadi.cylc b/src/CSET/cset_workflow/site/nci-gadi.cylc index 543ef4ea5..3e83a92f0 100644 --- a/src/CSET/cset_workflow/site/nci-gadi.cylc +++ b/src/CSET/cset_workflow/site/nci-gadi.cylc @@ -23,7 +23,7 @@ """ [[[ environment ]]] PROJECT = {{ PROJECT }} - PYTHONPATH = "{{VERPY_DIR}}" + PYTHONPATH = "{{VERPY_DIR|default("")}}" {% if RUN_METPLUS_GRID_STAT|default(False) or RUN_METPLUS_POINT_STAT|default(False) %} [[METPLUS]] diff --git a/src/CSET/loaders/verification.py b/src/CSET/loaders/verification.py index 6d1c05a8a..77c2cc7bd 100644 --- a/src/CSET/loaders/verification.py +++ b/src/CSET/loaders/verification.py @@ -23,7 +23,10 @@ def load(conf: Config): """Yield recipes from the given workflow configuration.""" # Load a list of model detail dictionaries. models = get_models(conf.asdict()) + if not models: + return # Models are listed in order, so model 1 is the first element. + base_model = models[0] if conf.SCORES_RMSE_SPATIAL or conf.SCORES_ALL: base_model = models[0] @@ -63,6 +66,25 @@ def load(conf: Config): model_ids=[base_model["id"], model["id"]], aggregation=False, ) + # including AGGREGATION_MODE in the variables dictionary to allow for clearer labeling of plots. + if conf.SCORES_RMSE_VERTICAL_PROFILES: + for model, field in itertools.product(models[1:], conf.PRESSURE_LEVEL_FIELDS): + yield RawRecipe( + recipe="generic_level_rmse_scores_profile.yaml", + variables={ + "VARNAME": field, + "BASE_MODEL": base_model["name"], + "OTHER_MODEL": model["name"], + "PRESERVED_COORDS": ["pressure"], + "AGGREGATION_MODE": "Case-study RMSE", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT + if conf.SELECT_SUBAREA + else None, + }, + model_ids=[base_model["id"], model["id"]], + aggregation=False, + ) if conf.SCORES_CRPS_FOR_ENSEMBLE: for model, field in itertools.product(models, conf.SURFACE_FIELDS): @@ -80,3 +102,22 @@ def load(conf: Config): model_ids=[model["id"]], aggregation=False, ) + # including AGGREGATION_MODE in the variables dictionary to allow for clearer labeling of plots. + if conf.SCORES_RMSE_VERTICAL_PROFILES_SEQUENCE: + for model, field in itertools.product(models[1:], conf.PRESSURE_LEVEL_FIELDS): + yield RawRecipe( + recipe="generic_level_rmse_scores_profile.yaml", + variables={ + "VARNAME": field, + "BASE_MODEL": base_model["name"], + "OTHER_MODEL": model["name"], + "PRESERVED_COORDS": ["time", "pressure"], + "AGGREGATION_MODE": "Time-step RMSE", + "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, + "SUBAREA_EXTENT": conf.SUBAREA_EXTENT + if conf.SELECT_SUBAREA + else None, + }, + model_ids=[base_model["id"], model["id"]], + aggregation=False, + ) diff --git a/src/CSET/operators/plot.py b/src/CSET/operators/plot.py index a8a72a51e..1fe8a758b 100644 --- a/src/CSET/operators/plot.py +++ b/src/CSET/operators/plot.py @@ -2299,28 +2299,74 @@ def plot_vertical_line_series( vmin = min(x_levels) vmax = max(x_levels) - # Matching the slices (matching by seq coord point; it may happen that - # evaluated models do not cover the same seq coord range, hence matching - # necessary) - cube_iterables = _find_matched_slices(cubes, sequence_coordinate) + # Check if the cube has a sequence coordinate (e.g. time). If not, plot + # a single profile directly without iterating over a sequence. + sequence_coords = [ + cube.coord(sequence_coordinate) + for cube in cubes + if cube.coords(sequence_coordinate) + ] + has_sequence_coord = len(sequence_coords) == len(cubes) and all( + np.size(coord.points) > 1 for coord in sequence_coords + ) + has_scalar_sequence_coord = len(sequence_coords) == len(cubes) and all( + np.size(coord.points) == 1 for coord in sequence_coords + ) - # Create a plot for each value of the sequence coordinate. - # Allowing for multiple cubes in a CubeList to be plotted in the same plot for - # similar sequence values. Passing a CubeList into the internal plotting function - # for similar values of the sequence coordinate. cube_slice can be an iris.cube.Cube - # or an iris.cube.CubeList. plot_index = [] - nplot = np.size(cubes[0].coord(sequence_coordinate).points) - for cubes_slice in cube_iterables: - # Format the coordinate value in a unit appropriate way. - seq_coord = cubes_slice[0].coord(sequence_coordinate) + if has_sequence_coord: + # Matching the slices (matching by seq coord point; it may happen that + # evaluated models do not cover the same seq coord range, hence matching + # necessary) + cube_iterables = _find_matched_slices(cubes, sequence_coordinate) + nplot = np.size(cubes[0].coord(sequence_coordinate).points) + for cubes_slice in cube_iterables: + # Format the coordinate value in a unit appropriate way. + seq_coord = cubes_slice[0].coord(sequence_coordinate) + plot_title, plot_filename = _set_title_and_filename( + seq_coord, nplot, recipe_title, filename + ) + + # Do the actual plotting. + _plot_and_save_vertical_line_series( + cubes_slice, + coords, + "realization", + plot_filename, + series_coordinate, + title=plot_title, + vmin=vmin, + vmax=vmax, + ) + plot_index.append(plot_filename) + elif has_scalar_sequence_coord: + # Scalar sequence coordinate (typically aggregated time bounds): + # make one plot and include sequence period in title/filename. plot_title, plot_filename = _set_title_and_filename( - seq_coord, nplot, recipe_title, filename + sequence_coords[0], 1, recipe_title, filename ) - # Do the actual plotting. _plot_and_save_vertical_line_series( - cubes_slice, + cubes, + coords, + "realization", + plot_filename, + series_coordinate, + title=plot_title, + vmin=vmin, + vmax=vmax, + ) + plot_index.append(plot_filename) + else: + # 1D case: no sequence coordinate, plot a single profile. + plot_title = recipe_title + if filename: + plot_filename = filename + else: + plot_filename = f"{slugify(plot_title)}.png" + + _plot_and_save_vertical_line_series( + cubes, coords, "realization", plot_filename, diff --git a/src/CSET/operators/scoreswrappers.py b/src/CSET/operators/scoreswrappers.py index 26e6cc449..92ce384f3 100644 --- a/src/CSET/operators/scoreswrappers.py +++ b/src/CSET/operators/scoreswrappers.py @@ -146,6 +146,53 @@ def _sort_cubes_for_verification(cubes: CubeList): return base, other +def _resolve_preserve_dims( + cube: Cube, + data_array: xr.DataArray, + preserved_coordinates: list[str] | str | None, +) -> list[str] | None: + """Resolve preserve coordinates to xarray dimension names. + + The ``scores`` package expects preserve dimensions to match xarray + dimension names. In Iris data, commonly used coordinates such as ``time`` + may be auxiliary coordinates attached to a differently named dimension + (e.g. ``dim0``). This helper maps coordinate names to their underlying + dimension names. + """ + if preserved_coordinates is None: + return None + + coord_names = ( + [preserved_coordinates] + if isinstance(preserved_coordinates, str) + else preserved_coordinates + ) + preserve_dims: list[str] = [] + + for coord_name in coord_names: + # Already an xarray dimension name. + if coord_name in data_array.dims: + if coord_name not in preserve_dims: + preserve_dims.append(coord_name) + continue + + # Otherwise, map coordinate name to dimension index/indices. + try: + dim_indices = cube.coord_dims(coord_name) + except iris.exceptions.CoordinateNotFoundError: + # Keep original name so scores raises a clear error for unknown keys. + if coord_name not in preserve_dims: + preserve_dims.append(coord_name) + continue + + for dim_index in dim_indices: + dim_name = data_array.dims[dim_index] + if dim_name not in preserve_dims: + preserve_dims.append(dim_name) + + return preserve_dims + + def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = None): r"""Calculate the Root Mean Square Error (RMSE) using scores. @@ -160,7 +207,7 @@ def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = A CubeList containing exactly two cubes: a base and an "other" model, this can be an analysis and the model. preserved_coordinates: list[str] | str | None, default is None. - The coordinates that you wish to preserve in the calculaiton of the + The coordinates (or xarray dimension names) that you wish to preserve in the calculaiton of the RMSE. For example if you want a map of each time you can preserve ["time","grid_latitude", "grid_longitude"] or if you want a time series you can preserve ["time"], if you want to collapse to a single value @@ -187,15 +234,49 @@ def scores_rmse(cubes: CubeList, preserved_coordinates: list[str] | str | None = forecasts, predictions or models (2.5.0)". Zenodo. doi: 10.5281/zenodo.18638494 """ base, other = _sort_cubes_for_verification(cubes) + + other_xr = xr.DataArray.from_iris(other) + base_xr = xr.DataArray.from_iris(base) + preserve_dims = _resolve_preserve_dims(other, other_xr, preserved_coordinates) + # Scores operators on xarray data arrays, so we transform the iris cube into an array, # apply scores, and then transform it back. RMSE = xr.DataArray.to_iris( scores.continuous.rmse( - xr.DataArray.from_iris(other), - xr.DataArray.from_iris(base), - preserve_dims=preserved_coordinates, + other_xr, + base_xr, + preserve_dims=preserve_dims, ) ) + + # If time is aggregated out, attach a scalar time coordinate with bounds + # so plotting can display the aggregated period in the title. + try: + if not RMSE.coords("time"): + base_time = base.coord("time") + time_vals = ( + base_time.bounds.flatten() + if base_time.has_bounds() + else base_time.points + ) + t_start = float(time_vals[0]) + t_end = float(time_vals[-1]) + t_mid = 0.5 * (t_start + t_end) + + RMSE.add_aux_coord( + iris.coords.AuxCoord( + t_mid, + standard_name=base_time.standard_name, + long_name=base_time.long_name, + var_name=base_time.var_name, + units=base_time.units, + bounds=np.array([t_start, t_end]), + attributes=base_time.attributes.copy(), + ) + ) + except iris.exceptions.CoordinateNotFoundError: + pass + RMSE.rename(f"RMSE_of_{base.name()}") return RMSE diff --git a/src/CSET/recipes/verification/generic_level_rmse_scores_profile.yaml b/src/CSET/recipes/verification/generic_level_rmse_scores_profile.yaml new file mode 100644 index 000000000..0c4a75517 --- /dev/null +++ b/src/CSET/recipes/verification/generic_level_rmse_scores_profile.yaml @@ -0,0 +1,40 @@ +category: Scores +title: "$VARNAME\nRMSE vertical profile between $OTHER_MODEL and $BASE_MODEL ($AGGREGATION_MODE)" +description: | + Extracts and plots the Root Mean Square Error (RMSE) of $VARNAME + as a vertical profile on pressure levels. The RMSE is calculated + pairwise between $BASE_MODEL and $OTHER_MODEL at each grid point + and vertical level, then aggregated over horizontal grid points. + + Aggregation mode: + - Case-study RMSE: Single RMSE profile calculated using all time points in the case study. + - Time-step RMSE: RMSE profiles calculated at each time point in a case study. + + A larger RMSE implies a greater error than a smaller RMSE. An + RMSE of zero implies the two fields match exactly. + +steps: + - operator: read.read_cubes + file_paths: $INPUT_PATHS + model_names: [$BASE_MODEL, $OTHER_MODEL] + constraint: + operator: constraints.combine_constraints + variable_constraint: + operator: constraints.generate_var_constraint + varname: $VARNAME + level_constraint: + operator: constraints.generate_level_constraint + coordinate: pressure + levels: "*" + subarea_type: $SUBAREA_TYPE + subarea_extent: $SUBAREA_EXTENT + + - operator: scoreswrappers.scores_rmse + preserved_coordinates: $PRESERVED_COORDS + + - operator: plot.plot_vertical_line_series + series_coordinate: pressure + sequence_coordinate: time + + - operator: write.write_cube_to_nc + overwrite: True