diff --git a/doc/release_notes.md b/doc/release_notes.md index 5cf6f67128..60cc96fd24 100644 --- a/doc/release_notes.md +++ b/doc/release_notes.md @@ -31,6 +31,8 @@ * Change dispatch of biomass and biogas generators in CBA by (a) setting the buses' marginal prices as the generators' marginal costs and (b) removing energy budget constraints ([#719](https://github.com/open-energy-transition/open-tyndp/pull/719)). +* Add CBA per horizon summary plots for each indicator benchmarking TYNDP and Open-TYNDP ([#753](https://github.com/open-energy-transition/open-tyndp/pull/753)). + **Bugfixes and Compatibility** * Rename bus for `t339` project (Tyrrhenian) from ITSI to ITVI ([#751](https://github.com/open-energy-transition/open-tyndp/pull/751)). diff --git a/rules/cba.smk b/rules/cba.smk index c90db6d1f6..19f2f67e6b 100644 --- a/rules/cba.smk +++ b/rules/cba.smk @@ -556,6 +556,36 @@ rule summarize_indicators_per_project: "../scripts/cba/summarize_indicators.py" +def summary_benchmark_indicators(w): + """ + Returns Indicator CSVs as inputs for the per-horizon summary benchmark plot. + If collection scenarios, returns the weighted-average ensemble indicators CSV as inputs for plotting. + """ + run = get_run_name(w) + if run in cba_collection_scenarios(w): + return expand( + rules.average_indicators_per_project_and_planning_horizon.output.indicators, + planning_horizons=[w.planning_horizons], + cba_project=cba_projects(w), + run=[run], + ) + return expand( + rules.collect_indicators.output.indicators, + planning_horizons=[w.planning_horizons], + run=[run], + ) + + +rule plot_summary_projects_benchmark: + input: + indicators=summary_benchmark_indicators, + output: + plot_file=RESULTS + + "cba/ensemble_plots/summary_benchmark_{planning_horizons}.png", + script: + "../scripts/cba/plot_benchmark_indicators.py" + + rule summarize_all_indicators: input: indicators=lambda w: expand( @@ -691,6 +721,13 @@ def collect_cba_scenario_inputs(w): run=cba_scenarios(w), ) ) + inputs.extend( + expand( + rules.plot_summary_projects_benchmark.output.plot_file, + planning_horizons=config_provider("cba", "planning_horizons")(w), + run=cba_scenarios(w), + ) + ) run = get_run_name(w) if run in cba_collection_scenarios(w): @@ -743,6 +780,13 @@ def cba_ensemble_inputs(w): run=runs, ) ) + inputs.extend( + expand( + rules.plot_summary_projects_benchmark.output.plot_file, + planning_horizons=config["cba"]["planning_horizons"], + run=runs, + ) + ) return inputs diff --git a/scripts/cba/plot_benchmark_indicators.py b/scripts/cba/plot_benchmark_indicators.py index d7af3c5e50..dc40924ffb 100644 --- a/scripts/cba/plot_benchmark_indicators.py +++ b/scripts/cba/plot_benchmark_indicators.py @@ -352,6 +352,152 @@ def plot_project_benchmarks( plt.close(fig) +def plot_summary_projects_benchmark( + df: pd.DataFrame, + output_path: Path, + planning_horizon: str | None = None, + area_subtitle: str | None = None, + eps: float = 1e-6, +) -> None: + """Plot one summary subplot per indicator for each horizon""" + model_df = df[df["source"] == "Open-TYNDP"] + if model_df.empty: + logger.info("No Open-TYNDP data to plot") + return + weighted_df = model_df[model_df["cyear"] == "weighted-average"] + if not weighted_df.empty: + model_df = weighted_df + + average_df = pd.concat( + [model_df, df[df["source"] == "TYNDP 2024"]], ignore_index=True + ) + + indicators = sorted( + average_df.loc[average_df["source"] == "Open-TYNDP", "indicator"] + .dropna() + .unique() + ) + project_ids = sorted( + average_df.loc[average_df["source"] == "Open-TYNDP", "project_id"] + .dropna() + .unique() + ) + + plot_items = {} + for indicator in indicators: + pairs = [] + for project_id in project_ids: + project_df = average_df[average_df["project_id"] == project_id] + if indicator == "B2a_societal_cost_variation": + model_val = select_value_by_subindex( + project_df, indicator, "Open-TYNDP", "central" + ) + bench_val = select_value_by_subindex( + project_df, indicator, "TYNDP 2024", "central" + ) + if model_val is not None and bench_val is not None: + pairs.append((bench_val, model_val)) + else: + model = benchmark_range(project_df, indicator, source="Open-TYNDP") + bench = benchmark_range(project_df, indicator, source="TYNDP 2024") + if model is not None and bench is not None: + pairs.append((bench[1], model[1])) + if pairs: + plot_items[indicator] = pairs + if not plot_items: + logger.info("Incomplete benchmark data to plot") + return + + ncols = min(4, len(plot_items)) + nrows = (len(plot_items) + ncols - 1) // ncols + fig, axes = plt.subplots( + nrows=nrows, + ncols=ncols, + figsize=(3.6 * ncols, 3.3 * nrows), + squeeze=False, + ) + + for ax, (indicator, pairs) in zip(axes.flatten(), plot_items.items()): + xs = [p[0] for p in pairs] + ys = [p[1] for p in pairs] + colors = "tab:blue" + ax.scatter( + xs, + ys, + s=20, + color=colors, + alpha=0.6, + edgecolor="white", + linewidth=0.5, + ) + ax.axline( + (0, 0), slope=1, color="black", linestyle="--", linewidth=1, alpha=0.5 + ) + + combined = [abs(v) for v in (*xs, *ys) if v != 0] + if combined: + abs_max = max(combined) + abs_min = min(combined) + if abs_min > 0 and abs_max / abs_min >= 1e3: + linthresh = max(abs_min, 1.0) + ax.set_xscale("symlog", linthresh=linthresh) + ax.set_yscale("symlog", linthresh=linthresh) + + units = average_df.loc[ + (average_df["indicator"] == indicator) + & (average_df["source"] == "Open-TYNDP"), + "units", + ].dropna() + unit_label = units.iloc[0] if not units.empty else "" + ax.set_title( + f"{indicator} ({unit_label})" if unit_label else indicator, fontsize=9 + ) + ax.set_xlabel("TYNDP 2024") + ax.set_ylabel("Open-TYNDP") + ax.axhline(0, color="gray", linewidth=0.5, alpha=0.4) + ax.axvline(0, color="gray", linewidth=0.5, alpha=0.4) + ax.grid(alpha=0.3) + + benchmark_df = pd.DataFrame(pairs, columns=["TYNDP 2024", "Open-TYNDP"]) + errors = ( + (benchmark_df["Open-TYNDP"] - benchmark_df["TYNDP 2024"]).abs() + / ( + (benchmark_df["Open-TYNDP"].abs() + benchmark_df["TYNDP 2024"].abs()) + / 2 + + eps + ) + * 100 + ) + values = ( + f"n = {len(pairs)}\n" + f"sMAPE = {errors.mean():.1f}%\n" + f"sMdAPE = {errors.median():.1f}%" + ) + ax.text( + 0.05, + 0.95, + values, + transform=ax.transAxes, + fontsize=6, + verticalalignment="top", + bbox=dict(boxstyle="round", facecolor="white", alpha=0.8), + ) + + for ax in axes.flatten()[len(plot_items) :]: + ax.axis("off") + + title = "Benchmark of indicators across all projects" + if planning_horizon: + title += f" ({planning_horizon})" + fig.suptitle(title, y=0.995) + if area_subtitle: + fig.text(0.5, 0.965, area_subtitle, ha="center", va="center", fontsize=9) + + fig.tight_layout(rect=[0, 0, 1, 0.94]) + fig.savefig(output_path, dpi=400) + plt.close(fig) + + def create_plots( indicators_file: str | Path, output_path: str | Path, @@ -435,6 +581,20 @@ def create_plots( set_scenario_config(snakemake) planning_horizon = snakemake.wildcards.get("planning_horizons") - output_target = snakemake.output.get("plot_file") or snakemake.output.plot_dir area = snakemake.config.get("cba", {}).get("area") - create_plots(snakemake.input.indicators, output_target, planning_horizon, area) + + if "cba_project" in snakemake.wildcards.keys(): + output_target = snakemake.output.get("plot_file") or snakemake.output.plot_dir + create_plots(snakemake.input.indicators, output_target, planning_horizon, area) + elif not snakemake.input.indicators: + logger.warning( + "No indicators input files for summary plot: %s", snakemake.output.plot_file + ) + else: + df = pd.concat(map(pd.read_csv, snakemake.input.indicators), ignore_index=True) + plot_summary_projects_benchmark( + df, + Path(snakemake.output.plot_file), + planning_horizon, + format_area_subtitle(area), + )