move building of distribution of existing heating to own script

This makes the distribution of existing heating to urban/rural,
residential/services and spatially more transparent.

rules/solve_myopic.smk

@@ -1,8 +1,40 @@
-# SPDX-FileCopyrightText: : 2023 The PyPSA-Eur Authors
+# SPDX-FileCopyrightText: : 2023-4 The PyPSA-Eur Authors
 #
 # SPDX-License-Identifier: MIT
 
 
+rule build_existing_heating_distribution:
+    params:
+        baseyear=config["scenario"]["planning_horizons"][0],
+        sector=config["sector"],
+        existing_capacities=config["existing_capacities"],
+    input:
+        existing_heating="data/existing_infrastructure/existing_heating_raw.csv",
+        clustered_pop_layout=RESOURCES + "pop_layout_elec_s{simpl}_{clusters}.csv",
+        clustered_pop_energy_layout=RESOURCES + "pop_weighted_energy_totals_s{simpl}_{clusters}.csv",
+        district_heat_share=RESOURCES + "district_heat_share_elec_s{simpl}_{clusters}_{planning_horizons}.csv",
+    output:
+        existing_heating_distribution=RESOURCES
+        + "existing_heating_distribution_elec_s{simpl}_{clusters}_{planning_horizons}.csv",
+    wildcard_constraints:
+        planning_horizons=config["scenario"]["planning_horizons"][0],  #only applies to baseyear
+    threads: 1
+    resources:
+        mem_mb=2000,
+    log:
+        LOGS
+        + "build_existing_heating_distribution_elec_s{simpl}_{clusters}_{planning_horizons}.log",
+    benchmark:
+        (
+            BENCHMARKS
+            + "build_existing_heating_distribution/elec_s{simpl}_{clusters}_{planning_horizons}"
+        )
+    conda:
+        "../envs/environment.yaml"
+    script:
+        "../scripts/build_existing_heating_distribution.py"
+
+
 rule add_existing_baseyear:
     params:
         baseyear=config["scenario"]["planning_horizons"][0],
@@ -19,7 +51,8 @@ rule add_existing_baseyear:
         costs="data/costs_{}.csv".format(config["scenario"]["planning_horizons"][0]),
         cop_soil_total=RESOURCES + "cop_soil_total_elec_s{simpl}_{clusters}.nc",
         cop_air_total=RESOURCES + "cop_air_total_elec_s{simpl}_{clusters}.nc",
-        existing_heating="data/existing_infrastructure/existing_heating_raw.csv",
+        existing_heating_distribution=RESOURCES
+        + "existing_heating_distribution_elec_s{simpl}_{clusters}_{planning_horizons}.csv",
         existing_solar="data/existing_infrastructure/solar_capacity_IRENA.csv",
         existing_onwind="data/existing_infrastructure/onwind_capacity_IRENA.csv",
         existing_offwind="data/existing_infrastructure/offwind_capacity_IRENA.csv",

scripts/add_existing_baseyear.py

