diff --git a/config/config.yaml b/config/config.yaml index 4b6d3a0..113d6ac 100644 --- a/config/config.yaml +++ b/config/config.yaml @@ -32,7 +32,7 @@ trade_chains: id: supplyconstraint cost_year: 2050 final_product: steel - wacc: regional #regional or uniform + wacc: uniform #regional or uniform labour_cost: True #Include labour cost tradeable_commodities: [iron_ore, hbi] stages: @@ -114,7 +114,7 @@ scenario: modifiers: cost_penalty: mga: - slack: [0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01] + slack: [0.001, 0.002, 0.003, 0.004, 0.005, 0.006, 0.007, 0.008, 0.009, 0.01, 0.011, 0.012] sense: "min" carrier: "hbi" indicator: "blocks" @@ -171,7 +171,7 @@ region_nice_names: "South_South_America": "South Cone" "Central_America": "Central America" "West_Asia": "South Asia" - "East_Asia": "East Asia (China-dominated)" + "East_Asia": "China, Mongolia, TWN & HKG" #East Asia (China-dominated)" "Pacific_Asia": "Pacific Asia" "East_East_Asia": "Japan & Korea" "Oceania": "Australia & New Zealand" @@ -287,7 +287,9 @@ colors: hydrogen direct iron reduction furnace: 'grey' battery inverter: '#6A000E' battery inverter (charging): '#6A000E' + battery (charging): '#6A000E' battery inverter (discharging): '#FF5733' + battery (discharging): '#FF5733' battery: 'purple' wind: 'blue' onshore wind: 'blue' diff --git a/rules/preparation.smk b/rules/preparation.smk index cfba22b..de94fd5 100644 --- a/rules/preparation.smk +++ b/rules/preparation.smk @@ -29,6 +29,7 @@ rule prepare_wacc: bus_locations="data/bus_locations.csv", output: wacc="resources/wacc-clustered.csv", + wacc_latex="resources/wacc-clustered.tex", threads: 2 resources: mem_mb=5000, diff --git a/rules/reporting.smk b/rules/reporting.smk index d3565e8..ce3ee09 100644 --- a/rules/reporting.smk +++ b/rules/reporting.smk @@ -5,6 +5,8 @@ main optimization workflow. """ rule plot_regions: + params: + region_nice_names=config["region_nice_names"], output: global_map_countries = "results/figures_general/global_map_countries.pdf", #workflow/notebooks/plot_countries.ipynb global_map_countries_png = "results/figures_general/global_map_countries.png", #workflow/notebooks/plot_countries.ipynb @@ -40,10 +42,13 @@ rule collect_figures: rule get_figures: input: [ + "results/figures_general/comparison/cost_year~2050/wacc~uniform/lcox_comparison_East_Asia_East_East_Asia.pdf", + "results/figures_general/comparison/cost_year~2050/wacc~uniform/lcox_comparison_South_America_Europe.pdf", "results/figures_general/comparison/cost_year~2050/wacc~regional/lcox_comparison_East_Asia_East_East_Asia.pdf", "results/figures_general/comparison/cost_year~2050/wacc~regional/lcox_comparison_South_America_Europe.pdf", "results/figures_general/mga/chain_id~supplyconstraint/wacc~regional/mga_analysis.pdf", "results/chain_id~supplyconstraint/cost_year~2050/interone~hbi/intertwo~eaf/wacc~regional/final~steel/scenario~default/map_hbi.pdf", + "results/chain_id~supplyconstraint/cost_year~2050/interone~hbi/intertwo~eaf/wacc~regional/final~steel/scenario~default/map_ironore.pdf", "results/figures_general/global_map_countries.pdf", "results/figures_general/chokepoints/map_chokepoints.pdf", "results/chain_id~supplyconstraint/cost_year~2050/interone~hbi/intertwo~eaf/wacc~regional/final~steel/scenario~mga-chokepoints/map_hbi_0.002.pdf", @@ -56,13 +61,24 @@ rule get_figures: "results/figures_general/global_supply_curve/chain_id~supplyconstraint/cost_year~2050/uniform/default/global_supply_curve_cost_global_hbi.pdf", "results/figures_general/global_supply_curve/chain_id~supplyconstraint/cost_year~2050/regional/default/global_supply_curve_cost_average_hbi.pdf", "results/figures_general/global_supply_curve/chain_id~supplyconstraint/cost_year~2050/regional/default/global_supply_curve_cost_global_hbi.pdf", + "results/figures_general/global_supply_curve/chain_id~supplyconstraint/cost_year~2050/regional/default/supply_curve_details_global_hbi.pdf", + "results/figures_general/pull/chain_id~supplyconstraint/cost_year~2050/regional/default/magnitude_pull_hbi.pdf", + "results/figures_general/mga/map_robust.pdf", + # Supplementary + "results/figures_general/trade-today/Iron_Ore_net_flow.pdf", + "results/figures_general/trade-today/DRI-HBI_net_flow.pdf", + "results/figures_general/trade-today/Steel_net_flow.pdf", + "results/figures_general/hourly/hourly-example.pdf", ] output: [ + "results/figures_streamlined/lcox-east-asia-homogenous.pdf", + "results/figures_streamlined/lcox-south-america-homogenous.pdf", "results/figures_streamlined/lcox-east-asia.pdf", "results/figures_streamlined/lcox-south-america.pdf", "results/figures_streamlined/mga-analysis.pdf", "results/figures_streamlined/map-hbi-opti.pdf", + "results/figures_streamlined/map-ironore-opti.pdf", "results/figures_streamlined/map-countries.pdf", "results/figures_streamlined/map-chokepoints.pdf", "results/figures_streamlined/map-hbi-chokepoints.pdf", @@ -75,6 +91,13 @@ rule get_figures: "results/figures_streamlined/supply-unsorted-homo.pdf", "results/figures_streamlined/supply-sorted-hetero.pdf", "results/figures_streamlined/supply-unsorted-hetero.pdf", + "results/figures_streamlined/supply-details.pdf", + "results/figures_streamlined/magnitude-pull.pdf", + "results/figures_streamlined/map-robust.pdf", + "results/figures_streamlined/today-ironore.pdf", + "results/figures_streamlined/today-dri-hbi.pdf", + "results/figures_streamlined/today-steel.pdf", + "results/figures_streamlined/hourly.pdf", ] threads: 1 run: @@ -84,6 +107,8 @@ rule get_figures: rule plot_mga: + params: + region_nice_names=config["region_nice_names"], input: network_mga_production = "results/chain_id~{trade_chain}/cost_year~2050/interone~hbi/intertwo~eaf/wacc~{wacc}/final~steel/scenario~mga-stability-weighted/network.nc", network_mga_chokepoints = "results/chain_id~{trade_chain}/cost_year~2050/interone~hbi/intertwo~eaf/wacc~{wacc}/final~steel/scenario~mga-chokepoints/network.nc", @@ -108,19 +133,19 @@ rule plot_mga_all: rule plot_trade_today: input: - baci_folder = ancient("../data/BACI_HS22_V202601"), + baci_folder = ancient("data/BACI_HS22_V202601"), output: - iron_ore = "../results/figures_general/trade-today/Iron_Ore_net_flow.pdf", - iron_ore_png = "../results/figures_general/trade-today/Iron_Ore_net_flow.png", - dri_hbi = "../results/figures_general/trade-today/DRI-HBI_net_flow.pdf", - dri_hbi_png = "../results/figures_general/trade-today/DRI-HBI_net_flow.png", - steel_raw = "../results/figures_general/trade-today/Steel_raw_net_flow.pdf", - steel_raw_png = "../results/figures_general/trade-today/Steel_raw_net_flow.png", - iron_ore_csv = "../results/figures_general/trade-today/Iron_Ore_trade_iso3.csv", - dri_hbi_csv = "../results/figures_general/trade-today/DRI-HBI_trade_iso3.csv", - steel_raw_csv = "../results/figures_general/trade-today/Steel_raw_trade_iso3.csv", + iron_ore = "results/figures_general/trade-today/Iron_Ore_net_flow.pdf", + iron_ore_png = "results/figures_general/trade-today/Iron_Ore_net_flow.png", + dri_hbi = "results/figures_general/trade-today/DRI-HBI_net_flow.pdf", + dri_hbi_png = "results/figures_general/trade-today/DRI-HBI_net_flow.png", + steel = "results/figures_general/trade-today/Steel_net_flow.pdf", + steel_png = "results/figures_general/trade-today/Steel_net_flow.png", + iron_ore_csv = "results/figures_general/trade-today/Iron_Ore_trade_iso3.csv", + dri_hbi_csv = "results/figures_general/trade-today/DRI-HBI_trade_iso3.csv", + steel_csv = "results/figures_general/trade-today/Steel_trade_iso3.csv", script: - "notebooks/plot_todays-trade.py" + str(NOTEBOOKS_DIR / "plot_todays-trade.py") rule plot_global_supply: input: @@ -131,10 +156,14 @@ rule plot_global_supply: "resources/supply_curves/cost_year~{cost_year}/wacc~{wacc}/{region}_marginal_cost_{interone}.csv", allow_missing=True, region=config["regions"]), steel_demand="resources/steel_demand_clustered_{cost_year}.csv", + population="data/owid-population/population.csv", output: network_curve="results/figures_general/global_supply_curve/chain_id~{trade_chain}/cost_year~{cost_year}/{wacc}/{scenario}/global_supply_curve_{sort}_{interone}.pdf", network_curve_png="results/figures_general/global_supply_curve/chain_id~{trade_chain}/cost_year~{cost_year}/{wacc}/{scenario}/global_supply_curve_{sort}_{interone}.png", # supply_curves + network_curve_details="results/figures_general/global_supply_curve/chain_id~{trade_chain}/cost_year~{cost_year}/{wacc}/{scenario}/supply_curve_details_{sort}_{interone}.pdf", + network_curve_details_png="results/figures_general/global_supply_curve/chain_id~{trade_chain}/cost_year~{cost_year}/{wacc}/{scenario}/supply_curve_details_{sort}_{interone}.png", + # supply_curves notebook: str(NOTEBOOKS_DIR / "analysis-globalsupplycurve.ipynb") @@ -144,9 +173,9 @@ rule plot_global_supply_all: rule plot_comparison: input: - default = "results/chain_id~labour_2050/cost_year~2050/interone~hbi/intertwo~eaf/wacc~regional/final~steel/scenario~default/network.nc", - stability = "results/chain_id~labour_2050/cost_year~2050/interone~hbi/intertwo~eaf/wacc~regional/final~steel/scenario~mga-stability-weighted/network_0.02.nc", - hightrans = "results/chain_id~hightrans_2050/cost_year~2050/interone~hbi/intertwo~eaf/wacc~regional/final~steel/scenario~default/network.nc", + default = "results/chain_id~supplyconstraint/cost_year~2050/interone~hbi/intertwo~eaf/wacc~regional/final~steel/scenario~default/network.nc", + constraint100 = "results/chain_id~supplyconstraint/cost_year~2050/interone~hbi/intertwo~eaf/wacc~regional/final~steel/scenario~constrain-supply/network_100.0.nc", + hightrans = "results/chain_id~hightrans/cost_year~2050/interone~hbi/intertwo~eaf/wacc~regional/final~steel/scenario~default/network.nc", output: cost_comparison="results/figures_general/comparison/cost_comparison.pdf", cost_comparison_png="results/figures_general/comparison/cost_comparison.png", @@ -154,6 +183,44 @@ rule plot_comparison: str(NOTEBOOKS_DIR / "compare-scenarios.ipynb") +rule plot_hourly: + input: + network="resources/lco-hbi/cost_year~2050/wacc~regional/Europe_allocated_share/network_1.nc" + output: + hourly_plot = "results/figures_general/hourly/hourly-example.pdf", + notebook: + str(NOTEBOOKS_DIR / "analysis-hourly.ipynb") + + +rule plot_magnitude_pull: + input: + steel_production = "data/brownfield-steel/country_crude_steel_production_population_per_capita_2024.csv", + supply_curves_interone = expand("resources/supply_curves/cost_year~{cost_year}/wacc~{wacc}/{region}_marginal_cost_{interone}.csv", allow_missing=True, region=config["regions"]), + steel_demand="resources/steel_demand_clustered_{cost_year}.csv", + population="data/owid-population/population.csv", + output: + magnitude_pull="results/figures_general/pull/chain_id~{trade_chain}/cost_year~{cost_year}/{wacc}/{scenario}/magnitude_pull_{interone}.pdf", + magnitude_pull_png="results/figures_general/pull/chain_id~{trade_chain}/cost_year~{cost_year}/{wacc}/{scenario}/magnitude_pull_{interone}.png", + notebook: + str(NOTEBOOKS_DIR / "plot-magnitude-pull.ipynb") + +rule plot_magnitude_pull_all: + input: + expand("results/figures_general/pull/chain_id~{trade_chain}/cost_year~{cost_year}/{wacc}/{scenario}/magnitude_pull_{interone}.pdf", trade_chain=[config["trade_chains"]["id"]], cost_year=[2050], wacc=["regional", "uniform"], interone=["hbi"], scenario=["default"], allow_missing=True) + + +rule plot_robust_map: + input: + chokepoints="results/chain_id~supplyconstraint/cost_year~2050/interone~hbi/intertwo~eaf/wacc~regional/final~steel/scenario~mga-chokepoints/network_0.002.nc", + blocs="results/chain_id~supplyconstraint/cost_year~2050/interone~hbi/intertwo~eaf/wacc~regional/final~steel/scenario~mga-blocs/network_0.001.nc", + constrain_supply="results/chain_id~supplyconstraint/cost_year~2050/interone~hbi/intertwo~eaf/wacc~regional/final~steel/scenario~constrain-supply/network_250.0.nc", + output: + robust_map_pdf="results/figures_general/mga/map_robust.pdf", + robust_map_png="results/figures_general/mga/map_robust.png", + notebook: + str(NOTEBOOKS_DIR / "plot-robust-map.ipynb") + + # Variables captured by the plot_compare_lcox input lambda (avoids two-argument lambda) # First case @@ -196,4 +263,4 @@ rule plot_compare_lcox: rule plot_compare_lcox_all: input: - expand("results/figures_general/comparison/cost_year~{cost_year}/wacc~{wacc}/lcox_comparison_" + f"{_lcox_low_cost[0]}" + "_" + f"{_lcox_high_cost[0]}" + ".pdf", cost_year=[2050], wacc=[config["trade_chains"]["wacc"]], low_cost=_lcox_low_cost, high_cost=_lcox_high_cost, allow_missing=True) + expand("results/figures_general/comparison/cost_year~{cost_year}/wacc~{wacc}/lcox_comparison_" + f"{_lcox_low_cost[0]}" + "_" + f"{_lcox_high_cost[0]}" + ".pdf", cost_year=[2050], wacc=["regional", "uniform"], low_cost=_lcox_low_cost, high_cost=_lcox_high_cost, allow_missing=True) diff --git a/workflow/notebooks/analysis-globalsupplycurve.ipynb b/workflow/notebooks/analysis-globalsupplycurve.ipynb index ad016c2..ffb7f4a 100644 --- a/workflow/notebooks/analysis-globalsupplycurve.ipynb +++ b/workflow/notebooks/analysis-globalsupplycurve.ipynb @@ -26,7 +26,7 @@ " \"plot_global_supply\",\n", " scenario=\"default\", # default, penalty-sa, penalty-ea, penalty-nwa, penalty-oc\n", " sort= \"cost_global\", #\"cost_global\", \"cost_average\"\n", - " wacc=\"uniform\",\n", + " wacc=\"regional\",\n", " cost_year=\"2050\",\n", " interone=\"hbi\",\n", " trade_chain=\"supplyconstraint\"\n", @@ -293,7 +293,7 @@ " alpha=0.6,\n", " )\n", "\n", - " ax.set_xlabel(\"Cumulative quantity (Mt)\")\n", + " ax.set_xlabel(\"Cumulative quantity in Mt/a\")\n", " ax.set_ylabel(\"Cost of HBI in €/t$_{hbi}$\")\n", " # ax.set_title(f\"Combined {process} supply curve\")\n", " ax.set_ylim(0, ylim)\n", @@ -483,6 +483,121 @@ " ylim=ylim,\n", ")" ] + }, + { + "cell_type": "markdown", + "id": "225f54ef", + "metadata": {}, + "source": [ + "### Plot country deep dive" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc82cbc3", + "metadata": {}, + "outputs": [], + "source": [ + "region_pos = {\"East_East_Asia\": (20,830),\n", + " \"East_Asia\": (355,630),\n", + " \"West_Asia\": (505,770),\n", + " \"Oceania\": (390,700),\n", + " \"Europe\": (520,890),\n", + " }" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "65773c92", + "metadata": {}, + "outputs": [], + "source": [ + "deep_dive_regions = [\"East_East_Asia\", \"East_Asia\", \"West_Asia\", \"Oceania\", \"Europe\"]\n", + "region_nice_names = snakemake.config.get(\"region_nice_names\", {})\n", + "\n", + "fig, ax = plt.subplots(figsize=(8, 4))\n", + "max_x = 0.0\n", + "demand_markers = []\n", + "europe_midpoint = None\n", + "\n", + "for region in deep_dive_regions:\n", + " region_df = df_all[df_all[\"region\"] == region].copy()\n", + " if region_df.empty:\n", + " print(f\"Skipping {region}: no supply-curve data found.\")\n", + " continue\n", + "\n", + " # Keep ordering only within each region (no cross-region ordering).\n", + " region_df = region_df.sort_values(\"lcox [EUR/t]\").reset_index(drop=True)\n", + "\n", + " x_mt = region_df[\"demand [t]\"].cumsum() / 1e6\n", + " y_cost = region_df[\"lcox [EUR/t]\"]\n", + " color = region_colors[region]\n", + " label = region_nice_names.get(region, region.replace(\"_\", \" \"))\n", + "\n", + " ax.step(x_mt, y_cost, where=\"pre\", color=color, linewidth=2.0, label=label)\n", + "\n", + " if region == \"Europe\" and len(x_mt) > 0:\n", + " i_mid = len(x_mt) // 2\n", + " europe_midpoint = (float(x_mt.iloc[i_mid]), float(y_cost.iloc[i_mid]))\n", + "\n", + " demand_mt = float(steel_demand.loc[region, \"SteelDemand_DRI_Mt\"]) if region in steel_demand.index else 0.0\n", + " if demand_mt > 0:\n", + " at_or_above = x_mt >= demand_mt\n", + " if at_or_above.any():\n", + " y_marker = y_cost.loc[at_or_above.idxmax()]\n", + " else:\n", + " y_marker = y_cost.iloc[-1]\n", + " ax.scatter(demand_mt, y_marker, color=color, s=52, edgecolor=\"black\", linewidth=0.4, zorder=6)\n", + " demand_markers.append((float(demand_mt), float(y_marker)))\n", + "\n", + " max_x = max(max_x, float(x_mt.iloc[-1]), demand_mt)\n", + "\n", + " # Include region labels\n", + " ax.annotate(\n", + " label,\n", + " xy=region_pos[region],\n", + " fontsize=8,\n", + " color=color,\n", + " )\n", + "\n", + "if demand_markers:\n", + " x_dm, y_dm = min(demand_markers, key=lambda p: p[1])\n", + " ax.annotate(\n", + " \"domestic demand\",\n", + " xy=(x_dm, y_dm),\n", + " xytext=(x_dm + -35, y_dm + -90),\n", + " fontsize=8,\n", + " color=\"black\",\n", + " arrowprops=dict(arrowstyle=\"-\", color=\"black\", lw=0.7),\n", + " )\n", + "\n", + "if europe_midpoint is not None:\n", + " x_eu, y_eu = europe_midpoint\n", + " ax.annotate(\n", + " \"domestic supply curve\",\n", + " xy=(x_eu, y_eu),\n", + " xytext=(x_eu + 20, y_eu + 100),\n", + " fontsize=8,\n", + " color=\"black\",\n", + " arrowprops=dict(arrowstyle=\"-\", color=\"black\", lw=0.7),\n", + " )\n", + "\n", + "ax.set_xlabel(\"Cumulative quantity in Mt/a\")\n", + "ax.set_ylabel(\"Cost of HBI in €/t$_{hbi}$\")\n", + "ax.set_ylim(540, 940)\n", + "\n", + "ax.set_xlim(0, 600) \n", + "\n", + "ax.grid(axis=\"y\", alpha=0.4, zorder=0)\n", + "ax.legend(title=\"Region\", frameon=False, bbox_to_anchor=(1.05, 1), loc=\"upper left\")\n", + "plt.tight_layout()\n", + "\n", + "plt.savefig(snakemake.output.network_curve_details)\n", + "plt.savefig(snakemake.output.network_curve_details_png, dpi=300)\n", + "plt.show()" + ] } ], "metadata": { diff --git a/workflow/notebooks/analysis-hourly.ipynb b/workflow/notebooks/analysis-hourly.ipynb index 4b93efa..19e0905 100644 --- a/workflow/notebooks/analysis-hourly.ipynb +++ b/workflow/notebooks/analysis-hourly.ipynb @@ -16,34 +16,24 @@ }, { "cell_type": "markdown", - "id": "bceaaa03", + "id": "9cf9b7f3", "metadata": {}, "source": [ - "### Paths" + "### Snakemake" ] }, { "cell_type": "code", "execution_count": null, - "id": "80cca193", + "id": "f6931849", "metadata": {}, "outputs": [], "source": [ - "product = \"steel\" # \"eaf\" or \"steel\" or \"eaf-grid\"\n", - "region = \"Europe\" # e.g. \"Europe\"" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "28b8bd9b", - "metadata": {}, - "outputs": [], - "source": [ - "single_region = (\n", - " f\"../../resources/lco-{product}/cost_year~2030/{region}/network_1-0partload.nc\"\n", - ")\n", - "# single_region = f\"../../resources/lcos/cost_year~2030/Europe/network_1.nc\"" + "if \"snakemake\" not in globals():\n", + " from _helpers_notebooks import mock_snakemake\n", + " snakemake = mock_snakemake(\n", + " \"plot_hourly\",\n", + " )" ] }, { @@ -61,7 +51,7 @@ "metadata": {}, "outputs": [], "source": [ - "n = pypsa.Network(single_region)" + "n = pypsa.Network(snakemake.input.network)" ] }, { @@ -76,100 +66,135 @@ }, { "cell_type": "markdown", - "id": "7896b568", + "id": "edbcf575", "metadata": {}, "source": [ - "### Plots in hourly resolution" + "### Hourly plot (clean)" ] }, { "cell_type": "code", "execution_count": null, - "id": "c81f7493", + "id": "a9247779", "metadata": {}, "outputs": [], "source": [ - "n.statistics()" - ] - }, - { - "cell_type": "markdown", - "id": "5fa2f28a", - "metadata": {}, - "source": [ - "#### Electricity balance" + "df = n.statistics.energy_balance(aggregate_time=False).loc[\n", + " :,\n", + " [\n", + " \"renewable_electricity\", # battery exchange with renewables\n", + " \"direct_reduction_furnace\",\n", + " \"electrolysis\",\n", + " \"renewable_solar\",\n", + " \"renewable_onwind\",\n", + " \"battery_elec\",\n", + " ],\n", + " \"renewable_electricity\",\n", + "]\n", + "\n", + "df = df.loc[:, pd.Timestamp(\"2050-07-21\") : pd.Timestamp(\"2050-07-26\")]\n", + "\n", + "df = (\n", + " df.reset_index(level=[\"component\", \"bus_carrier\"], drop=True)\n", + " .groupby(\"carrier\")\n", + " .sum()\n", + ")" ] }, { "cell_type": "code", "execution_count": null, - "id": "77b62c28", + "id": "281e49a2", "metadata": {}, "outputs": [], "source": [ - "n.statistics.energy_balance.iplot.area(bus_carrier=\"electricity\")" + "order = [\"renewable_onwind\",\"renewable_solar\",\"battery_elec\",\n", + " \"renewable_electricity\",\"direct_reduction_furnace\",\n", + " \"electrolysis\"]\n", + "df =df.loc[order, :]" ] }, { "cell_type": "code", "execution_count": null, - "id": "62a44662", + "id": "f1ef596d", "metadata": {}, "outputs": [], - "source": [] - }, - { - "cell_type": "markdown", - "id": "16522945", - "metadata": {}, "source": [ - "#### Electricity balance" + "colors = snakemake.config[\"colors\"]" ] }, { "cell_type": "code", "execution_count": null, - "id": "02b04d02", + "id": "4d4efb00", "metadata": {}, "outputs": [], "source": [ - "n.statistics.energy_balance.iplot.area(bus_carrier=\"hydrogen\")" - ] - }, - { - "cell_type": "markdown", - "id": "19169c8c", - "metadata": {}, - "source": [ - "#### HBI balance" + "df_plot = df.T\n", + "\n", + "df_plot.rename(\n", + " columns={\n", + " \"battery_elec\": \"battery (discharging)\",\n", + " \"electrolysis\": \"electrolysis\",\n", + " \"direct_reduction_furnace\": \"direct reduction furnace\",\n", + " \"renewable_solar\": \"solar\",\n", + " \"renewable_onwind\": \"onshore wind\",\n", + " \"renewable_electricity\": \"battery (charging)\",\n", + " },\n", + " inplace=True,\n", + ")" ] }, { "cell_type": "code", "execution_count": null, - "id": "8650abf7", + "id": "75fde3bd", "metadata": {}, "outputs": [], "source": [ - "n.statistics.energy_balance.iplot.area(bus_carrier=\"hot briquetted iron\")" + "fig, ax = plt.subplots(figsize=(9, 4))\n", + "\n", + "df_plot.plot.area(ax=ax, stacked=True, color=colors, linewidth=0)\n", + "\n", + "# Optional: overlay hydrogen consumption (as before)\n", + "n.statistics.energy_balance(aggregate_time=False).loc[\n", + " \"Link\", \"direct_reduction_furnace\", \"hydrogen\"\n", + "].loc[pd.Timestamp(\"2050-07-21\") : pd.Timestamp(\"2050-07-26\")].plot(\n", + " ax=ax, color=\"black\", linewidth=1.5, label=\"DRI H$_2$ consumption\"\n", + ")\n", + "\n", + "# Styling (matching your original)\n", + "ax.grid(True, alpha=0.5, linestyle=\"--\")\n", + "ax.set_xlabel(\"Day\") #, fontsize=9)\n", + "ax.set_ylabel(\"Electricity Balance in MW\")\n", + "\n", + "ax.legend(\n", + " loc=\"center left\",\n", + " bbox_to_anchor=(1.02, 0.5),\n", + ")\n", + "\n", + "plt.tight_layout()\n", + "plt.savefig(snakemake.output.hourly_plot, dpi=300, bbox_inches=\"tight\")\n", + "plt.show()" ] }, { "cell_type": "markdown", - "id": "3b90f701", + "id": "7896b568", "metadata": {}, "source": [ - "#### Steel balance" + "### Plots in hourly resolution" ] }, { "cell_type": "code", "execution_count": null, - "id": "68bc351e", + "id": "02b04d02", "metadata": {}, "outputs": [], "source": [ - "n.statistics.energy_balance.iplot.area(bus_carrier=\"steel\")" + "# n.statistics.energy_balance.iplot.area(bus_carrier=\"hydrogen\")" ] }, { @@ -187,104 +212,7 @@ "metadata": {}, "outputs": [], "source": [ - "n.statistics.energy_balance.iplot.area() # Ely als negativen Stromverbraucher, H2 Verbrauch der DRI als linie" - ] - }, - { - "cell_type": "markdown", - "id": "6ec60deb", - "metadata": {}, - "source": [ - "### Hourly plot (clean)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e6320390", - "metadata": {}, - "outputs": [], - "source": [ - "df = n.statistics.energy_balance(aggregate_time=False).loc[\n", - " :,\n", - " [\n", - " \"Photovoltaics\",\n", - " \"Wind energy\",\n", - " \"battery inverter (charging)\",\n", - " \"battery inverter (discharging)\",\n", - " \"electrolysis\",\n", - " \"electric arc furnace\",\n", - " \"direct reduction furnace\",\n", - " ],\n", - " \"Electricity\",\n", - "]\n", - "\n", - "df = df.loc[:, pd.Timestamp(\"2013-07-23\") : pd.Timestamp(\"2013-07-28\")]\n", - "\n", - "df = (\n", - " df.reset_index(level=[\"component\", \"bus_carrier\"], drop=True)\n", - " .groupby(\"carrier\")\n", - " .sum()\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1d7a77d0", - "metadata": {}, - "outputs": [], - "source": [ - "config_fn = \"../../config/config.yaml\"\n", - "colors = load_config(config_fn)[\"colors\"]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7f757dfe", - "metadata": {}, - "outputs": [], - "source": [ - "show_carriers = [\n", - " \"electrolysis\",\n", - " \"Photovoltaics\",\n", - " \"Wind energy\",\n", - " \"battery inverter (charging)\",\n", - " \"battery inverter (discharging)\",\n", - " \"electric arc furnace\",\n", - " \"direct reduction furnace\",\n", - "]\n", - "\n", - "# transpose so time is on x-axis and divide by 4 to have smaller system\n", - "system_scale = 4\n", - "\n", - "df_plot = df.loc[show_carriers].T.div(system_scale)\n", - "\n", - "fig, ax = plt.subplots(figsize=(9, 4))\n", - "\n", - "df_plot.plot.area(ax=ax, stacked=True, color=colors)\n", - "\n", - "# Optional: overlay hydrogen consumption (as before)\n", - "n.statistics.energy_balance(aggregate_time=False).loc[\n", - " \"Link\", \"direct reduction furnace\", \"Hydrogen\"\n", - "].loc[pd.Timestamp(\"2013-07-23\") : pd.Timestamp(\"2013-07-28\")].div(system_scale).plot(\n", - " ax=ax, color=\"black\", linewidth=1.5, label=\"DRI H$_2$ consumption\"\n", - ")\n", - "\n", - "# Styling (matching your original)\n", - "ax.grid(True, alpha=0.5, linestyle=\"--\")\n", - "ax.set_xlabel(\"Day\", fontsize=9)\n", - "ax.set_ylabel(\"Electricity Balance in MW\", fontsize=9)\n", - "ax.set_xlim(pd.Timestamp(\"2013-07-23\"), pd.Timestamp(\"2013-07-28\"))\n", - "\n", - "ax.legend(\n", - " loc=\"center left\",\n", - " bbox_to_anchor=(1.02, 0.5),\n", - ")\n", - "\n", - "plt.tight_layout()\n", - "plt.show()" + "# n.statistics.energy_balance.iplot.area() # Ely als negativen Stromverbraucher, H2 Verbrauch der DRI als linie" ] }, { @@ -302,7 +230,7 @@ "metadata": {}, "outputs": [], "source": [ - "n.statistics.capacity_factor.plot.bar()" + "# n.statistics.capacity_factor.plot.bar()" ] }, { @@ -348,40 +276,20 @@ "metadata": {}, "outputs": [], "source": [ - "# Electricity to steel ratio\n", - "steel_demand = (\n", - " n.statistics.energy_balance().loc[:, \"Steel\", \"Steel\"].values[0]\n", - ") # in t of steel\n", - "pv_supply = (\n", - " n.statistics.energy_balance().loc[:, \"Photovoltaics\", \"Electricity\"].values[0]\n", - ") # in MWh\n", - "wind_supply = (\n", - " n.statistics.energy_balance().loc[:, \"Wind energy\", \"Electricity\"].values[0]\n", - ") # in MWh\n", + "# # Electricity to steel ratio\n", + "# steel_demand = (\n", + "# n.statistics.energy_balance().loc[:, \"Steel\", \"Steel\"].values[0]\n", + "# ) # in t of steel\n", + "# pv_supply = (\n", + "# n.statistics.energy_balance().loc[:, \"Photovoltaics\", \"Electricity\"].values[0]\n", + "# ) # in MWh\n", + "# wind_supply = (\n", + "# n.statistics.energy_balance().loc[:, \"Wind energy\", \"Electricity\"].values[0]\n", + "# ) # in MWh\n", "\n", - "steel_demand\n", - "elec_to_steel = (pv_supply + wind_supply) * (-1) / steel_demand # in MWh/t\n", - "elec_to_steel" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "155c909e", - "metadata": {}, - "outputs": [], - "source": [ - "pv_supply" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "1029014d", - "metadata": {}, - "outputs": [], - "source": [ - "wind_supply" + "# steel_demand\n", + "# elec_to_steel = (pv_supply + wind_supply) * (-1) / steel_demand # in MWh/t\n", + "# elec_to_steel" ] }, { @@ -399,7 +307,7 @@ "metadata": {}, "outputs": [], "source": [ - "n.statistics.energy_balance().loc[\"Load\", \"Steel\", :].iloc[0] # tonnes of steel" + "# n.statistics.energy_balance().loc[\"Load\", \"Steel\", :].iloc[0] # tonnes of steel" ] }, { @@ -409,9 +317,9 @@ "metadata": {}, "outputs": [], "source": [ - "n.statistics.optimal_capacity().mul(1e3) / n.statistics.energy_balance().loc[\n", - " \"Load\", \"Steel\", :\n", - "].iloc[0]" + "# n.statistics.optimal_capacity().mul(1e3) / n.statistics.energy_balance().loc[\n", + "# \"Load\", \"Steel\", :\n", + "# ].iloc[0]" ] } ], diff --git a/workflow/notebooks/compare-scenarios.ipynb b/workflow/notebooks/compare-scenarios.ipynb index cbf005f..796c084 100644 --- a/workflow/notebooks/compare-scenarios.ipynb +++ b/workflow/notebooks/compare-scenarios.ipynb @@ -239,6 +239,27 @@ "nc.statistics.market_value()" ] }, + { + "cell_type": "markdown", + "id": "e212e772", + "metadata": {}, + "source": [ + "### Cost increase in self-sufficient scenario " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "da75ff94", + "metadata": {}, + "outputs": [], + "source": [ + "# Cost increase from default to hightrans in %\n", + "\n", + "increase = (plot_df.sum(axis=1).loc[\"hightrans\"] - plot_df.sum(axis=1).loc[\"default\"]) / plot_df.sum(axis=1).loc[\"default\"] * 100\n", + "print(f\"Cost increase from default to hightrans: {increase:.2f}%\")" + ] + }, { "cell_type": "markdown", "id": "9e0d9359", diff --git a/workflow/notebooks/plot-compare-lcox.ipynb b/workflow/notebooks/plot-compare-lcox.ipynb index ac6023c..0926d4e 100644 --- a/workflow/notebooks/plot-compare-lcox.ipynb +++ b/workflow/notebooks/plot-compare-lcox.ipynb @@ -19,14 +19,15 @@ "metadata": {}, "outputs": [], "source": [ - "from _helpers_notebooks import mock_snakemake\n", + "if \"snakemake\" not in globals():\n", + " from _helpers_notebooks import mock_snakemake\n", "\n", - "snakemake = mock_snakemake(\n", - " \"plot_compare_lcox\",\n", - " wacc=\"regional\", # uniform or regional\n", - " cost_year=\"2050\",\n", - " trade_chain=\"supplyconstraint\", # default, mga-stability-weighted, mga-chokepoints, mga-blocs, constrain-supply\n", - ")" + " snakemake = mock_snakemake(\n", + " \"plot_compare_lcox\",\n", + " wacc=\"regional\", # uniform or regional\n", + " cost_year=\"2050\",\n", + " trade_chain=\"supplyconstraint\", # default, mga-stability-weighted, mga-chokepoints, mga-blocs, constrain-supply\n", + " )" ] }, { @@ -392,14 +393,16 @@ "\n", "# --- Inner x-axis: quantity labels ---\n", "ax.set_xticks(x_positions)\n", - "qty_labels = [f\"{row['quantity']:.4g} Mt\" for _, row in plot_df.iterrows()]\n", + "qty_labels = [f\"{row['quantity']:.4g}\" for _, row in plot_df.iterrows()]\n", "ax.set_xticklabels(qty_labels, rotation=0, fontsize=9)\n", "ax.tick_params(axis=\"x\", length=0)\n", "\n", - "# --- Outer x-axis: region group labels with bracket ---\n", + "# --- Top labels: region group bracket with text above ---\n", "xaxis_transform = ax.get_xaxis_transform()\n", - "y_bracket = -0.11\n", - "y_label_pts = -7\n", + "\n", + "\n", + "y_bracket = 1.01\n", + "y_label_pts = 6\n", "\n", "for g in groups_ordered:\n", " center = group_centers[g]\n", @@ -421,7 +424,7 @@ " xytext=(0, y_label_pts),\n", " textcoords=\"offset points\",\n", " ha=\"center\",\n", - " va=\"top\",\n", + " va=\"bottom\",\n", " fontsize=10,\n", " annotation_clip=False,\n", " )\n", @@ -478,6 +481,8 @@ " loc=\"upper left\",\n", ")\n", "\n", + "ax.set_xlabel(\"Quantity of HBI in Mt$_{hbi}$/a\", labelpad=10)\n", + "\n", "plt.ylim(0,830)\n", "plt.tight_layout()\n", "plt.savefig(snakemake.output.lcox_comparison, dpi=300, bbox_inches=\"tight\")\n", diff --git a/workflow/notebooks/plot-magnitude-pull.ipynb b/workflow/notebooks/plot-magnitude-pull.ipynb new file mode 100644 index 0000000..67ca214 --- /dev/null +++ b/workflow/notebooks/plot-magnitude-pull.ipynb @@ -0,0 +1,602 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "79525d77", + "metadata": {}, + "outputs": [], + "source": [ + "import pypsa\n", + "import pandas as pd\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9a62fe91", + "metadata": {}, + "outputs": [], + "source": [ + "if \"snakemake\" not in globals():\n", + " from _helpers_notebooks import mock_snakemake\n", + " snakemake = mock_snakemake(\n", + " \"plot_magnitude_pull\",\n", + " scenario=\"default\", # default, penalty-sa, penalty-ea, penalty-nwa, penalty-oc\n", + " wacc=\"regional\", # regional, uniform\n", + " cost_year=\"2050\",\n", + " interone=\"hbi\",\n", + " trade_chain=\"supplyconstraint\"\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "03daa6b5", + "metadata": {}, + "outputs": [], + "source": [ + "add_iron_ore_cost = True\n", + "ironorecost = 97 * 1.59" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0e0f78f6", + "metadata": {}, + "outputs": [], + "source": [ + "config = snakemake.config" + ] + }, + { + "cell_type": "markdown", + "id": "49f352ef", + "metadata": {}, + "source": [ + "### Get steel demand" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4d634892", + "metadata": {}, + "outputs": [], + "source": [ + "steel_demand = pd.read_csv(snakemake.input.steel_demand, index_col=\"region\")" + ] + }, + { + "cell_type": "markdown", + "id": "d05e9f9a", + "metadata": {}, + "source": [ + "### Get costs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9e3a3904", + "metadata": {}, + "outputs": [], + "source": [ + "df_all = pd.DataFrame()\n", + "\n", + "# Select input files directly from snakemake (avoids hardcoded paths)\n", + "input_files = snakemake.input.supply_curves_interone\n", + "\n", + "for fn in input_files:\n", + " # Extract region from filename pattern: {region}_{process}.csv\n", + " region = fn.split(\"/\")[-1].split(\"_marginal_cost\")[0]\n", + "\n", + " df = pd.read_csv(fn)\n", + " df[\"region\"] = region\n", + "\n", + " # Calculate the difference in demand to get the quantity supplied at each cost step\n", + " diff = df[\"demand [t]\"].diff()\n", + " diff[0] = diff[1] # set first value to second to avoid NaN\n", + " df[\"demand [t]\"] = diff\n", + "\n", + " df_all = pd.concat([df_all, df], ignore_index=True)\n", + "\n", + "# Remove rows where the lcox [EUR/t] value is NaN (e.g. due to missing data for some regions)\n", + "df_all = df_all.dropna(subset=[\"lcox [EUR/t]\"])\n", + "\n", + "\n", + "if add_iron_ore_cost:\n", + " df_all[\"lcox [EUR/t]\"] += ironorecost\n", + "\n", + "df_all.set_index(\"region\", inplace=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4c67e779", + "metadata": {}, + "outputs": [], + "source": [ + "costs = df_all[[\"lcox [EUR/t]\", \"product_demand_mt\"]]" + ] + }, + { + "cell_type": "markdown", + "id": "079aa89a", + "metadata": {}, + "source": [ + "#### Get costs at demand level" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "23a3a732", + "metadata": {}, + "outputs": [], + "source": [ + "def interpolate_cost(region, demand, costs):\n", + " if region not in costs.index:\n", + " return np.nan\n", + "\n", + " sub = costs.loc[[region]].sort_values(\"product_demand_mt\")\n", + " demands = sub[\"product_demand_mt\"].values\n", + " lcox = sub[\"lcox [EUR/t]\"].values\n", + "\n", + " # exact match\n", + " if demand in demands:\n", + " return lcox[demands == demand].mean()\n", + "\n", + " # below smallest or above largest -> take the closest edge value\n", + " if demand <= demands[0]:\n", + " return lcox[0]\n", + " if demand >= demands[-1]:\n", + " return lcox[-1]\n", + "\n", + " # find bracketing points\n", + " idx_upper = np.searchsorted(demands, demand)\n", + " idx_lower = idx_upper - 1\n", + "\n", + " return (lcox[idx_lower] + lcox[idx_upper]) / 2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4c78b3db", + "metadata": {}, + "outputs": [], + "source": [ + "steel_demand[\"lcox_EUR_per_t\"] = steel_demand.apply(\n", + " lambda row: interpolate_cost(row.name, row[\"SteelDemand_DRI_Mt\"], costs),\n", + " axis=1\n", + ")\n", + "steel_cost_demand = steel_demand" + ] + }, + { + "cell_type": "markdown", + "id": "2f0e2cd5", + "metadata": {}, + "source": [ + "### Steel production capacities" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41ca4363", + "metadata": {}, + "outputs": [], + "source": [ + "steel_production = pd.read_csv(snakemake.input.steel_production)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1c930bbe", + "metadata": {}, + "outputs": [], + "source": [ + "# Build ISO3 -> region mapping from config\n", + "regions = config[\"regions\"]\n", + "iso3_to_region = {code: region for region, codes in regions.items() for code in codes}\n", + "\n", + "# Map each country row to its region\n", + "steel_production_new = steel_production.copy()\n", + "steel_production_new[\"region\"] = steel_production_new[\"iso3\"].map(\n", + " lambda x: iso3_to_region.get(x, \"Other\")\n", + ")\n", + "\n", + "# Sum crude steel production and population per region, then compute per-capita\n", + "steel_production_by_region = (\n", + " steel_production_new\n", + " .groupby(\"region\")[[\"crude_steel_production_t\", \"population\"]]\n", + " .sum()\n", + ")\n", + "steel_production_by_region[\"SteelProductionMt\"] = (\n", + " steel_production_by_region[\"crude_steel_production_t\"] / 1e6\n", + ")\n", + "steel_production_by_region[\"SteelProductionMtPerCapita\"] = (\n", + " steel_production_by_region[\"crude_steel_production_t\"]\n", + " / steel_production_by_region[\"population\"]\n", + ")\n", + "steel_production_by_region = steel_production_by_region.sort_values(\n", + " \"SteelProductionMt\", ascending=False\n", + ")\n", + "\n", + "steel_production_by_region[[\"SteelProductionMt\", \"SteelProductionMtPerCapita\"]]\n" + ] + }, + { + "cell_type": "markdown", + "id": "81da539b", + "metadata": {}, + "source": [ + "## Plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "48c7dcad", + "metadata": {}, + "outputs": [], + "source": [ + "import seaborn as sns\n", + "from matplotlib.lines import Line2D\n", + "\n", + "# --- user-facing plotting controls ---\n", + "plot_style = {\n", + " \"neutral_color\": \"#797777\",\n", + " \"axis_color\": \"#000000\",\n", + " \"show_axes_box\": True, # toggle outer box/frame\n", + "}\n", + "\n", + "# Base label position: always above bubble. Per-region x/y are additional offsets.\n", + "base_label_offset = {\"x\": 0, \"y\": 7}\n", + "\n", + "# --- merge production capacity + cost/demand on region ---\n", + "plot_df = steel_cost_demand.join(\n", + " steel_production_by_region[[\"SteelProductionMtPerCapita\"]], how=\"inner\"\n", + ")\n", + "plot_df = plot_df.reset_index() # region: index -> column\n", + "\n", + "nice_names = config.get(\"region_nice_names\", {})\n", + "plot_df[\"label\"] = plot_df[\"region\"].map(lambda r: nice_names.get(r, r))\n", + "\n", + "# --- per-region manual label controls (hard-coded block) ---\n", + "# x, y: additional offset points relative to base top position; line: draw connector line\n", + "label_layout = {\n", + " \"Central_America\": {\"x\": 30, \"y\": 3, \"line\": False},\n", + " \"East_Asia\": {\"x\": -28, \"y\": 10, \"line\": False},\n", + " \"East_East_Asia\": {\"x\": 0, \"y\": 0, \"line\": False},\n", + " \"Eurasia\": {\"x\": 20, \"y\": 0, \"line\": False},\n", + " \"Europe\": {\"x\": 0, \"y\": 5, \"line\": False},\n", + " \"Far_West_Europe\": {\"x\": 30, \"y\": -40, \"line\": True},\n", + " \"Middle_East\": {\"x\": 0, \"y\": 0, \"line\": False},\n", + " \"North_America\": {\"x\": 0, \"y\": 0, \"line\": False},\n", + " \"North_West_Africa\": {\"x\": 16, \"y\": -30, \"line\": True},\n", + " \"Oceania\": {\"x\": 0, \"y\": 0, \"line\": False},\n", + " \"Other\": {\"x\": 0, \"y\": 0, \"line\": False},\n", + " \"Pacific_Asia\": {\"x\": -20, \"y\": 20, \"line\": True},\n", + " \"South_America\": {\"x\": 0, \"y\": 5, \"line\": False},\n", + " \"South_South_America\": {\"x\": -10, \"y\": -30, \"line\": True},\n", + " \"Subsaharan_Africa\": {\"x\": 15, \"y\": 5, \"line\": False},\n", + " \"West_Asia\": {\"x\": 35, \"y\": 15, \"line\": False},\n", + "}\n", + "\n", + "# --- quadrant thresholds (medians) ---\n", + "x_mid = plot_df[\"SteelProductionMtPerCapita\"].median()\n", + "y_mid = plot_df[\"lcox_EUR_per_t\"].median()\n", + "\n", + "# --- bubble sizes: area proportional to DRI demand ---\n", + "max_bubble = 2500\n", + "demand_max = plot_df[\"SteelDemand_DRI_Mt\"].max()\n", + "plot_df[\"bubble\"] = (plot_df[\"SteelDemand_DRI_Mt\"] / demand_max * max_bubble).clip(lower=30)\n", + "\n", + "# --- region colours from config.yaml ---\n", + "region_color_map = {\n", + " r: config.get(\"colors\", {}).get(r, \"#666666\")\n", + " for r in plot_df[\"region\"].tolist()\n", + "}\n", + "\n", + "fig, ax = plt.subplots(figsize=(8.8, 6.2))\n", + "\n", + "plot_df_ordered = plot_df.assign(_plot_first=(plot_df[\"region\"] == \"West_Asia\").astype(int))\n", + "plot_df_ordered = plot_df_ordered.sort_values(\"_plot_first\", ascending=False)\n", + "\n", + "for i, (_, row) in enumerate(plot_df_ordered.iterrows()):\n", + " x_val = row[\"SteelProductionMtPerCapita\"]\n", + " y_val = row[\"lcox_EUR_per_t\"]\n", + " bubble_area = row[\"bubble\"]\n", + " label_text = row[\"label\"]\n", + " region_key = row[\"region\"]\n", + " point_color = region_color_map.get(region_key, \"#666666\")\n", + "\n", + " ax.scatter(\n", + " x_val,\n", + " y_val,\n", + " s=bubble_area,\n", + " color=point_color,\n", + " alpha=0.8,\n", + " edgecolors=\"white\",\n", + " linewidths=0.8,\n", + " zorder=3,\n", + " )\n", + "\n", + " region_cfg = label_layout.get(region_key, {\"x\": 0, \"y\": 0, \"line\": False})\n", + " label_offset = (\n", + " base_label_offset[\"x\"] + region_cfg.get(\"x\", 0),\n", + " base_label_offset[\"y\"] + region_cfg.get(\"y\", 0),\n", + " )\n", + " draw_line = bool(region_cfg.get(\"line\", False))\n", + "\n", + " annotate_kwargs = {\n", + " \"xy\": (x_val, y_val),\n", + " \"xytext\": label_offset,\n", + " \"textcoords\": \"offset points\",\n", + " \"fontsize\": 8.5,\n", + " \"color\": point_color,\n", + " \"ha\": \"center\",\n", + " \"va\": \"bottom\",\n", + " \"zorder\": 4,\n", + " }\n", + " if draw_line:\n", + " annotate_kwargs[\"arrowprops\"] = {\n", + " \"arrowstyle\": \"-\",\n", + " \"lw\": 0.8,\n", + " \"color\": plot_style[\"neutral_color\"],\n", + " \"alpha\": 0.9,\n", + " \"shrinkA\": 2,\n", + " \"shrinkB\": 7,\n", + " }\n", + "\n", + " ax.annotate(label_text, **annotate_kwargs)\n", + "\n", + "# --- quadrant dividers ---\n", + "# ax.axvline(x_mid, color=plot_style[\"neutral_color\"], linestyle=\"--\", linewidth=1.0, alpha=0.6, zorder=2)\n", + "# ax.axhline(y_mid, color=plot_style[\"neutral_color\"], linestyle=\"--\", linewidth=1.0, alpha=0.6, zorder=2)\n", + "\n", + "ax.set_ylim(600, 850)\n", + "ax.set_xlim(-0.09, 0.8)\n", + "\n", + "# --- quadrant text annotations ---\n", + "xl, xr = ax.get_xlim()\n", + "yb, yt = ax.get_ylim()\n", + "px = (xr - xl) * 0.02\n", + "py = (yt - yb) * 0.02\n", + "\n", + "quadrant_labels = [\n", + " (xl + px, yt - py-20, \"top\", \"left\", \"Demand-anchored\\nemerging production centres\"),\n", + " (xr - px, yt - py, \"top\", \"right\", \"Exposed incumbent centres\"),\n", + " (xl + px, yb + py, \"bottom\", \"left\", \"Potential emerging\\ngreen-iron producers\"),\n", + " (xr - px, yb + py, \"bottom\", \"right\", \"Well-positioned\\nincumbent centres\"),\n", + "]\n", + "for xq, yq, va, ha, txt in quadrant_labels:\n", + " ax.text(\n", + " xq,\n", + " yq,\n", + " txt,\n", + " fontsize=8,\n", + " va=va,\n", + " ha=ha,\n", + " color=plot_style[\"neutral_color\"],\n", + " style=\"italic\",\n", + " alpha=0.9,\n", + " zorder=2,\n", + " )\n", + "\n", + "# --- inner-axis directional annotations ---\n", + "# Manual fine-tuning positions in axes coordinates.\n", + "directional_guides = {\n", + " \"renewables\": {\n", + " \"label_x\": 0.006, \"label_y\": 0.5,\n", + " \"high_x\": 0.008, \"high_y\": 0.71,\n", + " \"low_x\": 0.008, \"low_y\": 0.29,\n", + " \"arrow_x\": 0.015,\n", + " \"arrow_high_start\": 0.74, \"arrow_high_end\": 0.80,\n", + " \"arrow_low_start\": 0.25, \"arrow_low_end\": 0.19,\n", + " },\n", + " \"industry\": {\n", + " \"label_x\": 0.5, \"label_y\": 0.005,\n", + " \"low_x\": 0.38, \"low_y\": 0.01,\n", + " \"high_x\": 0.62, \"high_y\": 0.01,\n", + " \"arrow_y\": 0.02,\n", + " \"arrow_low_start\": 0.34, \"arrow_low_end\": 0.28,\n", + " \"arrow_high_start\": 0.66, \"arrow_high_end\": 0.72,\n", + " },\n", + "}\n", + "\n", + "# RENEWABLES PULL (vertical): LOW arrow down, HIGH arrow up\n", + "ax.text(\n", + " directional_guides[\"renewables\"][\"label_x\"],\n", + " directional_guides[\"renewables\"][\"label_y\"],\n", + " \"RENEWABLES PULL\",\n", + " transform=ax.transAxes,\n", + " rotation=90,\n", + " fontsize=14,\n", + " color=plot_style[\"neutral_color\"],\n", + " ha=\"left\",\n", + " va=\"center\",\n", + " alpha=0.7,\n", + ")\n", + "ax.text(\n", + " directional_guides[\"renewables\"][\"high_x\"],\n", + " directional_guides[\"renewables\"][\"high_y\"],\n", + " \"HIGH\",\n", + " transform=ax.transAxes,\n", + " rotation=90,\n", + " fontsize=7,\n", + " color=plot_style[\"neutral_color\"],\n", + " ha=\"left\",\n", + " va=\"center\",\n", + " alpha=0.7,\n", + ")\n", + "ax.text(\n", + " directional_guides[\"renewables\"][\"low_x\"],\n", + " directional_guides[\"renewables\"][\"low_y\"],\n", + " \"LOW\",\n", + " transform=ax.transAxes,\n", + " rotation=90,\n", + " fontsize=7,\n", + " color=plot_style[\"neutral_color\"],\n", + " ha=\"left\",\n", + " va=\"center\",\n", + " alpha=0.7,\n", + ")\n", + "ax.annotate(\n", + " \"\",\n", + " xy=(directional_guides[\"renewables\"][\"arrow_x\"], directional_guides[\"renewables\"][\"arrow_high_end\"]),\n", + " xytext=(directional_guides[\"renewables\"][\"arrow_x\"], directional_guides[\"renewables\"][\"arrow_high_start\"]),\n", + " xycoords=ax.transAxes,\n", + " arrowprops=dict(arrowstyle=\"->\", lw=0.8, color=plot_style[\"neutral_color\"], alpha=0.7),\n", + " \n", + ")\n", + "ax.annotate(\n", + " \"\",\n", + " xy=(directional_guides[\"renewables\"][\"arrow_x\"], directional_guides[\"renewables\"][\"arrow_low_end\"]),\n", + " xytext=(directional_guides[\"renewables\"][\"arrow_x\"], directional_guides[\"renewables\"][\"arrow_low_start\"]),\n", + " xycoords=ax.transAxes,\n", + " arrowprops=dict(arrowstyle=\"->\", lw=0.8, color=plot_style[\"neutral_color\"]),\n", + ")\n", + "\n", + "# INDUSTRY (horizontal): LOW arrow left, HIGH arrow right\n", + "ax.text(\n", + " directional_guides[\"industry\"][\"label_x\"],\n", + " directional_guides[\"industry\"][\"label_y\"],\n", + " \"INDUSTRY TODAY\",\n", + " transform=ax.transAxes,\n", + " fontsize=14,\n", + " color=plot_style[\"neutral_color\"],\n", + " ha=\"center\",\n", + " va=\"bottom\",\n", + " alpha=0.7,\n", + ")\n", + "ax.text(\n", + " directional_guides[\"industry\"][\"low_x\"],\n", + " directional_guides[\"industry\"][\"low_y\"],\n", + " \"LOW\",\n", + " transform=ax.transAxes,\n", + " fontsize=7,\n", + " color=plot_style[\"neutral_color\"],\n", + " ha=\"right\",\n", + " va=\"bottom\",\n", + " alpha=0.7,\n", + ")\n", + "ax.annotate(\n", + " \"\",\n", + " xy=(directional_guides[\"industry\"][\"arrow_low_end\"], directional_guides[\"industry\"][\"arrow_y\"]),\n", + " xytext=(directional_guides[\"industry\"][\"arrow_low_start\"], directional_guides[\"industry\"][\"arrow_y\"]),\n", + " xycoords=ax.transAxes,\n", + " arrowprops=dict(arrowstyle=\"->\", lw=0.8, color=plot_style[\"neutral_color\"],alpha=0.7,),\n", + ")\n", + "ax.text(\n", + " directional_guides[\"industry\"][\"high_x\"],\n", + " directional_guides[\"industry\"][\"high_y\"],\n", + " \"HIGH\",\n", + " transform=ax.transAxes,\n", + " fontsize=7,\n", + " color=plot_style[\"neutral_color\"],\n", + " ha=\"left\",\n", + " va=\"bottom\",\n", + " alpha=0.7,\n", + ")\n", + "ax.annotate(\n", + " \"\",\n", + " xy=(directional_guides[\"industry\"][\"arrow_high_end\"], directional_guides[\"industry\"][\"arrow_y\"]),\n", + " xytext=(directional_guides[\"industry\"][\"arrow_high_start\"], directional_guides[\"industry\"][\"arrow_y\"]),\n", + " xycoords=ax.transAxes,\n", + " arrowprops=dict(arrowstyle=\"->\", lw=0.8, color=plot_style[\"neutral_color\"],alpha=0.7,),\n", + ")\n", + "\n", + "# --- axis labels (no title) ---\n", + "ax.set_xlabel(\n", + " \"Today's Steel Production Capacity in t$_{steel}$ per capita\",\n", + " fontsize=12,\n", + " color=plot_style[\"axis_color\"],\n", + ")\n", + "ax.set_ylabel(\n", + " \"Domestic HBI costs in 2050 in €/t$_{hbi}$ \\nas a proxy for the magnitude of Renewables Pull\",\n", + " fontsize=12,\n", + " color=plot_style[\"axis_color\"],\n", + ")\n", + "\n", + "# --- use neutral colour for all non-bubble elements ---\n", + "ax.tick_params(axis=\"both\", colors=plot_style[\"axis_color\"])\n", + "for spine in ax.spines.values():\n", + " spine.set_visible(plot_style[\"show_axes_box\"])\n", + " spine.set_color(plot_style[\"axis_color\"])\n", + "\n", + "# --- bubble-size legend outside right; spaced markers; no frame ---\n", + "legend_values = [10, 50, 100]\n", + "legend_handles = []\n", + "for mt_val in legend_values:\n", + " s = (mt_val / demand_max) * max_bubble\n", + " marker_size = np.sqrt(s) * 0.72 # convert area-like scale to marker size\n", + " legend_handles.append(\n", + " Line2D(\n", + " [0],\n", + " [0],\n", + " marker=\"o\",\n", + " linestyle=\"\",\n", + " markerfacecolor=plot_style[\"neutral_color\"],\n", + " markeredgecolor=plot_style[\"neutral_color\"],\n", + " alpha=0.55,\n", + " markersize=marker_size,\n", + " label=f\"{mt_val} Mt\",\n", + " )\n", + " )\n", + "\n", + "legend = ax.legend(\n", + " handles=legend_handles,\n", + " title=\"HBI Demand in 2050\",\n", + " loc=\"center left\",\n", + " bbox_to_anchor=(0.75, 0.5),\n", + " frameon=True,\n", + " fontsize=9,\n", + " title_fontsize=9,\n", + " labelspacing=1.4,\n", + " handletextpad=1.2,\n", + " borderaxespad=1.2,\n", + " borderpad=1.0, \n", + ")\n", + "legend.get_title().set_color(plot_style[\"axis_color\"])\n", + "for t in legend.get_texts():\n", + " t.set_color(plot_style[\"neutral_color\"])\n", + "\n", + "plt.tight_layout()\n", + "plt.savefig(snakemake.output.magnitude_pull_png, dpi=150, bbox_inches=\"tight\")\n", + "plt.savefig(snakemake.output.magnitude_pull, bbox_inches=\"tight\")\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "shift", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/workflow/notebooks/plot-mga.ipynb b/workflow/notebooks/plot-mga.ipynb index 7f122b7..1422ade 100644 --- a/workflow/notebooks/plot-mga.ipynb +++ b/workflow/notebooks/plot-mga.ipynb @@ -1481,7 +1481,9 @@ "# Combined plot — v2\n", "# ═══════════════════════════════════════════════════════════════════════════════\n", "Y_MAX_V2 = 100\n", - "GAP_FRAC_V2 = 0.04\n", + "GAP_FRAC_V2 = 0.0\n", + "ARROW_LEN_V2 = 14 # fixed arrow length in y-axis units for consistent visuals\n", + "ARROW_Y_START_V2 = 10 # default arrow start y-position (manual control)\n", "\n", "plots_v2 = [\n", " dict(\n", @@ -1492,7 +1494,9 @@ " y_mid=cp_total_vals_v2_pct[0],\n", " color=\"#000000\",\n", " arrow_dn=[\"Minimize\", \"high-risk routes\"],\n", - " arrow_x_frac=0.40,\n", + " arrow_x_frac=0.20,\n", + " arrow_y_start=20,\n", + " arrow_label_y=15,\n", " y_label_add=\"trade through maritime chokepoints \\nin % of market share\",\n", " regions=cp_region_lines_v2_pct,\n", " extra_curves=[],\n", @@ -1506,6 +1510,8 @@ " color=\"#000000\",\n", " arrow_dn=[\"Minimize trade with\", \"non-friendshored regions\"],\n", " arrow_x_frac=0.40,\n", + " arrow_y_start=25,\n", + " arrow_label_y=20,\n", " y_label_add=\"trade between groups \\nin % of market share\",\n", " regions=bl_region_lines_v2_pct,\n", " extra_curves=[\n", @@ -1521,7 +1527,9 @@ " y_mid=sc_cap_vals_pct[0],\n", " color=\"#000000\",\n", " arrow_dn=[\"Constrain production\", \"per supplier\"],\n", - " arrow_x_frac=0.30,\n", + " arrow_x_frac=0.20,\n", + " arrow_y_start=40,\n", + " arrow_label_y=35,\n", " y_label_add=\"production by region \\nin % of market share\",\n", " regions=sc_region_lines_pct,\n", " extra_curves=[],\n", @@ -1536,8 +1544,7 @@ "for ax, p in zip(axes_v2, plots_v2):\n", " eps = p[\"epsilon\"]\n", " curve = p[\"main_curve\"]\n", - " y_mid = p[\"y_mid\"]\n", - " gap = Y_MAX_V2 * GAP_FRAC_V2\n", + "\n", "\n", " # ── Main curve ───────────────────────────────────────────────────────────\n", " ax.plot(eps, curve, color=p[\"color\"], linewidth=1.5, marker=\"o\",markersize=3, label=p[\"main_label\"], zorder=0, linestyle=\"-\")\n", @@ -1581,14 +1588,14 @@ " arrow_x_frac = p.get(\"arrow_x_frac\", 0.70)\n", " arrow_x = eps[0] + arrow_x_frac * (eps[-1] - eps[0])\n", " label_x = eps[0] + (arrow_x_frac + 0.05) * (eps[-1] - eps[0])\n", - " curve_at_arrow = np.interp(arrow_x, eps, curve)\n", - " y_start = y_mid - gap\n", - " y_end = curve_at_arrow\n", + " y_start = p.get(\"arrow_y_start\", ARROW_Y_START_V2)\n", + " y_end = max(0, y_start - ARROW_LEN_V2)\n", " if y_end < y_start:\n", " ax.annotate(\"\", xy=(arrow_x, y_end), xytext=(arrow_x, y_start),\n", " xycoords=\"data\", textcoords=\"data\",\n", " arrowprops=arrowprops_v2, clip_on=False)\n", - " ax.text(label_x, (y_start + y_end) / 2, \"\\n\".join(p[\"arrow_dn\"]),\n", + " label_y = p.get(\"arrow_label_y\", (y_start + y_end) / 2)\n", + " ax.text(label_x, label_y, \"\\n\".join(p[\"arrow_dn\"]),\n", " fontsize=8, va=\"center\", ha=\"left\", color=\"#444441\", clip_on=False)\n", "\n", "plt.tight_layout()\n", diff --git a/workflow/notebooks/plot-robust-map.ipynb b/workflow/notebooks/plot-robust-map.ipynb new file mode 100644 index 0000000..7545d52 --- /dev/null +++ b/workflow/notebooks/plot-robust-map.ipynb @@ -0,0 +1,222 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "070b1b23", + "metadata": {}, + "source": [ + "# Robust HBI Trade Map\n", + "Build a robust HBI map from three scenarios by taking minimum HBI trade capacity and minimum realized HBI production supply dispatch, while keeping demand unchanged." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "273e7cb1", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import re\n", + "import sys\n", + "from pathlib import Path\n", + "\n", + "import pandas as pd\n", + "import pypsa\n", + "\n", + "if \"snakemake\" not in globals():\n", + " from _helpers_notebooks import mock_snakemake\n", + " snakemake = mock_snakemake(\"plot_robust_map\")\n", + "\n", + "scripts_dir = Path.cwd() / \"workflow\" / \"scripts\"\n", + "if str(scripts_dir) not in sys.path:\n", + " sys.path.insert(0, str(scripts_dir))\n", + "\n", + "import model_trade as mt\n", + "\n", + "mt.snakemake = snakemake\n", + "trade_chain = snakemake.config[\"trade_chains\"]\n", + "mt.final = trade_chain[\"final_product\"]\n", + "tradeable = trade_chain[\"tradeable_commodities\"]\n", + "mt.interone = next(c for c in tradeable if c != \"iron_ore\")\n", + "stages = trade_chain[\"stages\"]\n", + "mt.intertwo = stages[max(stages.keys())][\"process_label\"]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "796dbd22", + "metadata": {}, + "outputs": [], + "source": [ + "scenario_input_paths = {\n", + " \"chokepoints\": snakemake.input.chokepoints,\n", + " \"blocs\": snakemake.input.blocs,\n", + " \"constrain_supply\": snakemake.input.constrain_supply,\n", + "}\n", + "\n", + "\n", + "scenario_network_paths = {\n", + " name: path for name, path in scenario_input_paths.items()\n", + "}\n", + "networks = {k: pypsa.Network(path) for k, path in scenario_network_paths.items()}\n", + "\n", + "for name, path in scenario_network_paths.items():\n", + " print(f\"Loaded {name}: {path}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e9f4593e", + "metadata": {}, + "outputs": [], + "source": [ + "def robust_min_series(series_list):\n", + " return pd.concat(series_list, axis=1).fillna(0.0).min(axis=1)\n", + "\n", + "def robust_min_dispatch(frames):\n", + " merged = pd.concat(frames, axis=1, keys=range(len(frames))).fillna(0.0)\n", + " return merged.T.groupby(level=1).min().T\n", + "\n", + "trade_carrier = \"shipping_hbi\"\n", + "production_carrier = \"hbi\"\n", + "\n", + "robust_network = networks[\"chokepoints\"].copy()\n", + "\n", + "trade_series = [\n", + " net.links.loc[net.links.carrier == trade_carrier, \"p_nom_opt\"]\n", + " for net in networks.values()\n", + "]\n", + "robust_trade = robust_min_series(trade_series)\n", + "robust_network.links.loc[robust_trade.index, \"p_nom_opt\"] = robust_trade.values\n", + "\n", + "hbi_link_names = robust_network.links.index[robust_network.links.carrier == production_carrier]\n", + "has_dispatch = all(\n", + " hasattr(net, \"links_t\") and (\"p1\" in net.links_t) and (not net.links_t[\"p1\"].empty)\n", + " for net in networks.values()\n", + ")\n", + "\n", + "if has_dispatch:\n", + " production_dispatch_frames = [\n", + " (-net.links_t[\"p1\"].loc[:, net.links.index[net.links.carrier == production_carrier]]).clip(lower=0)\n", + " for net in networks.values()\n", + " ]\n", + " robust_dispatch = robust_min_dispatch(production_dispatch_frames)\n", + "\n", + " robust_network.links_t[\"p1\"].loc[:, robust_dispatch.columns] = -robust_dispatch.values\n", + "\n", + " if \"p0\" in robust_network.links_t:\n", + " eff = robust_network.links.loc[robust_dispatch.columns, \"efficiency\"].replace(0, pd.NA)\n", + " robust_network.links_t[\"p0\"].loc[:, robust_dispatch.columns] = (\n", + " robust_dispatch.div(eff, axis=1).fillna(0.0).values\n", + " )\n", + "\n", + " robust_network.links.loc[robust_dispatch.columns, \"p_nom_opt\"] = robust_dispatch.max(axis=0).values\n", + " robust_production = robust_dispatch.sum(axis=0)\n", + "else:\n", + " production_series = [\n", + " net.links.loc[net.links.carrier == production_carrier, \"p_nom_opt\"]\n", + " for net in networks.values()\n", + " ]\n", + " robust_production = robust_min_series(production_series)\n", + " robust_network.links.loc[robust_production.index, \"p_nom_opt\"] = robust_production.values\n", + "\n", + "print(f\"Robust trade links: {len(robust_trade)}\")\n", + "print(f\"Robust production links: {len(hbi_link_names)}\")\n", + "print(f\"Total robust HBI trade [Mt]: {robust_trade.sum() / 1e6:.2f}\")\n", + "print(f\"Total robust HBI production [Mt]: {robust_production.sum() / 1e6:.2f}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cdfa54e2", + "metadata": {}, + "outputs": [], + "source": [ + "regions_to_check = [\"East_Asia\", \"Oceania\"]\n", + "\n", + "def region_hbi_capacity_mt(n, region):\n", + " return (\n", + " n.links.loc[\n", + " (n.links.carrier == \"hbi\") & (n.links.bus0 == f\"{region}_ore\"),\n", + " \"p_nom_opt\",\n", + " ].sum()\n", + " / 1e6\n", + " )\n", + "\n", + "def region_hbi_supply_mt(n, region):\n", + " s = (\n", + " n.statistics.supply(comps=[\"Link\"], groupby=[\"bus\", \"carrier\"])\n", + " .loc[:, :, \"hbi\"]\n", + " .droplevel(0)\n", + " )\n", + " return float(s.get(f\"{region}_hbi\", 0.0)) / 1e6\n", + "\n", + "rows = []\n", + "for region in regions_to_check:\n", + " cap_vals = {name: region_hbi_capacity_mt(net, region) for name, net in networks.items()}\n", + " flow_vals = {name: region_hbi_supply_mt(net, region) for name, net in networks.items()}\n", + "\n", + " rows.append(\n", + " {\n", + " \"region\": region,\n", + " \"min_capacity_input_scenarios\": min(cap_vals.values()),\n", + " \"robust_capacity\": region_hbi_capacity_mt(robust_network, region),\n", + " \"min_supply_input_scenarios\": min(flow_vals.values()),\n", + " \"robust_supply\": region_hbi_supply_mt(robust_network, region),\n", + " \"chokepoints_supply\": region_hbi_supply_mt(networks[\"chokepoints\"], region),\n", + " }\n", + " )\n", + "\n", + "pd.DataFrame(rows)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "49a651fb", + "metadata": {}, + "outputs": [], + "source": [ + "region_gdf = mt.build_region_geodataframe(snakemake.config)\n", + "\n", + "mt.plot_trade_network(\n", + " robust_network,\n", + " product=\"hbi\",\n", + " alpha_supply=0.5,\n", + " output_path=snakemake.output.robust_map_pdf,\n", + " output_path_png=snakemake.output.robust_map_png,\n", + " region_gdf=region_gdf,\n", + ")\n", + "\n", + "print(\"Saved robust HBI map:\")\n", + "print(f\" PDF: {snakemake.output.robust_map_pdf}\")\n", + "print(f\" PNG: {snakemake.output.robust_map_png}\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "shift", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.12.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/workflow/notebooks/plot_todays-trade.py b/workflow/notebooks/plot_todays-trade.py index 14e43c6..844c875 100644 --- a/workflow/notebooks/plot_todays-trade.py +++ b/workflow/notebooks/plot_todays-trade.py @@ -1,11 +1,11 @@ """ Today's steel trade: -Data: https://www.cepii.fr/CEPII/en/bdd_modele/bdd_modele_item.asp?id=37 --> https://www.cepii.fr/DATA_DOWNLOAD/baci/data/BACI_HS22_V202601.zip (2022-2024) - Latest update from BACI. +Data: https://www.cepii.fr/CEPII/en/bdd_modele/bdd_modele_item.asp?id=37 +-> https://www.cepii.fr/DATA_DOWNLOAD/baci/data/BACI_HS22_V202601.zip (2022-2024) - Latest update from BACI. -> workflow/notebooks/BACI_HS22_V202601 -Process (in general): +Process (in general): Iron ore -> DRI/HBI -> steel (raw) Relevant HS product codes for iron ore trade: @@ -26,7 +26,6 @@ - Source: product_codes_HS22_V202601.csv """ - import os import sys import warnings @@ -48,57 +47,165 @@ # ============================================================================== # Resolved at runtime from snakemake (see __main__ block below) PRODUCT_OPTIONS = { - "Iron Ore": [260111, 260112], - "DRI-HBI": [720310], # 720390 excluded: covers electrolytic/carbonyl iron (≥99.94% Fe), not commercial DRI/HBI + "Iron Ore": [260111, 260112], + "DRI-HBI": [ + 720310 + ], # 720390 excluded: covers electrolytic/carbonyl iron (≥99.94% Fe), not commercial DRI/HBI # Semi-finished (7207) + flat-rolled (7208-7212) + long products (7213-7217) # + stainless (7218-7223) + other alloy steel (7224-7229). # Excludes pig iron/ferro-alloys/scrap (7201-7205) and DRI/HBI (7203, already above). - "Steel (raw)": [ + "Steel": [ # Semi-finished billets / slabs - 720711, 720712, 720719, 720720, + 720711, + 720712, + 720719, + 720720, # Hot-rolled flat, wide strip - 720810, 720825, 720826, 720827, 720836, 720837, 720838, 720839, - 720840, 720851, 720852, 720853, 720854, 720890, + 720810, + 720825, + 720826, + 720827, + 720836, + 720837, + 720838, + 720839, + 720840, + 720851, + 720852, + 720853, + 720854, + 720890, # Cold-rolled flat, wide strip - 720915, 720916, 720917, 720918, 720925, 720926, 720927, 720928, 720990, + 720915, + 720916, + 720917, + 720918, + 720925, + 720926, + 720927, + 720928, + 720990, # Coated / painted flat, wide strip - 721011, 721012, 721020, 721030, 721041, 721049, 721050, - 721061, 721069, 721070, 721090, + 721011, + 721012, + 721020, + 721030, + 721041, + 721049, + 721050, + 721061, + 721069, + 721070, + 721090, # Flat-rolled, narrow strip (< 600 mm) - 721113, 721114, 721119, 721123, 721129, 721190, - 721210, 721220, 721230, 721240, + 721113, + 721114, + 721119, + 721123, + 721129, + 721190, + 721210, + 721220, + 721230, + 721240, # Wire rod - 721310, 721320, 721391, 721399, + 721310, + 721320, + 721391, + 721399, # Bars and rods - 721410, 721420, 721430, 721491, 721499, - 721510, 721550, 721590, + 721410, + 721420, + 721430, + 721491, + 721499, + 721510, + 721550, + 721590, # Sections (angles, beams, channels) - 721610, 721621, 721622, 721631, 721632, 721633, - 721640, 721650, 721661, 721669, 721691, 721699, + 721610, + 721621, + 721622, + 721631, + 721632, + 721633, + 721640, + 721650, + 721661, + 721669, + 721691, + 721699, # Wire - 721710, 721720, 721730, 721790, + 721710, + 721720, + 721730, + 721790, # Stainless steel (flat, long, wire) - 721810, 721891, 721899, - 721911, 721912, 721913, 721914, 721921, 721922, 721923, 721924, - 721931, 721932, 721933, 721934, 721935, - 722011, 722012, 722020, 722090, - 722110, 722190, 722210, 722220, 722230, 722240, 722290, + 721810, + 721891, + 721899, + 721911, + 721912, + 721913, + 721914, + 721921, + 721922, + 721923, + 721924, + 721931, + 721932, + 721933, + 721934, + 721935, + 722011, + 722012, + 722020, + 722090, + 722110, + 722190, + 722210, + 722220, + 722230, + 722240, + 722290, 722300, # Other alloy steel (flat, long, wire) - 722410, 722490, - 722511, 722519, 722520, 722530, 722540, 722550, - 722611, 722619, 722620, 722691, 722699, - 722710, 722720, 722790, - 722810, 722820, 722830, 722840, 722850, 722860, 722870, 722880, 722890, - 722910, 722920, 722990, + 722410, + 722490, + 722511, + 722519, + 722520, + 722530, + 722540, + 722550, + 722611, + 722619, + 722620, + 722691, + 722699, + 722710, + 722720, + 722790, + 722810, + 722820, + 722830, + 722840, + 722850, + 722860, + 722870, + 722880, + 722890, + 722910, + 722920, + 722990, ], } # Safe filename stems matching Snakemake output keys PRODUCT_SAFE_NAMES = { - "Iron Ore": "Iron_Ore", - "DRI-HBI": "DRI-HBI", - "Steel (raw)": "Steel_raw", + "Iron Ore": "Iron_Ore", + "DRI-HBI": "DRI-HBI", + "Steel": "Steel", } # How many bilateral net-flow arrows to draw on the map. @@ -111,7 +218,6 @@ VALUE_COLUMN = "quantity" - # ============================================================================== # 1. PREPARE TRADE DATA # ============================================================================== @@ -128,8 +234,7 @@ def prepare_trade(baci_folder, product_name, product_codes, out_csv, year): print("Aggregating bilateral trade") trade_pair = ( - df - .groupby(["i", "j"], as_index=False) + df.groupby(["i", "j"], as_index=False) .agg(trade_value=("v", "sum"), quantity=("q", "sum")) .rename(columns={"i": "exporter", "j": "importer"}) ) @@ -143,12 +248,14 @@ def prepare_trade(baci_folder, product_name, product_codes, out_csv, year): cc[["country_code", "country_iso3"]].rename( columns={"country_code": "exporter", "country_iso3": "exporter_iso3"} ), - on="exporter", how="left" + on="exporter", + how="left", ).merge( cc[["country_code", "country_iso3"]].rename( columns={"country_code": "importer", "country_iso3": "importer_iso3"} ), - on="importer", how="left" + on="importer", + how="left", ) out_iso3 = out_csv @@ -167,36 +274,41 @@ def get_net_pairwise_flows(df_trade, exp_col, imp_col, val_col): pairs.columns = ["node1", "node2", "val"] pairs["pair_key"] = pairs.apply( - lambda x: "-".join(sorted([str(x["node1"]), str(x["node2"])])), - axis=1 + lambda x: "-".join(sorted([str(x["node1"]), str(x["node2"])])), axis=1 ) net_flows = [] for _, group in pairs.groupby("pair_key"): if len(group) == 1: row = group.iloc[0] - net_flows.append({ - "exporter": row["node1"], - "importer": row["node2"], - "quantity": row["val"] - }) + net_flows.append( + { + "exporter": row["node1"], + "importer": row["node2"], + "quantity": row["val"], + } + ) else: row_a = group.iloc[0] row_b = group.iloc[1] diff = row_a["val"] - row_b["val"] if diff > 0: - net_flows.append({ - "exporter": row_a["node1"], - "importer": row_a["node2"], - "quantity": diff - }) + net_flows.append( + { + "exporter": row_a["node1"], + "importer": row_a["node2"], + "quantity": diff, + } + ) elif diff < 0: - net_flows.append({ - "exporter": row_b["node1"], - "importer": row_b["node2"], - "quantity": abs(diff) - }) + net_flows.append( + { + "exporter": row_b["node1"], + "importer": row_b["node2"], + "quantity": abs(diff), + } + ) return pd.DataFrame(net_flows) @@ -227,12 +339,16 @@ def process_data_for_plot(trade_file, world_gdf, value_column): ) net_pairwise = get_net_pairwise_flows(df_trade, exp_col, imp_col, val_col) - top_net_flows = net_pairwise.sort_values("quantity", ascending=False).head(TOP_N_FLOWS) + top_net_flows = net_pairwise.sort_values("quantity", ascending=False).head( + TOP_N_FLOWS + ) temp_world = world_gdf.copy() temp_world = temp_world.merge(total_exp, on="iso3", how="left") temp_world = temp_world.merge(total_imp, on="iso3", how="left") - temp_world[["total_exp", "total_imp"]] = temp_world[["total_exp", "total_imp"]].fillna(0) + temp_world[["total_exp", "total_imp"]] = temp_world[ + ["total_exp", "total_imp"] + ].fillna(0) temp_world["net_export"] = temp_world["total_exp"] - temp_world["total_imp"] temp_world["centroid"] = temp_world.geometry.centroid @@ -242,11 +358,21 @@ def process_data_for_plot(trade_file, world_gdf, value_column): # ============================================================================== # 4. PLOT # ============================================================================== -def plot(data_plot, year, product_name, value_column, colors, output_path=None, output_path_pdf=None): +def plot( + data_plot, + year, + product_name, + value_column, + colors, + output_path=None, + output_path_pdf=None, +): fig, ax = plt.subplots(1, 1, figsize=(18, 9)) world_df, flow_df = data_plot - world_df.plot(ax=ax, color=colors["land"], edgecolor=colors["border"], linewidth=0.4) + world_df.plot( + ax=ax, color=colors["land"], edgecolor=colors["border"], linewidth=0.4 + ) global_max_bubble = world_df["net_export"].abs().max() global_max_line = flow_df["quantity"].max() if len(flow_df) > 0 else 0 @@ -260,19 +386,26 @@ def plot(data_plot, year, product_name, value_column, colors, output_path=None, def get_radius(val): if abs(val) <= 0 or global_max_bubble == 0: return 0 - return min_r + (max_r - min_r) * (np.sqrt(abs(val)) / np.sqrt(global_max_bubble)) + return min_r + (max_r - min_r) * ( + np.sqrt(abs(val)) / np.sqrt(global_max_bubble) + ) # --- trade flow lines (thickness ∝ quantity) --- for _, row in flow_df.iterrows(): - if row["exporter"] in loc_map and row["importer"] in loc_map and global_max_line > 0: + if ( + row["exporter"] in loc_map + and row["importer"] in loc_map + and global_max_line > 0 + ): p1, p2 = loc_map[row["exporter"]], loc_map[row["importer"]] lw = (row["quantity"] / global_max_line) * 5.5 ax.plot( - [p1.x, p2.x], [p1.y, p2.y], + [p1.x, p2.x], + [p1.y, p2.y], color=colors["flow_line"], linewidth=max(lw, 0.3), alpha=0.35, - zorder=2 + zorder=2, ) # --- net-export / net-import bubbles (area ∝ quantity) --- @@ -289,28 +422,41 @@ def get_radius(val): alpha=0.7, edgecolor="white", linewidth=0.5, - zorder=5 + zorder=5, ) ) ax.set_title( - f"Global {product_name} Trade ({year})", - fontsize=22, - fontweight="bold", - pad=20 + f"Global {product_name} Trade ({year})", fontsize=22, fontweight="bold", pad=20 ) ax.axis("off") # ------------------------------------------------------------------ legend # 1) Category legend (surplus / deficit / flow line) cat_elements = [ - Line2D([0], [0], marker='o', color='w', label='Net exporter', - markerfacecolor=colors["surplus"], markeredgecolor='white', - markersize=12, alpha=0.7), - Line2D([0], [0], marker='o', color='w', label='Net importer', - markerfacecolor=colors["deficit"], markeredgecolor='white', - markersize=12, alpha=0.7), - Line2D([0], [0], color=colors["flow_line"], lw=2, label='Trade flow'), + Line2D( + [0], + [0], + marker="o", + color="w", + label="Net exporter", + markerfacecolor=colors["surplus"], + markeredgecolor="white", + markersize=12, + alpha=0.7, + ), + Line2D( + [0], + [0], + marker="o", + color="w", + label="Net importer", + markerfacecolor=colors["deficit"], + markeredgecolor="white", + markersize=12, + alpha=0.7, + ), + Line2D([0], [0], color=colors["flow_line"], lw=2, label="Trade flow"), ] leg1 = ax.legend( handles=cat_elements, @@ -329,78 +475,104 @@ def get_radius(val): # Stacked vertically (graduated-circles style) in the upper-left corner, # below the category legend. if global_max_bubble > 0: - unit = "Mt" if value_column == "quantity" else "M USD" - denom = 1e6 if value_column == "quantity" else 1e6 - scale_fracs = [1.0, 0.50, 0.25] # largest on top - scale_labels = [f"{frac * global_max_bubble / denom:.0f} {unit}" - for frac in scale_fracs] + unit = "Mt" if value_column == "quantity" else "M USD" + denom = 1e6 if value_column == "quantity" else 1e6 + scale_fracs = [1.0, 0.50, 0.25] # largest on top + scale_labels = [ + f"{frac * global_max_bubble / denom:.0f} {unit}" for frac in scale_fracs + ] map_height = bounds[3] - bounds[1] # anchor: upper-left, leave a small margin - margin_x = map_width * 0.03 + margin_x = map_width * 0.03 margin_y = map_height * 0.03 cx = bounds[0] + margin_x + max_r * 1.6 # centre x of all circles # stack circles top-to-bottom: largest first at the top - title_h = max_r * 1.1 - pad_y = max_r * 0.4 - y_start = bounds[3] - margin_y - title_h # top of first circle centre + title_h = max_r * 1.1 + pad_y = max_r * 0.4 + y_start = bounds[3] - margin_y - title_h # top of first circle centre centres_y = [] - y_cursor = y_start + y_cursor = y_start for frac in scale_fracs: r_i = get_radius(frac * global_max_bubble) - y_cursor -= r_i # move down by radius + y_cursor -= r_i # move down by radius centres_y.append(y_cursor) - y_cursor -= r_i + max_r * 0.25 # gap between circles + y_cursor -= r_i + max_r * 0.25 # gap between circles # extra room to the right for labels — use 3× max_r so values fit - pad_x_left = max_r * 0.6 + pad_x_left = max_r * 0.6 pad_x_right = max_r * 3.2 - box_left = cx - max_r - pad_x_left - box_top = bounds[3] - margin_y + box_left = cx - max_r - pad_x_left + box_top = bounds[3] - margin_y box_bottom = centres_y[-1] - max_r - pad_y - box_width = max_r + pad_x_left + pad_x_right + box_width = max_r + pad_x_left + pad_x_right box_height = box_top - box_bottom from matplotlib.patches import FancyBboxPatch - ax.add_patch(FancyBboxPatch( - (box_left, box_bottom), box_width, box_height, - boxstyle="round,pad=0.02", - facecolor=colors["legend_bg"], - edgecolor=colors["legend_edge"], - linewidth=0.8, zorder=8, - )) + + ax.add_patch( + FancyBboxPatch( + (box_left, box_bottom), + box_width, + box_height, + boxstyle="round,pad=0.02", + facecolor=colors["legend_bg"], + edgecolor=colors["legend_edge"], + linewidth=0.8, + zorder=8, + ) + ) # title at the top of the box - ax.text(box_left + box_width / 2, - box_top - pad_y * 0.5, - "Net trade\nvolume", - ha="center", va="top", fontsize=9, zorder=9) + ax.text( + box_left + box_width / 2, + box_top - pad_y * 0.5, + "Net trade\nvolume", + ha="center", + va="top", + fontsize=9, + zorder=9, + ) for frac, label, cy_i in zip(scale_fracs, scale_labels, centres_y): r_i = get_radius(frac * global_max_bubble) - ax.add_patch(Circle( - (cx, cy_i), r_i, - facecolor=colors["surplus"], alpha=0.7, - edgecolor="white", linewidth=0.5, zorder=9, - )) + ax.add_patch( + Circle( + (cx, cy_i), + r_i, + facecolor=colors["surplus"], + alpha=0.7, + edgecolor="white", + linewidth=0.5, + zorder=9, + ) + ) # label to the right of each circle - ax.text(cx + max_r + pad_x_left * 0.5, cy_i, label, - ha="left", va="center", fontsize=8, zorder=9) + ax.text( + cx + max_r + pad_x_left * 0.5, + cy_i, + label, + ha="left", + va="center", + fontsize=8, + zorder=9, + ) # 3) Line-width scale: show 3 reference flow sizes if global_max_line > 0: - unit = "Mt" if value_column == "quantity" else "M USD" - denom = 1e6 if value_column == "quantity" else 1e6 + unit = "Mt" if value_column == "quantity" else "M USD" + denom = 1e6 if value_column == "quantity" else 1e6 line_fracs = [0.25, 0.50, 1.0] line_handles = [] for frac in line_fracs: - lw = max(frac * 5.5, 0.3) + lw = max(frac * 5.5, 0.3) label = f"{frac * global_max_line / denom:.0f} {unit}" line_handles.append( - Line2D([0], [0], color=colors["flow_line"], lw=lw, - alpha=0.7, label=label) + Line2D( + [0], [0], color=colors["flow_line"], lw=lw, alpha=0.7, label=label + ) ) leg3 = ax.legend( handles=line_handles, @@ -441,9 +613,21 @@ def get_radius(val): # Map product names to snakemake output paths output_map = { - "Iron Ore": {"pdf": str(snakemake.output.iron_ore), "png": str(snakemake.output.iron_ore_png), "csv": str(snakemake.output.iron_ore_csv)}, - "DRI-HBI": {"pdf": str(snakemake.output.dri_hbi), "png": str(snakemake.output.dri_hbi_png), "csv": str(snakemake.output.dri_hbi_csv)}, - "Steel (raw)": {"pdf": str(snakemake.output.steel_raw), "png": str(snakemake.output.steel_raw_png), "csv": str(snakemake.output.steel_raw_csv)}, + "Iron Ore": { + "pdf": str(snakemake.output.iron_ore), + "png": str(snakemake.output.iron_ore_png), + "csv": str(snakemake.output.iron_ore_csv), + }, + "DRI-HBI": { + "pdf": str(snakemake.output.dri_hbi), + "png": str(snakemake.output.dri_hbi_png), + "csv": str(snakemake.output.dri_hbi_csv), + }, + "Steel": { + "pdf": str(snakemake.output.steel), + "png": str(snakemake.output.steel_png), + "csv": str(snakemake.output.steel_csv), + }, } print("Loading world geometry") diff --git a/workflow/notebooks/prepare-chokepoints.ipynb b/workflow/notebooks/prepare-chokepoints.ipynb index bbf35ae..9d093d6 100644 --- a/workflow/notebooks/prepare-chokepoints.ipynb +++ b/workflow/notebooks/prepare-chokepoints.ipynb @@ -83,7 +83,8 @@ "region_coords = {\n", " row[\"region_name\"]: [row[\"long\"], row[\"lat\"]]\n", " for _, row in bus_locations.iterrows()\n", - "}\n" + "}\n", + "region_display_names = snakemake.config.get(\"region_nice_names\", {})\n" ] }, { @@ -388,6 +389,11 @@ "fig = plt.figure(figsize=(10, 7))\n", "ax = plt.axes(projection=ccrs.Robinson())\n", "\n", + "cp_display_names = {\n", + " \"ormuz\": \"Hormuz\",\n", + " \"babalmandab\": \"Bab el-Mandeb\",\n", + "}\n", + "\n", "# Map features\n", "ax.add_feature(cfeature.LAND, facecolor=\"#f0f0f0\", edgecolor=\"none\")\n", "ax.add_feature(cfeature.OCEAN, facecolor=\"#dbe9f4\")\n", @@ -415,7 +421,8 @@ "for name, (lon, lat) in region_coords.items():\n", " ax.plot(lon, lat, marker=\"o\", color=\"#1f77b4\", markersize=5,\n", " transform=ccrs.PlateCarree(), zorder=10)\n", - " ax.text(lon + 2, lat + 2, name.replace(\"_\", \" \"), fontsize=5,\n", + " label = region_display_names.get(name, name.replace(\"_\", \" \"))\n", + " ax.text(lon + 2, lat + 2, label, fontsize=5,\n", " transform=ccrs.PlateCarree(), zorder=10,\n", " ha=\"left\", va=\"bottom\", color=\"#333333\")\n", "\n", @@ -424,7 +431,7 @@ " ax.plot(lon, lat, marker=\"*\", color=\"#d62728\", markersize=12,\n", " markeredgecolor=\"white\", markeredgewidth=0.5,\n", " transform=ccrs.PlateCarree(), zorder=20)\n", - " label = cp_name.replace(\"_\", \" \").title()\n", + " label = cp_display_names.get(cp_name, cp_name.replace(\"_\", \" \").title())\n", " ax.text(lon + 3, lat - 3, label, fontsize=6, fontweight=\"bold\",\n", " transform=ccrs.PlateCarree(), zorder=20,\n", " ha=\"left\", va=\"top\", color=\"#d62728\",\n", @@ -440,13 +447,13 @@ "]\n", "ax.legend(handles=legend_elements, loc=\"lower left\", fontsize=8, framealpha=0.9)\n", "\n", - "ax.set_title(\"Global Shipping Routes and Maritime Chokepoints\")\n", + "# ax.set_title(\"Global Shipping Routes and Maritime Chokepoints\")\n", "\n", "plt.tight_layout()\n", "plt.savefig(snakemake.output.map_chokepoints, dpi=300, bbox_inches=\"tight\")\n", "plt.savefig(snakemake.output.map_chokepoints_png, dpi=300, bbox_inches=\"tight\")\n", "plt.show()\n", - "print(f\"\\nRoutes through chokepoints: {sum(1 for r in routes_geo if r['has_chokepoint'])} / {len(routes_geo)}\")\n" + "print(f\"\\nRoutes through chokepoints: {sum(1 for r in routes_geo if r['has_chokepoint'])} / {len(routes_geo)}\")" ] } ], diff --git a/workflow/notebooks/prepare-wacc.ipynb b/workflow/notebooks/prepare-wacc.ipynb index 8663705..469cd2b 100644 --- a/workflow/notebooks/prepare-wacc.ipynb +++ b/workflow/notebooks/prepare-wacc.ipynb @@ -45,6 +45,7 @@ "\n", "# Output\n", "wacc_clustered_fn = snakemake.output.wacc\n", + "wacc_latex_fn = snakemake.output.wacc_latex\n", "\n", "# config\n", "config = snakemake.config" @@ -301,6 +302,77 @@ "wacc_by_region.to_csv(wacc_clustered_fn)" ] }, + { + "cell_type": "markdown", + "id": "5a1444f1", + "metadata": {}, + "source": [ + "### Save latex table" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cc7209fe", + "metadata": {}, + "outputs": [], + "source": [ + "wacc_table = (\n", + " wacc_by_region[[\"wacc_real\"]]\n", + " .rename(columns={\"wacc_real\": \"Weighted average cost of capital in \\%\"})\n", + " .rename(index=lambda region: config[\"region_nice_names\"].get(region, region)).mul(100)\n", + " .round(1)\n", + ")\n", + "\n", + "# wacc_table *= 100 # convert to percentage\n", + "wacc_table = wacc_table.rename(index={\"global_weighted_avg\": \"Global weighted average\"})\n", + "\n", + "# Flatten index so Region and value header are on the same header row in LaTeX.\n", + "wacc_table_out = wacc_table.rename_axis(\"Region\").reset_index()\n", + "\n", + "\n", + "def escape_latex(text):\n", + " return str(text).replace(\"&\", r\"\\&\")\n", + "\n", + "\n", + "value_col = wacc_table_out.columns[1]\n", + "column_name = escape_latex(value_col)\n", + "newline = r\"\\\\\"\n", + "latex_lines = [\n", + " r\"\\begin{table}[htbp]\",\n", + " r\"\\centering\",\n", + " r\"\\caption{Regional weighted average cost of capital (real), in \\%.}\",\n", + " r\"\\begin{tabular}{lr}\",\n", + " r\"\\toprule\",\n", + " f\"\\\\textbf{{Region}} & \\\\textbf{{{column_name}}} {newline}\",\n", + " r\"\\midrule\",\n", + "]\n", + "\n", + "for _, row in wacc_table_out.iterrows():\n", + " region = escape_latex(row[\"Region\"])\n", + " value = row[value_col]\n", + " latex_lines.append(f\"{region} & {value:.1f} {newline}\")\n", + "\n", + "latex_lines.extend([\n", + " r\"\\bottomrule\",\n", + " r\"\\end{tabular}\",\n", + " r\"\\end{table}\",\n", + "])\n", + "\n", + "with open(wacc_latex_fn, \"w\", encoding=\"utf-8\") as f:\n", + " f.write(\"\\n\".join(latex_lines) + \"\\n\")\n", + "\n", + "wacc_table_out" + ] + }, + { + "cell_type": "markdown", + "id": "6b852fc9", + "metadata": {}, + "source": [ + "### Further analysis" + ] + }, { "cell_type": "code", "execution_count": null, @@ -308,7 +380,17 @@ "metadata": {}, "outputs": [], "source": [ - "wacc_by_region" + "wacc_by_region.loc[:,[\"wacc_real\", \"wacc_spread\"]].round(3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d0af51fc", + "metadata": {}, + "outputs": [], + "source": [ + "wacc.loc[wacc.region.isin([\"Middle_East\",\"North_West_Africa\"])].loc[:,[\"country_name\", \"wacc_real\", \"wacc\", \"region\"]].sort_values(\"region\").round(3)" ] } ], diff --git a/workflow/scripts/model_trade.py b/workflow/scripts/model_trade.py index de90570..20ffd0e 100644 --- a/workflow/scripts/model_trade.py +++ b/workflow/scripts/model_trade.py @@ -20,6 +20,41 @@ plt.style.use("bmh") +def add_link_cost_column(df, quantity_col, cost_col, link_cost_col): + """ + Convert cumulative quantity / average cost data into stepwise link costs. + + Example: + - 50 Mt at 500 EUR/t + - 100 Mt at 1000 EUR/t + + becomes: + - first 50 Mt at 500 EUR/t + - next 50 Mt at 1500 EUR/t + """ + + df = df.copy() + + quantity = df[quantity_col].astype(float) + average_cost = df[cost_col].astype(float) + + previous_quantity = quantity.shift(fill_value=0.0) + total_cost = quantity * average_cost + previous_total_cost = total_cost.shift(fill_value=0.0) + + incremental_quantity = quantity - previous_quantity + incremental_total_cost = total_cost - previous_total_cost + + if (incremental_quantity <= 0).any(): + raise ValueError( + f"Column '{quantity_col}' must be strictly increasing to derive link_cost." + ) + + df[link_cost_col] = incremental_total_cost / incremental_quantity + + return df + + def build_region_geodataframe(config): """ Build a dissolved GeoDataFrame of model regions from the config country lists. @@ -110,11 +145,23 @@ def building_model( .dropna(subset=[f"{cost_descriptor} [EUR/{unit}]"]) .reset_index(drop=True) ) # Filter out rows where "lcox [EUR/t]" is NaN + region_data_interone = add_link_cost_column( + region_data_interone, + quantity_col=f"demand [{unit}]", + cost_col=f"{cost_descriptor} [EUR/{unit}]", + link_cost_col=f"link_cost [EUR/{unit}]", + ) region_data_intertwo = ( pd.read_csv(region_file_intertwo, header=0) .dropna(subset=[f"{cost_descriptor} [EUR/{unit}]"]) .reset_index(drop=True) ) # Filter out rows where "lcox [EUR/t]" is NaN + region_data_intertwo = add_link_cost_column( + region_data_intertwo, + quantity_col=f"demand [{unit}]", + cost_col=f"{cost_descriptor} [EUR/{unit}]", + link_cost_col=f"link_cost [EUR/{unit}]", + ) region_name = os.path.basename(region_file_interone).rsplit( "_marginal_cost_", 1 @@ -190,7 +237,7 @@ def building_model( region_data_interone[f"demand [{unit}]"][s] ) - float(region_data_interone[f"demand [{unit}]"][s - 1]) - m_cost = float(region_data_interone[f"{cost_descriptor} [EUR/{unit}]"][s]) + m_cost = float(region_data_interone[f"link_cost [EUR/{unit}]"][s]) if not is_material_chain: # Pure energy chain (e.g. hydrogen): Generator directly on final bus @@ -249,9 +296,7 @@ def building_model( region_data_intertwo ) # split evenly across supply steps - m_cost = float( - region_data_intertwo[f"{cost_descriptor} [EUR/{unit}]"][s] - ) + m_cost = float(region_data_intertwo[f"link_cost [EUR/{unit}]"][s]) n.add( "Link", @@ -555,6 +600,7 @@ def plot_trade_network( bus_colors=supply_color, bus_alpha=alpha_supply, link_widths=trade * plot_config["link_width"], + link_alpha=0.7, branch_components=["Link"], link_colors=link_colors, geomap=False, @@ -1255,7 +1301,7 @@ def solve_pareto(n, pareto_config, interone, final, is_material_chain): interone="hbi", intertwo="eaf", final="steel", - scenario="constrain-supply", + scenario="default", wacc="regional", chain_id="supplyconstraint", )