Expand analysis of results

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'

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,

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,16 +173,54 @@ 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",
     notebook:
         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)

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": {