@@ -409,97 +409,20 @@ def add_heating_capacities_installed_before_baseyear(
     # file: "WP2_DataAnnex_1_BuildingTechs_ForPublication_201603.xls" -> "existing_heating_raw.csv".
     # TODO start from original file
 
-    # retrieve existing heating capacities
+    existing_heating = pd.read_csv(snakemake.input.existing_heating_distribution,
-    techs = [
+                                   header=[0,1],
-        "gas boiler",
+                                   index_col=0)
-        "oil boiler",
-        "resistive heater",
-        "air heat pump",
-        "ground heat pump",
-    ]
-    df = pd.read_csv(snakemake.input.existing_heating, index_col=0, header=0)
 
-    # data for Albania, Montenegro and Macedonia not included in database
-    df.loc["Albania"] = np.nan
-    df.loc["Montenegro"] = np.nan
-    df.loc["Macedonia"] = np.nan
 
-    df.fillna(0.0, inplace=True)
+    techs = existing_heating.columns.get_level_values(1).unique()
 
-    # convert GW to MW
+    for name in existing_heating.columns.get_level_values(0).unique():
-    df *= 1e3
 
-    df.index = cc.convert(df.index, to="iso2")
-
-    # coal and oil boilers are assimilated to oil boilers
-    df["oil boiler"] = df["oil boiler"] + df["coal boiler"]
-    df.drop(["coal boiler"], axis=1, inplace=True)
-
-    # distribute technologies to nodes by population
-    pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0)
-
-    nodal_df = df.loc[pop_layout.ct]
-    nodal_df.index = pop_layout.index
-    nodal_df = nodal_df.multiply(pop_layout.fraction, axis=0)
-
-    # split existing capacities between residential and services
-    # proportional to energy demand
-    p_set_sum = n.loads_t.p_set.sum()
-    ratio_residential = pd.Series(
-        [
-            (
-                p_set_sum[f"{node} residential rural heat"]
-                / (
-                    p_set_sum[f"{node} residential rural heat"]
-                    + p_set_sum[f"{node} services rural heat"]
-                )
-            )
-            # if rural heating demand for one of the nodes doesn't exist,
-            # then columns were dropped before and heating demand share should be 0.0
-            if all(
-                f"{node} {service} rural heat" in p_set_sum.index
-                for service in ["residential", "services"]
-            )
-            else 0.0
-            for node in nodal_df.index
-        ],
-        index=nodal_df.index,
-    )
-
-    for tech in techs:
-        nodal_df["residential " + tech] = nodal_df[tech] * ratio_residential
-        nodal_df["services " + tech] = nodal_df[tech] * (1 - ratio_residential)
-
-    names = [
-        "residential rural",
-        "services rural",
-        "residential urban decentral",
-        "services urban decentral",
-        "urban central",
-    ]
-
-    nodes = {}
-    p_nom = {}
-    for name in names:
         name_type = "central" if name == "urban central" else "decentral"
-        nodes[name] = pd.Index(
-            [
-                n.buses.at[index, "location"]
-                for index in n.buses.index[
-                    n.buses.index.str.contains(name)
-                    & n.buses.index.str.contains("heat")
-                ]
-            ]
-        )
-        heat_pump_type = "air" if "urban" in name else "ground"
-        heat_type = "residential" if "residential" in name else "services"
 
-        if name == "urban central":
+        nodes = pd.Index(n.buses.location[n.buses.index.str.contains(f"{name} heat")])
-            p_nom[name] = nodal_df["air heat pump"][nodes[name]]
+
-        else:
+        heat_pump_type = "air" if "urban" in name else "ground"
-            p_nom[name] = nodal_df[f"{heat_type} {heat_pump_type} heat pump"][
-                nodes[name]
-            ]
 
         # Add heat pumps
         costs_name = f"decentral {heat_pump_type}-sourced heat pump"
