diff --git a/src/CSET/loaders/verification.py b/src/CSET/loaders/verification.py index c0b136a30..627cc4b3c 100644 --- a/src/CSET/loaders/verification.py +++ b/src/CSET/loaders/verification.py @@ -124,7 +124,7 @@ def load(conf: Config): variables={ "VARNAME": field, "BASE_MODEL": base_model["name"], - "OTHER_MODEL": model["name"], + "OTHER_MODELS": model["name"], "METHOD": method, "PRESERVED_COORDS": preserved_coords, "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", @@ -140,22 +140,22 @@ def load(conf: Config): scores_timeseries_methods = _get_scores_timeseries_methods(conf) if scores_timeseries_methods: # Produce timeseries plots of scores metrics averaged over the domain for each case study. - for model, field, scores_method in itertools.product( - models[1:], conf.SURFACE_FIELDS, scores_timeseries_methods + for field, scores_method in itertools.product( + conf.SURFACE_FIELDS, scores_timeseries_methods ): yield RawRecipe( recipe=f"timeseries_surface_difference_scores_{scores_method}.yaml", variables={ "VARNAME": field, "BASE_MODEL": base_model["name"], - "OTHER_MODEL": model["name"], + "OTHER_MODELS": [model["name"] for model in models[1:]], "SUBAREA_NAME": conf.SUBAREA_NAME if conf.SELECT_SUBAREA else "", "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"]], + model_ids=[base_model["id"]] + [model["id"] for model in models[1:]], aggregation=False, ) @@ -251,7 +251,7 @@ def load(conf: Config): variables={ "VARNAME": field, "BASE_MODEL": base_model["name"], - "OTHER_MODEL": model["name"], + "OTHER_MODELS": model["name"], "PRESERVED_COORDS": ["pressure"], "AGGREGATION_MODE": "Case-study RMSE", "SUBAREA_TYPE": conf.SUBAREA_TYPE if conf.SELECT_SUBAREA else None, diff --git a/src/CSET/operators/read.py b/src/CSET/operators/read.py index 96791d383..1a9f49ea7 100644 --- a/src/CSET/operators/read.py +++ b/src/CSET/operators/read.py @@ -167,6 +167,16 @@ def read_cubes( paths = iter_maybe(file_paths) model_names = iter_maybe(model_names) + # flattens model_names if needed into one dimensional list. + if model_names != (None,): + flat = [] + for item in model_names: + if isinstance(item, list): + flat.extend(item) + else: + flat.append(item) + model_names = flat + # Check we have appropriate number of model names. if model_names != (None,) and len(model_names) != len(paths): raise ValueError( diff --git a/src/CSET/operators/scoreswrappers.py b/src/CSET/operators/scoreswrappers.py index 681f7f5be..499de1695 100644 --- a/src/CSET/operators/scoreswrappers.py +++ b/src/CSET/operators/scoreswrappers.py @@ -41,43 +41,45 @@ logger = logging.getLogger(__name__) -def _sort_cubes_for_verification(cubes: CubeList): - """Prepare cubes ready for verification in scores. +def _sort_cube_into_base_and_other(cubes): + """Sorts cube into base and other models. Parameters ---------- cubes: iris.cube.CubeList - A CubeList of exact 2 cubes, one from each model. + A CubeList of multiple cubes. One base cube and other model cubes. Returns ------- base: iris.cube.Cube The cube from the "analysis" in the same format as the other model. - other: iris.cube.Cube - The cube from the model in the same format as the base model. + others: iris.cube.CubeList + The cube list of containing the cube(s) from the model in the same format as the base model. - Raises - ------ - ValueError: "cubes should contain exactly 2 cubes." - If any other number of cubes are present. - - Notes - ----- - This operator is used for sorting the data into the correct format. It - is likely going to need to be refactored out of CSET and perhaps moved into - `CSET._utils` given common code between here and `misc.difference`. """ - # Set cubes into correct format using code from difference operator - if len(cubes) != 2: - raise ValueError("cubes should contain exactly 2 cubes.") base: Cube = cubes.extract_cube(iris.AttributeConstraint(cset_comparison_base=1)) - other: Cube = cubes.extract_cube( + others: CubeList = cubes.extract( iris.Constraint( cube_func=lambda cube: "cset_comparison_base" not in cube.attributes ) ) - # If cubes contain a pressure coordinate, ensure it is increasing. + return base, others + + +def _ensure_increasing_pressure_coordinates(cubes): + """Ensure the pressure coordinate is increasing. + + Parameters + ---------- + cubes: iris.cube.CubeList + A CubeList of n cubes + + Returns + ------- + Cubes: iris.cube.CubeList + The original cube list but where each cube is ensured to have an increasing pressure coordinate. + """ for cube in cubes: try: if len(cube.coord("pressure").points) > 2 and not is_increasing( @@ -88,7 +90,38 @@ def _sort_cubes_for_verification(cubes: CubeList): except iris.exceptions.CoordinateNotFoundError: pass + +def _process_cubes_for_verification(base: Cube, other: Cube): + """Prepare cubes ready for verification in scores. + + Parameters + ---------- + cubes: iris.cube.CubeList + A CubeList of exact 2 cubes, one from each model. + + Returns + ------- + base: iris.cube.Cube + The cube from the "analysis" in the same format as the other model. + other: iris.cube.Cube + The cube from the model in the same format as the base model. + + Raises + ------ + ValueError: "cubes should contain exactly 2 cubes." + If any other number of cubes are present. + + Notes + ----- + This operator is used for sorting the data into the correct format. It + is likely going to need to be refactored out of CSET and perhaps moved into + `CSET._utils` given common code between here and `misc.difference`. + """ + # Set cubes into correct format using code from difference operator + # Extract just common time points. + other_model_name = other.attributes["model_name"] + base, other = _extract_common_time_points(base, other) # Get spatial coord names. @@ -144,6 +177,8 @@ def _sort_cubes_for_verification(cubes: CubeList): # Equalise attributes so we can merge. fully_equalise_attributes(CubeList([base, other])) + + other.attributes["model_name"] = other_model_name logger.debug("Base: %s\nOther: %s", base, other) return base, other @@ -219,8 +254,8 @@ def scores_rmse_model_obs( Returns ------- - scores_cube: iris.cube.Cube - A cube containing the RMSE between the models and observation cube. + scores_cubelist: iris.cube.CubeList + A cubelist containing the RMSE between the models and observation cube(s). """ rmse_cubes = CubeList() model_list = CubeList() @@ -268,68 +303,77 @@ def scores_rmse( ["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 use `None`. The default is `None`. + obs_model_comparison: bool, default False + Set true if doing model-obs comparison. Returns ------- - scores_cube: iris.cube.Cube - A cube containing the RMSE between the base and other cube. + scores_cubelist: iris.cube.CubeList + A cubelist containing the RMSE between the base and other cube. """ + scores_cubelist = CubeList() if obs_model_comparison: for cb in cubes: if "observed" in cb.long_name: base = cb else: - other = cb + others = [cb] else: - base, other = _sort_cubes_for_verification(cubes) - - # Copy the coordinates of the input 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 operates on xarray data arrays, so we transform the iris cube into an array, - # apply scores, and then transform it back. - scores_cube = xr.DataArray.to_iris( - scores.continuous.rmse( - other_xr, - base_xr, - preserve_dims=preserve_dims, - ) - ) + base, others = _sort_cube_into_base_and_other(cubes) - # 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 scores_cube.coords("time"): - base_time = base.coord("time") - time_vals = ( - base_time.bounds.flatten() - if base_time.has_bounds() - else base_time.points + for other in others: + base, other = _process_cubes_for_verification(base, other) + + # Copy the coordinates of the input 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 operates on xarray data arrays, so we transform the iris cube into an array, + # apply scores, and then transform it back. + scores_cube = xr.DataArray.to_iris( + scores.continuous.rmse( + other_xr, + base_xr, + preserve_dims=preserve_dims, ) - t_start = float(time_vals[0]) - t_end = float(time_vals[-1]) - t_mid = 0.5 * (t_start + t_end) - - scores_cube.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(), + ) + + # 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 scores_cube.coords("time"): + base_time = base.coord("time") + time_vals = ( + base_time.bounds.flatten() + if base_time.has_bounds() + else base_time.points ) - ) - except iris.exceptions.CoordinateNotFoundError: - pass + t_start = float(time_vals[0]) + t_end = float(time_vals[-1]) + t_mid = 0.5 * (t_start + t_end) + + scores_cube.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 - scores_cube.rename(f"RMSE_of_{base.name()}") - # if preserved_coordinates == ["grid_latitude", "grid_longitude"]: - # scores_cube.add_aux_coord(time_coord) - return scores_cube + scores_cube.rename(f"RMSE_of_{base.name()}") + scores_cubelist.append(scores_cube) + + model_name = other.attributes["model_name"] + scores_cube.attributes["model_name"] = model_name + + return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = None): @@ -351,56 +395,63 @@ def scores_mae(cubes: CubeList, preserved_coordinates: list[str] | str | None = Returns ------- - scores_cube: iris.cube.Cube - A cube containing the MAE between the base and other cube. + scores_cubelist: iris.cube.CubeList + A cubelist containing the MAE between the base and other cube(s). """ - base, other = _sort_cubes_for_verification(cubes) - - # Copy the coordinates of the input 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 operates on xarray data arrays, so we transform the iris cube into an array, - # apply scores, and then transform it back. - scores_cube = xr.DataArray.to_iris( - scores.continuous.mae( - other_xr, - base_xr, - preserve_dims=preserve_dims, + base, others = _sort_cube_into_base_and_other(cubes) + scores_cubelist = CubeList() + for other in others: + base, other = _process_cubes_for_verification(base, other) + + # Copy the coordinates of the input 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 operates on xarray data arrays, so we transform the iris cube into an array, + # apply scores, and then transform it back. + scores_cube = xr.DataArray.to_iris( + scores.continuous.mae( + 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 scores_cube.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) - - scores_cube.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(), + # 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 scores_cube.coords("time"): + base_time = base.coord("time") + time_vals = ( + base_time.bounds.flatten() + if base_time.has_bounds() + else base_time.points ) - ) - except iris.exceptions.CoordinateNotFoundError: - pass + t_start = float(time_vals[0]) + t_end = float(time_vals[-1]) + t_mid = 0.5 * (t_start + t_end) + + scores_cube.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 + + scores_cube.rename(f"MAE_of_{base.name()}") + scores_cubelist.append(scores_cube) + model_name = other.attributes["model_name"] + scores_cube.attributes["model_name"] = model_name - scores_cube.rename(f"MAE_of_{base.name()}") - return scores_cube + return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist def scores_additive_bias( @@ -424,55 +475,62 @@ def scores_additive_bias( Returns ------- - scores_cube: iris.cube.Cube - A cube containing the ME between the base and other cube. + scores_cubelist: iris.cube.CubeList + A cubelist containing the ME between the base and other cube(s). """ - base, other = _sort_cubes_for_verification(cubes) - - # Copy the coordinates of the input 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 operates on xarray data arrays, so we transform the iris cube into an array, - # apply scores, and then transform it back. - scores_cube = xr.DataArray.to_iris( - scores.continuous.additive_bias( - other_xr, - base_xr, - preserve_dims=preserve_dims, + base, others = _sort_cube_into_base_and_other(cubes) + scores_cubelist = CubeList() + for other in others: + base, other = _process_cubes_for_verification(base, other) + + # Copy the coordinates of the input 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 operates on xarray data arrays, so we transform the iris cube into an array, + # apply scores, and then transform it back. + scores_cube = xr.DataArray.to_iris( + scores.continuous.additive_bias( + 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 scores_cube.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) - - scores_cube.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(), + # 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 scores_cube.coords("time"): + base_time = base.coord("time") + time_vals = ( + base_time.bounds.flatten() + if base_time.has_bounds() + else base_time.points ) - ) - except iris.exceptions.CoordinateNotFoundError: - pass - scores_cube.rename(f"Additive_Bias_of_{base.name()}") - return scores_cube + t_start = float(time_vals[0]) + t_end = float(time_vals[-1]) + t_mid = 0.5 * (t_start + t_end) + + scores_cube.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 + scores_cube.rename(f"Additive_Bias_of_{base.name()}") + scores_cubelist.append(scores_cube) + model_name = other.attributes["model_name"] + scores_cube.attributes["model_name"] = model_name + + return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist def scores_correlation_pearsonr( @@ -496,56 +554,62 @@ def scores_correlation_pearsonr( Returns ------- - scores_cube: iris.cube.Cube - A cube containing the PC between the base and other cube. + scores_cubelist: iris.cube.CubeList + A cubelist containing the PC between the base and other cube(s). """ - base, other = _sort_cubes_for_verification(cubes) - - # Copy the coordinates of the input 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 operates on xarray data arrays, so we transform the iris cube into an array, - # apply scores, and then transform it back. - scores_cube = xr.DataArray.to_iris( - scores.continuous.correlation.pearsonr( - other_xr, - base_xr, - preserve_dims=preserve_dims, + base, others = _sort_cube_into_base_and_other(cubes) + scores_cubelist = CubeList() + for other in others: + base, other = _process_cubes_for_verification(base, other) + + # Copy the coordinates of the input 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 operates on xarray data arrays, so we transform the iris cube into an array, + # apply scores, and then transform it back. + scores_cube = xr.DataArray.to_iris( + scores.continuous.correlation.pearsonr( + 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 scores_cube.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) - - scores_cube.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(), + # 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 scores_cube.coords("time"): + base_time = base.coord("time") + time_vals = ( + base_time.bounds.flatten() + if base_time.has_bounds() + else base_time.points ) - ) - except iris.exceptions.CoordinateNotFoundError: - pass + t_start = float(time_vals[0]) + t_end = float(time_vals[-1]) + t_mid = 0.5 * (t_start + t_end) + + scores_cube.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 - scores_cube.rename(f"Pearson_Correlation_of_{base.name()}") - return scores_cube + scores_cube.rename(f"Pearson_Correlation_of_{base.name()}") + scores_cubelist.append(scores_cube) + model_name = other.attributes["model_name"] + scores_cube.attributes["model_name"] = model_name + return scores_cubelist[0] if len(scores_cubelist) == 1 else scores_cubelist def scores_crps_for_ensemble( diff --git a/src/CSET/recipes/verification/generic_level_rmse_scores_profile.yaml b/src/CSET/recipes/verification/generic_level_rmse_scores_profile.yaml index d4b2788d7..db0a735bd 100644 --- a/src/CSET/recipes/verification/generic_level_rmse_scores_profile.yaml +++ b/src/CSET/recipes/verification/generic_level_rmse_scores_profile.yaml @@ -1,5 +1,5 @@ category: Scores -title: "$VARNAME\nRMSE vertical profile between $OTHER_MODEL and $BASE_MODEL ($AGGREGATION_MODE)" +title: "$VARNAME\nRMSE vertical profile between $OTHER_MODELS 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 @@ -24,7 +24,7 @@ description: | steps: - operator: read.read_cubes file_paths: $INPUT_PATHS - model_names: [$BASE_MODEL, $OTHER_MODEL] + model_names: [$BASE_MODEL, $OTHER_MODELS] constraint: operator: constraints.combine_constraints variable_constraint: diff --git a/src/CSET/recipes/verification/surface_difference_scores_MAE.yaml b/src/CSET/recipes/verification/surface_difference_scores_MAE.yaml index cc77db011..a20e69095 100644 --- a/src/CSET/recipes/verification/surface_difference_scores_MAE.yaml +++ b/src/CSET/recipes/verification/surface_difference_scores_MAE.yaml @@ -29,7 +29,7 @@ description: | steps: - operator: read.read_cubes file_paths: $INPUT_PATHS - model_names: [$BASE_MODEL, $OTHER_MODEL] + model_names: [$BASE_MODEL, $OTHER_MODELS] constraint: operator: constraints.combine_constraints varname_constraint: diff --git a/src/CSET/recipes/verification/surface_difference_scores_RMSE.yaml b/src/CSET/recipes/verification/surface_difference_scores_RMSE.yaml index 9b4642cf4..65d7b5e42 100644 --- a/src/CSET/recipes/verification/surface_difference_scores_RMSE.yaml +++ b/src/CSET/recipes/verification/surface_difference_scores_RMSE.yaml @@ -26,7 +26,7 @@ description: | steps: - operator: read.read_cubes file_paths: $INPUT_PATHS - model_names: [$BASE_MODEL, $OTHER_MODEL] + model_names: [$BASE_MODEL, $OTHER_MODELS] constraint: operator: constraints.combine_constraints varname_constraint: diff --git a/src/CSET/recipes/verification/surface_difference_scores_additive_bias.yaml b/src/CSET/recipes/verification/surface_difference_scores_additive_bias.yaml index 9cc6b22ca..4b5d57bbb 100644 --- a/src/CSET/recipes/verification/surface_difference_scores_additive_bias.yaml +++ b/src/CSET/recipes/verification/surface_difference_scores_additive_bias.yaml @@ -1,5 +1,5 @@ category: Scores -title: "$METHOD $VARNAME\nScores Mean Error between $OTHER_MODEL and $BASE_MODEL $SUBAREA_NAME" +title: "$METHOD $VARNAME\nScores Mean Error between $OTHER_MODELS and $BASE_MODEL $SUBAREA_NAME" description: | Extracts and plots the Mean Error (also known as the Additive Bias) in $METHOD of $VARNAME. The ME is calculated based on that used in the @@ -24,7 +24,7 @@ description: | steps: - operator: read.read_cubes file_paths: $INPUT_PATHS - model_names: [$BASE_MODEL, $OTHER_MODEL] + model_names: [$BASE_MODEL, $OTHER_MODELS] constraint: operator: constraints.combine_constraints varname_constraint: diff --git a/src/CSET/recipes/verification/timeseries_surface_difference_scores_MAE.yaml b/src/CSET/recipes/verification/timeseries_surface_difference_scores_MAE.yaml index d5973fbbb..dc25e8176 100644 --- a/src/CSET/recipes/verification/timeseries_surface_difference_scores_MAE.yaml +++ b/src/CSET/recipes/verification/timeseries_surface_difference_scores_MAE.yaml @@ -1,5 +1,5 @@ category: Scores -title: "$VARNAME\nScores Mean Absolute Error timeseries between $OTHER_MODEL and $BASE_MODEL $SUBAREA_NAME" +title: "$VARNAME\nScores Mean Absolute Error timeseries between models and $BASE_MODEL $SUBAREA_NAME" description: | Extracts and plots the Mean Absolute Error in $VARNAME computed over the domain for each timestep. The MAE is calculated based on that used in the @@ -17,7 +17,7 @@ description: | steps: - operator: read.read_cubes file_paths: $INPUT_PATHS - model_names: [$BASE_MODEL, $OTHER_MODEL] + model_names: [$BASE_MODEL, $OTHER_MODELS] constraint: operator: constraints.combine_constraints varname_constraint: diff --git a/src/CSET/recipes/verification/timeseries_surface_difference_scores_RMSE.yaml b/src/CSET/recipes/verification/timeseries_surface_difference_scores_RMSE.yaml index 7915c332d..0306fee24 100644 --- a/src/CSET/recipes/verification/timeseries_surface_difference_scores_RMSE.yaml +++ b/src/CSET/recipes/verification/timeseries_surface_difference_scores_RMSE.yaml @@ -1,5 +1,5 @@ category: Scores -title: "$VARNAME\nScores RMSE timeseries between $OTHER_MODEL and $BASE_MODEL $SUBAREA_NAME" +title: "$VARNAME\nScores RMSE timeseries between models and $BASE_MODEL and $SUBAREA_NAME" description: | Extracts and plots the Root Mean Square Error in $VARNAME computed over the domain for each timestep. The RMSE is calculated based on that used in the @@ -16,7 +16,7 @@ description: | steps: - operator: read.read_cubes file_paths: $INPUT_PATHS - model_names: [$BASE_MODEL, $OTHER_MODEL] + model_names: [$BASE_MODEL, $OTHER_MODELS] constraint: operator: constraints.combine_constraints varname_constraint: diff --git a/src/CSET/recipes/verification/timeseries_surface_difference_scores_additive_bias.yaml b/src/CSET/recipes/verification/timeseries_surface_difference_scores_additive_bias.yaml index df9ee7c15..cf0d6a911 100644 --- a/src/CSET/recipes/verification/timeseries_surface_difference_scores_additive_bias.yaml +++ b/src/CSET/recipes/verification/timeseries_surface_difference_scores_additive_bias.yaml @@ -1,5 +1,5 @@ category: Scores -title: "$VARNAME\nScores Mean Error timeseries between $OTHER_MODEL and $BASE_MODEL $SUBAREA_NAME" +title: "$VARNAME\nScores Mean Error timeseries between models and $BASE_MODEL $SUBAREA_NAME" description: | Extracts and plots the Mean Error (also known as the Additive Bias) in $VARNAME computed over the domain for each timestep. The ME is calculated based on that used in the @@ -15,7 +15,7 @@ description: | steps: - operator: read.read_cubes file_paths: $INPUT_PATHS - model_names: [$BASE_MODEL, $OTHER_MODEL] + model_names: [$BASE_MODEL, $OTHER_MODELS] constraint: operator: constraints.combine_constraints varname_constraint: diff --git a/src/CSET/recipes/verification/timeseries_surface_difference_scores_correlation_pearsonr.yaml b/src/CSET/recipes/verification/timeseries_surface_difference_scores_correlation_pearsonr.yaml index 179a735f9..539fe0cbe 100644 --- a/src/CSET/recipes/verification/timeseries_surface_difference_scores_correlation_pearsonr.yaml +++ b/src/CSET/recipes/verification/timeseries_surface_difference_scores_correlation_pearsonr.yaml @@ -1,5 +1,5 @@ category: Scores -title: "$VARNAME\nScores Pearson's Correlation timeseries between $OTHER_MODEL and $BASE_MODEL $SUBAREA_NAME" +title: "$VARNAME\nScores Pearson's Correlation timeseries between models and $BASE_MODEL $SUBAREA_NAME" description: | Extracts and plots the Pearson's Correlation coefficient in $VARNAME computed over the domain for each timestep. The PC is calculated based on that used in the @@ -14,7 +14,7 @@ description: | steps: - operator: read.read_cubes file_paths: $INPUT_PATHS - model_names: [$BASE_MODEL, $OTHER_MODEL] + model_names: [$BASE_MODEL, $OTHER_MODELS] constraint: operator: constraints.combine_constraints varname_constraint: diff --git a/src/CSET/recipes/verification/timeseries_surface_rmse_scores.yaml b/src/CSET/recipes/verification/timeseries_surface_rmse_scores.yaml deleted file mode 100644 index db1ed972a..000000000 --- a/src/CSET/recipes/verification/timeseries_surface_rmse_scores.yaml +++ /dev/null @@ -1,43 +0,0 @@ -category: Scores -title: "$VARNAME\nRMSE timeseries between $OTHER_MODEL and $BASE_MODEL" -description: | - Extracts and plots the Root Mean Square Error in $VARNAME - for all times. The RMSE is calculated based on that used in the - package [`scores`](https://scores.readthedocs.io/en/stable/api.html#scores.continuous.rmse). - This recipe allows the preservation of the time coordinate to produce a timeseries. Therefore, the RMSE is - collapsed over all other coordinates in the cube and calculated for every timestep. - - A larger RMSE implies a greater error than a smaller RMSE. An - RMSE of zero implies the two fields match. The RMSE is calculated - on the grid point and thus a spatial view of the RMSE provides useful - information about where the differences are, or if placement errors - are domininating the score (usually indicated by dipoles). - -steps: - - operator: read.read_cubes - file_paths: $INPUT_PATHS - model_names: [$BASE_MODEL, $OTHER_MODEL] - constraint: - operator: constraints.combine_constraints - varname_constraint: - operator: constraints.generate_var_constraint - varname: $VARNAME - cell_methods_constraint: - operator: constraints.generate_cell_methods_constraint - cell_methods: [] - varname: $VARNAME - pressure_level_constraint: - operator: constraints.generate_level_constraint - coordinate: "pressure" - levels: [] - subarea_type: $SUBAREA_TYPE - subarea_extent: $SUBAREA_EXTENT - - - operator: scoreswrappers.scores_rmse - preserved_coordinates: ["time"] - - - operator: plot.plot_line_series - series_coordinate: time - - - operator: write.write_cube_to_nc - overwrite: True diff --git a/tests/operators/test_scoreswrappers.py b/tests/operators/test_scoreswrappers.py index a4430a27a..7c2661d67 100644 --- a/tests/operators/test_scoreswrappers.py +++ b/tests/operators/test_scoreswrappers.py @@ -38,6 +38,8 @@ def test_scores_correlation_pearsonr(cube: Cube): # Data preparation. other_cube = cube.copy() del other_cube.attributes["cset_comparison_base"] + cube.attributes["model_name"] = "model1" + other_cube.attributes["model_name"] = "model2" cubes = CubeList([cube, other_cube]) # Take difference. @@ -61,6 +63,8 @@ def test_scores_additive_bias(cube: Cube): # Data preparation. other_cube = cube.copy() del other_cube.attributes["cset_comparison_base"] + cube.attributes["model_name"] = "model1" + other_cube.attributes["model_name"] = "model2" cubes = CubeList([cube, other_cube]) # Take difference. @@ -80,6 +84,8 @@ def test_scores_mae(cube: Cube): # Data preparation. other_cube = cube.copy() del other_cube.attributes["cset_comparison_base"] + cube.attributes["model_name"] = "model1" + other_cube.attributes["model_name"] = "model2" cubes = CubeList([cube, other_cube]) # Take difference. @@ -97,6 +103,8 @@ def test_scores_rmse(cube: Cube): # Data preparation. other_cube = cube.copy() del other_cube.attributes["cset_comparison_base"] + cube.attributes["model_name"] = "model1" + other_cube.attributes["model_name"] = "model2" cubes = CubeList([cube, other_cube]) # Take difference. @@ -121,6 +129,8 @@ def test_scores_rmse_nonzero(): ) other_cube = cube.copy(data=np.ones((2, 2))) cube.attributes["cset_comparison_base"] = 1 + cube.attributes["model_name"] = "model1" + other_cube.attributes["model_name"] = "model2" different_cubes = CubeList((cube, other_cube)) # Take difference. rmse_cube = scoreswrappers.scores_rmse(different_cubes) @@ -138,6 +148,8 @@ def test_scores_rmse_no_time_coord(cube): c1.remove_coord("time") c2 = c1.copy() del c2.attributes["cset_comparison_base"] + c1.attributes["model_name"] = "model1" + c2.attributes["model_name"] = "model2" cubes = CubeList([c1, c2]) rmse_cube = scoreswrappers.scores_rmse(cubes) assert isinstance(rmse_cube, Cube) @@ -152,26 +164,13 @@ def test_scores_rmse_no_common_points(cube): new_times += 6 other_cube.coord("time").points = new_times del other_cube.attributes["cset_comparison_base"] + cube.attributes["model_name"] = "model1" + other_cube.attributes["model_name"] = "model2" cubes = CubeList([cube, other_cube]) with pytest.raises(ValueError, match="No common time points found!"): scoreswrappers.scores_rmse(cubes) -def test_scores_rmse_incorrect_number_of_cubes(cube): - """Test exception when incorrect number of cubes provided.""" - no_cubes = CubeList([]) - with pytest.raises(ValueError, match="cubes should contain exactly 2 cubes."): - scoreswrappers.scores_rmse(no_cubes) - - one_cube = CubeList([cube]) - with pytest.raises(ValueError, match="cubes should contain exactly 2 cubes."): - scoreswrappers.scores_rmse(one_cube) - - three_cubes = CubeList([cube, cube, cube]) - with pytest.raises(ValueError, match="cubes should contain exactly 2 cubes."): - scoreswrappers.scores_rmse(three_cubes) - - def test_scores_rmse_different_data_shape_regrid(cube): """Test when data shape differs, but gets regridded. @@ -180,6 +179,8 @@ def test_scores_rmse_different_data_shape_regrid(cube): rearranged_cube = cube.copy() rearranged_cube = rearranged_cube[:, :, 1:] del cube.attributes["cset_comparison_base"] + rearranged_cube.attributes["model_name"] = "model1" + cube.attributes["model_name"] = "model2" cubes = CubeList([rearranged_cube, cube]) # Need to preserve coordinates to test shape. rmse = scoreswrappers.scores_rmse( @@ -195,6 +196,8 @@ def test_rmse_grid_staggering_regrid(cube): rearranged_cube = cube.copy() rearranged_cube.rename("eastward_wind_at_10m") del cube.attributes["cset_comparison_base"] + rearranged_cube.attributes["model_name"] = "model1" + cube.attributes["model_name"] = "model2" cubes = CubeList([rearranged_cube, cube]) # Need to preserve coordinates to test shape. rmse = scoreswrappers.scores_rmse( @@ -209,6 +212,8 @@ def test_difference_different_model_types(cube): flipped = cube.copy() reverse(flipped, "grid_latitude") del flipped.attributes["cset_comparison_base"] + flipped.attributes["model_name"] = "model1" + cube.attributes["model_name"] = "model2" cubes = CubeList([cube, flipped]) # Take rmse. @@ -224,6 +229,8 @@ def test_difference_flip_pressure_order(transect_source_cube_readonly): flipped = transect_source_cube_readonly.copy() reverse(flipped, "pressure") del flipped.attributes["cset_comparison_base"] + flipped.attributes["model_name"] = "model1" + transect_source_cube_readonly.attributes["model_name"] = "model2" cubes = CubeList([transect_source_cube_readonly, flipped]) # Take rmse.