@@ -507,7 +430,7 @@ def add_heating_capacities_installed_before_baseyear(
         cop = {"air": ashp_cop, "ground": gshp_cop}
 
         if time_dep_hp_cop:
-            efficiency = cop[heat_pump_type][nodes[name]]
+            efficiency = cop[heat_pump_type][nodes]
         else:
             efficiency = costs.at[costs_name, "efficiency"]
 
@@ -520,27 +443,26 @@ def add_heating_capacities_installed_before_baseyear(
 
             n.madd(
                 "Link",
-                nodes[name],
+                nodes,
                 suffix=f" {name} {heat_pump_type} heat pump-{grouping_year}",
-                bus0=nodes[name],
+                bus0=nodes,
-                bus1=nodes[name] + " " + name + " heat",
+                bus1=nodes + " " + name + " heat",
                 carrier=f"{name} {heat_pump_type} heat pump",
                 efficiency=efficiency,
                 capital_cost=costs.at[costs_name, "efficiency"]
                 * costs.at[costs_name, "fixed"],
-                p_nom=p_nom[name] * ratio / costs.at[costs_name, "efficiency"],
+                p_nom=existing_heating[(name, f"{heat_pump_type} heat pump")][nodes] * ratio / costs.at[costs_name, "efficiency"],
                 build_year=int(grouping_year),
                 lifetime=costs.at[costs_name, "lifetime"],
             )
 
             # add resistive heater, gas boilers and oil boilers
-            # (50% capacities to rural buses, 50% to urban buses)
             n.madd(
                 "Link",
-                nodes[name],
+                nodes,
                 suffix=f" {name} resistive heater-{grouping_year}",
-                bus0=nodes[name],
+                bus0=nodes,
-                bus1=nodes[name] + " " + name + " heat",
+                bus1=nodes + " " + name + " heat",
                 carrier=name + " resistive heater",
                 efficiency=costs.at[f"{name_type} resistive heater", "efficiency"],
                 capital_cost=(
@@ -548,21 +470,20 @@ def add_heating_capacities_installed_before_baseyear(
                     * costs.at[f"{name_type} resistive heater", "fixed"]
                 ),
                 p_nom=(
-                    0.5
+                    existing_heating[(name, "resistive heater")][nodes]
-                    * nodal_df[f"{heat_type} resistive heater"][nodes[name]]
                     * ratio
                     / costs.at[f"{name_type} resistive heater", "efficiency"]
                 ),
                 build_year=int(grouping_year),
-                lifetime=costs.at[costs_name, "lifetime"],
+                lifetime=costs.at[f"{name_type} resistive heater", "lifetime"],
             )
 
             n.madd(
                 "Link",
-                nodes[name],
+                nodes,
                 suffix=f" {name} gas boiler-{grouping_year}",
                 bus0=spatial.gas.nodes,
-                bus1=nodes[name] + " " + name + " heat",
+                bus1=nodes + " " + name + " heat",
                 bus2="co2 atmosphere",
                 carrier=name + " gas boiler",
                 efficiency=costs.at[f"{name_type} gas boiler", "efficiency"],
@@ -572,8 +493,7 @@ def add_heating_capacities_installed_before_baseyear(
                     * costs.at[f"{name_type} gas boiler", "fixed"]
                 ),
                 p_nom=(
-                    0.5
+                    existing_heating[(name, "gas boiler")][nodes]
-                    * nodal_df[f"{heat_type} gas boiler"][nodes[name]]
                     * ratio
                     / costs.at[f"{name_type} gas boiler", "efficiency"]
                 ),
@@ -583,20 +503,20 @@ def add_heating_capacities_installed_before_baseyear(
 
             n.madd(
                 "Link",
-                nodes[name],
+                nodes,
                 suffix=f" {name} oil boiler-{grouping_year}",
                 bus0=spatial.oil.nodes,
-                bus1=nodes[name] + " " + name + " heat",
+                bus1=nodes + " " + name + " heat",
                 bus2="co2 atmosphere",
                 carrier=name + " oil boiler",
                 efficiency=costs.at["decentral oil boiler", "efficiency"],
                 efficiency2=costs.at["oil", "CO2 intensity"],
                 capital_cost=costs.at["decentral oil boiler", "efficiency"]
                 * costs.at["decentral oil boiler", "fixed"],
-                p_nom=0.5
+                p_nom= (
-                * nodal_df[f"{heat_type} oil boiler"][nodes[name]]
+                    existing_heating[(name, "oil boiler")][nodes]
-                * ratio
+                    * ratio
-                / costs.at["decentral oil boiler", "efficiency"],
+                    / costs.at["decentral oil boiler", "efficiency"]),
                 build_year=int(grouping_year),
                 lifetime=costs.at[f"{name_type} gas boiler", "lifetime"],
             )
@@ -624,6 +544,8 @@ def add_heating_capacities_installed_before_baseyear(
 
             # drop assets which are at the end of their lifetime
             links_i = n.links[(n.links.build_year + n.links.lifetime <= baseyear)].index
+            logger.info("Removing following links because at end of their lifetime:")
+            logger.info(links_i)
             n.mremove("Link", links_i)
 
 

scripts/build_existing_heating_distribution.py

@@ -0,0 +1,108 @@
+# -*- coding: utf-8 -*-
+# SPDX-FileCopyrightText: : 2020-2023 The PyPSA-Eur Authors
+#
+# SPDX-License-Identifier: MIT
+"""
+Builds table of existing heat generation capacities for initial planning
+horizon.
+"""
+import pandas as pd
+import sys
+from pypsa.descriptors import Dict
+import numpy as np
+import country_converter as coco
+
+cc = coco.CountryConverter()
+
+
+def build_existing_heating():
+    # retrieve existing heating capacities
+    techs = [
+        "gas boiler",
+        "oil boiler",
+        "resistive heater",
+        "air heat pump",
+        "ground heat pump",
+    ]
+
+    existing_heating = pd.read_csv(snakemake.input.existing_heating,
+                                   index_col=0,
+                                   header=0)
+
+    # data for Albania, Montenegro and Macedonia not included in database                                                                          existing_heating.loc["Albania"] = np.nan
+    existing_heating.loc["Montenegro"] = np.nan
+    existing_heating.loc["Macedonia"] = np.nan
+
+    existing_heating.fillna(0.0, inplace=True)
+
+    # convert GW to MW
+    existing_heating *= 1e3
+
+    existing_heating.index = cc.convert(existing_heating.index, to="iso2")
+
+    # coal and oil boilers are assimilated to oil boilers
+    existing_heating["oil boiler"] = existing_heating["oil boiler"] + existing_heating["coal boiler"]
+    existing_heating.drop(["coal boiler"], axis=1, inplace=True)
+
+    # distribute technologies to nodes by population
+    pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout,
+                             index_col=0)
+
+    nodal_heating = existing_heating.loc[pop_layout.ct]
+    nodal_heating.index = pop_layout.index
+    nodal_heating = nodal_heating.multiply(pop_layout.fraction, axis=0)
+
+    district_heat_info = pd.read_csv(snakemake.input.district_heat_share,
+                                     index_col=0)
+    dist_fraction = district_heat_info["district fraction of node"]
+    urban_fraction = district_heat_info["urban fraction"]
+
+    energy_layout = pd.read_csv(snakemake.input.clustered_pop_energy_layout,
+                                index_col=0)
+
+    uses = ["space", "water"]
+    sectors = ["residential", "services"]
+
+    nodal_sectoral_totals = pd.DataFrame(dtype=float)
+
+    for sector in sectors:
+        nodal_sectoral_totals[sector] = energy_layout[[f"total {sector} {use}" for use in uses]].sum(axis=1)
+
+    nodal_sectoral_fraction = nodal_sectoral_totals.div(nodal_sectoral_totals.sum(axis=1),
+                                                        axis=0)
+
+
+    nodal_heat_name_fraction = pd.DataFrame(dtype=float)
+
+    nodal_heat_name_fraction["urban central"] = dist_fraction
+
+    for sector in sectors:
+
+        nodal_heat_name_fraction[f"{sector} rural"] = nodal_sectoral_fraction[sector]*(1 - urban_fraction)
+        nodal_heat_name_fraction[f"{sector} urban decentral"] = nodal_sectoral_fraction[sector]*(urban_fraction - dist_fraction)
+
+
+    nodal_heat_name_tech = pd.concat({name : nodal_heating .multiply(nodal_heat_name_fraction[name],
+                                                                     axis=0) for name in nodal_heat_name_fraction.columns},
+                                     axis=1,
+                                     names=["heat name","technology"])
+
+
+    #move all ground HPs to rural, all air to urban
+
+    for sector in sectors:
+        nodal_heat_name_tech[(f"{sector} rural","ground heat pump")] += (nodal_heat_name_tech[("urban central","ground heat pump")]*nodal_sectoral_fraction[sector]
+                                                                         + nodal_heat_name_tech[(f"{sector} urban decentral","ground heat pump")])
+        nodal_heat_name_tech[(f"{sector} urban decentral","ground heat pump")] = 0.
+
+        nodal_heat_name_tech[(f"{sector} urban decentral","air heat pump")] += nodal_heat_name_tech[(f"{sector} rural","air heat pump")]
+        nodal_heat_name_tech[(f"{sector} rural","air heat pump")] = 0.
+
+    nodal_heat_name_tech[("urban central","ground heat pump")] = 0.
+
+    nodal_heat_name_tech.to_csv(snakemake.output.existing_heating_distribution)
+
+
+if __name__ == "__main__":
+
+    build_existing_heating()