Merge remote-tracking branch 'upstream/master' into add_methanol_techs

2024-08-12 09:52:25 +02:00 · 2024-08-12 09:52:25 +02:00 · 77385721b6
commit 77385721b6
parent 280e8d6d01 f8d0efbe99
29 changed files with 1802 additions and 478 deletions
--- a/0
+++ b/0
--- a/config/config.default.yaml
+++ b/config/config.default.yaml
@ -67,7 +67,6 @@ snapshots:
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#enable
 enable:
  retrieve: auto
  prepare_links_p_nom: false
  retrieve_databundle: true
  retrieve_cost_data: true
  build_cutout: false
@ -355,7 +354,6 @@ biomass:
    - Secondary Forestry residues - woodchips
    - Sawdust
    - Residues from landscape care
    - Municipal waste
    not included:
    - Sugar from sugar beet
    - Rape seed
@ -369,6 +367,25 @@ biomass:
    biogas:
    - Manure solid, liquid
    - Sludge
    municipal solid waste:
    - Municipal waste
  share_unsustainable_use_retained:
    2020: 1
    2025: 0.66
    2030: 0.33
    2035: 0
    2040: 0
    2045: 0
    2050: 0
  share_sustainable_potential_available:
    2020: 0
    2025: 0.33
    2030: 0.66
    2035: 1
    2040: 1
    2045: 1
    2050: 1
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#solar-thermal
 solar_thermal:
@ -409,6 +426,22 @@ sector:
      2045: 0.8
      2050: 1.0
    district_heating_loss: 0.15
    forward_temperature: 90 #C
    return_temperature: 50 #C
    heat_source_cooling: 6 #K
    heat_pump_cop_approximation:
      refrigerant: ammonia
      heat_exchanger_pinch_point_temperature_difference: 5 #K
      isentropic_compressor_efficiency: 0.8
      heat_loss: 0.0
  heat_pump_sources:
    urban central:
    - air
    urban decentral:
    - air
    rural:
    - air
    - ground
  cluster_heat_buses: true
  heat_demand_cutout: default
  bev_dsm_restriction_value: 0.75
@ -491,7 +524,7 @@ sector:
  aviation_demand_factor: 1.
  HVC_demand_factor: 1.
  time_dep_hp_cop: true
-  heat_pump_sink_T: 55.
+  heat_pump_sink_T_individual_heating: 55.
  reduce_space_heat_exogenously: true
  reduce_space_heat_exogenously_factor:
    2020: 0.10  # this results in a space heat demand reduction of 10%
@ -613,6 +646,7 @@ sector:
  biomass_to_liquid: false
  electrobiofuels: false
  biosng: false
  municipal_solid_waste: false
  limit_max_growth:
    enable: false
    # allowing 30% larger than max historic growth
@ -640,6 +674,7 @@ sector:
    max_amount: 1390 # TWh
    upstream_emissions_factor: .1 #share of solid biomass CO2 emissions at full combustion
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#industry
 industry:
  St_primary_fraction:
@ -733,7 +768,7 @@ industry:
 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#costs
 costs:
  year: 2030
-  version: v0.9.0
+  version: v0.9.1
  social_discountrate: 0.02
  fill_values:
    FOM: 0
@ -1036,6 +1071,7 @@ plotting:
    biogas: '#e3d37d'
    biomass: '#baa741'
    solid biomass: '#baa741'
    municipal solid waste: '#91ba41'
    solid biomass import: '#d5ca8d'
    solid biomass transport: '#baa741'
    solid biomass for industry: '#7a6d26'
@ -1050,6 +1086,7 @@ plotting:
    services rural biomass boiler: '#c6cf98'
    services urban decentral biomass boiler: '#dde5b5'
    biomass to liquid: '#32CD32'
    unsustainable bioliquids: '#32CD32'
    electrobiofuels: 'red'
    BioSNG: '#123456'
    # power transmission
--- a/data/links_p_nom.csv
+++ b/data/links_p_nom.csv
@ -5,7 +5,7 @@ Cross-Channel,France - Echingen 50°41′48″N 1°38′21″E / 50.69667
 Volgograd-Donbass,Russia - Volzhskaya 48°49′34″N 44°40′20″E / 48.82611°N 44.67222°E,Ukraine - Mikhailovskaya 48°39′13″N 38°33′56″E / 48.65361°N 38.56556°E,475(0/475),400,750.0,1965,Merc/Thyr,Shut down in 2014,[1],44.672222222222224,48.82611111111111,38.565555555555555,48.65361111111111
 Konti-Skan 1,Denmark - Vester Hassing 57°3′46″N 10°5′24″E / 57.06278°N 10.09000°E,Sweden - Stenkullen 57°48′15″N 12°19′13″E / 57.80417°N 12.32028°E,176(87/89),250,250.0,1965,Merc,Replaced in August 2006 by modern converters using thyristors,[1],10.09,57.062777777777775,12.320277777777777,57.80416666666667
 SACOI 1a,Italy - Suvereto 43°3′10″N 10°41′42″E / 43.05278°N 10.69500°E ( before 1992: Italy - San Dalmazio 43°15′43″N 10°55′05″E / 43.26194°N 10.91806°E),"France- Lucciana 42°31′40″N 9°26′59″E / 42.52778°N 9.44972°E",483(365/118),200,200.0,1965,Merc,"Replaced in 1986 by Thyr- multiterminal scheme",[1],10.695,43.05277777777778,9.449722222222222,42.52777777777778
-SACOI 1b,"France- Lucciana 42°31′40″N 9°26′59″E / 42.52778°N 9.44972°E", "Codrongianos- Italy 40°39′7″N 8°42′48″E / 40.65194°N 8.71333°E",483(365/118),200,200.0,1965,Merc,"Replaced in 1986 by Thyr- multiterminal scheme",[1],9.449722222222222,42.52777777777778,8.679351,40.65765
+SACOI 1b,"France- Lucciana 42°31′40″N 9°26′59″E / 42.52778°N 9.44972°E","Codrongianos- Italy 40°39′7″N 8°42′48″E / 40.65194°N 8.71333°E",483(365/118),200,200.0,1965,Merc,"Replaced in 1986 by Thyr- multiterminal scheme",[1],9.449722222222222,42.52777777777778,8.679351,40.65765
 Kingsnorth,UK - Kingsnorth 51°25′11″N 0°35′46″E / 51.41972°N 0.59611°E,UK - London-Beddington 51°22′23″N 0°7′38″W / 51.37306°N 0.12722°W,85(85/0),266,320.0,1975,Merc,Bipolar scheme Supplier: English Electric Shut down in 1987,[33],0.5961111111111111,51.41972222222222,-0.1272222222222222,51.37305555555555
 Skagerrak 1 + 2,Denmark - Tjele 56°28′44″N 9°34′1″E / 56.47889°N 9.56694°E,Norway - Kristiansand 58°15′36″N 7°53′55″E / 58.26000°N 7.89861°E,230(130/100),250,500.0,1977,Thyr,Supplier: STK(Nexans) Control system upgrade by ABB in 2007,[34][35][36],9.566944444444445,56.47888888888889,7.898611111111111,58.26
 Gotland 2,Sweden - Västervik 57°43′41″N 16°38′51″E / 57.72806°N 16.64750°E,Sweden - Yigne 57°35′13″N 18°11′44″E / 57.58694°N 18.19556°E,99.5(92.9/6.6),150,130.0,1983,Thyr,Supplier: ABB,,16.6475,57.72805555555556,18.195555555555554,57.58694444444444
@ -23,7 +23,7 @@ Visby-Nas,Sweden - Nas 57°05′58″N 18°14′27″E / 57.09944°N 18.24
 SwePol,Poland - Wierzbięcin 54°30′8″N 16°53′28″E / 54.50222°N 16.89111°E,Sweden - Stärnö 56°9′11″N 14°50′29″E / 56.15306°N 14.84139°E,245(245/0),450,600.0,2000,Thyr,Supplier: ABB,[38],16.891111111111112,54.50222222222222,14.841388888888888,56.153055555555554
 Tjæreborg,Denmark - Tjæreborg/Enge 55°26′52″N 8°35′34″E / 55.44778°N 8.59278°E,Denmark - Tjæreborg/Substation 55°28′07″N 8°33′36″E / 55.46861°N 8.56000°E,4.3(4.3/0),9,7.0,2000,IGBT,Interconnection to wind power generating stations,,8.592777777777778,55.44777777777778,8.56,55.46861111111111
 Italy-Greece,Greece - Arachthos 39°11′00″N 20°57′48″E / 39.18333°N 20.96333°E,Italy - Galatina 40°9′53″N 18°7′49″E / 40.16472°N 18.13028°E,310(200/110),400,500.0,2001,Thyr,,,20.963333333333335,39.18333333333333,18.130277777777778,40.164722222222224
-Moyle,UK - Auchencrosh 55°04′10″N 4°58′50″W / 55.06944°N 4.98056°W,UK - N. Ireland- Ballycronan More 54°50′34″N 5°46′11″W / 54.84278°N 5.76972°W,63.5(63.5/0),250,2501.0,2001,Thyr,"Supplier: Siemens- Nexans",[39],-4.980555555555556,55.06944444444444,-5.769722222222223,54.842777777777776
+Moyle,UK - Auchencrosh 55°04′10″N 4°58′50″W / 55.06944°N 4.98056°W,UK - N. Ireland- Ballycronan More 54°50′34″N 5°46′11″W / 54.84278°N 5.76972°W,63.5(63.5/0),250,500.0,2001,Thyr,"Supplier: Siemens- Nexans",[39],-4.980555555555556,55.06944444444444,-5.769722222222223,54.842777777777776
 HVDC Troll,Norway - Kollsnes 60°33′01″N 4°50′26″E / 60.55028°N 4.84056°E,Norway - Offshore platform Troll A 60°40′00″N 3°40′00″E / 60.66667°N 3.66667°E,70(70/0),60,80.0,2004,IGBT,Power supply for offshore gas compressor Supplier: ABB,[40],4.8405555555555555,60.55027777777778,3.6666666666666665,60.666666666666664
 Estlink,Finland - Espoo 60°12′14″N 24°33′06″E / 60.20389°N 24.55167°E,Estonia - Harku 59°23′5″N 24°33′37″E / 59.38472°N 24.56028°E,105(105/0),150,350.0,2006,IGBT,Supplier: ABB,[40],24.551666666666666,60.20388888888889,24.560277777777777,59.38472222222222
 NorNed,Netherlands - Eemshaven 53°26′4″N 6°51′57″E / 53.43444°N 6.86583°E,Norway - Feda 58°16′58″N 6°51′55″E / 58.28278°N 6.86528°E,580(580/0),450,700.0,2008,Thyr,"Supplier: ABB- Nexans",[40],6.865833333333334,53.434444444444445,6.865277777777778,58.28277777777778
--- a/doc/conf.py
+++ b/doc/conf.py
@ -53,7 +53,6 @@ extensions = [
 autodoc_mock_imports = [
    "atlite",
    "snakemake",
    "pycountry",
    "rioxarray",
    "country_converter",
    "tabula",
@ -341,4 +340,6 @@ texinfo_documents = [
 # Example configuration for intersphinx: refer to the Python standard library.
-intersphinx_mapping = {"https://docs.python.org/": None}
+intersphinx_mapping = {
    "https://docs.python.org/": ("https://docs.python.org/3", None),
 }
--- a/doc/configtables/biomass.csv
+++ b/doc/configtables/biomass.csv
@ -5,3 +5,5 @@ classes ,,,
 -- solid biomass,--,Array of biomass comodity,The comodity that are included as solid biomass
 -- not included,--,Array of biomass comodity,The comodity that are not included as a biomass potential
 -- biogas,--,Array of biomass comodity,The comodity that are included as biogas
 share_unsustainable_use_retained,--,Dictionary with planning horizons as keys., Share of unsustainable biomass use retained using primary production of Eurostat data as reference
 share_sustainable_potential_available,--,Dictionary with planning horizons as keys., Share determines phase-in of ENSPRESO biomass potentials
--- a/doc/configtables/enable.csv
+++ b/doc/configtables/enable.csv
@ -1,6 +1,5 @@
 ,Unit,Values,Description
 enable,str or bool,"{auto, true, false}","Switch to include (true) or exclude (false) the retrieve_* rules of snakemake into the workflow; 'auto' sets true|false based on availability of an internet connection to prevent issues with snakemake failing due to lack of internet connection."
 prepare_links_p_nom,bool,"{true, false}","Switch to retrieve current HVDC projects from `Wikipedia <https://en.wikipedia.org/wiki/List_of_HVDC_projects>`_"
 retrieve_databundle,bool,"{true, false}","Switch to retrieve databundle from zenodo via the rule :mod:`retrieve_databundle` or whether to keep a custom databundle located in the corresponding folder."
 retrieve_cost_data,bool,"{true, false}","Switch to retrieve technology cost data from `technology-data repository <https://github.com/PyPSA/technology-data>`_."
 build_cutout,bool,"{true, false}","Switch to enable the building of cutouts via the rule :mod:`build_cutout`."
--- a/doc/configtables/sector.csv
+++ b/doc/configtables/sector.csv
@ -9,6 +9,18 @@ district_heating,--,,`prepare_sector_network.py <https://github.com/PyPSA/pypsa-
 -- potential,--,float,maximum fraction of urban demand which can be supplied by district heating
 -- progress,--,Dictionary with planning horizons as keys., Increase of today's district heating demand to potential maximum district heating share. Progress = 0 means today's district heating share. Progress = 1 means maximum fraction of urban demand is supplied by district heating
 -- district_heating_loss,--,float,Share increase in district heat demand in urban central due to heat losses
 -- forward_temperature,°C,float,Forward temperature in district heating
 -- return_temperature,°C,float,Return temperature in district heating. Must be lower than forward temperature
 -- heat_source_cooling,K,float,Cooling of heat source for heat pumps
 -- heat_pump_cop_approximation,,,
 -- -- refrigerant,--,"{ammonia, isobutane}",Heat pump refrigerant assumed for COP approximation
 -- -- heat_exchanger_pinch_point_temperature_difference,K,float,Heat pump pinch point temperature difference in heat exchangers assumed for approximation.
 -- -- isentropic_compressor_efficiency,--,float,Isentropic efficiency of heat pump compressor assumed for approximation. Must be between 0 and 1.
 -- -- heat_loss,--,float,Heat pump heat loss assumed for approximation. Must be between 0 and 1.
 -- heat_pump_sources,--,,
 -- -- urban central,--,List of heat sources for heat pumps in urban central heating,
 -- -- urban decentral,--,List of heat sources for heat pumps in urban decentral heating,
 -- -- rural,--,List of heat sources for heat pumps in rural heating,
 cluster_heat_buses,--,"{true, false}",Cluster residential and service heat buses in `prepare_sector_network.py <https://github.com/PyPSA/pypsa-eur-sec/blob/master/scripts/prepare_sector_network.py>`_  to one to save memory.
 ,,,
 bev_dsm_restriction _value,--,float,Adds a lower state of charge (SOC) limit for battery electric vehicles (BEV) to manage its own energy demand (DSM). Located in `build_transport_demand.py <https://github.com/PyPSA/pypsa-eur-sec/blob/master/scripts/build_transport_demand.py>`_. Set to 0 for no restriction on BEV DSM
@ -72,7 +84,7 @@ boilers,--,"{true, false}",Add option for transforming gas into heat using gas b
 resistive_heaters,--,"{true, false}",Add option for transforming electricity into heat using resistive heaters (independently from gas boilers)
 oil_boilers,--,"{true, false}",Add option for transforming oil into heat using boilers
 biomass_boiler,--,"{true, false}",Add option for transforming biomass into heat using boilers
-overdimension_individual_heating,--,"float",Add option for overdimensioning individual heating systems by a certain factor. This allows them to cover heat demand peaks e.g. 10% higher than those in the data with a setting of 1.1.
+overdimension_individual_heating,--,float,Add option for overdimensioning individual heating systems by a certain factor. This allows them to cover heat demand peaks e.g. 10% higher than those in the data with a setting of 1.1.
 chp,--,"{true, false}",Add option for using Combined Heat and Power (CHP)
 micro_chp,--,"{true, false}",Add option for using Combined Heat and Power (CHP) for decentral areas.
 solar_thermal,--,"{true, false}",Add option for using solar thermal to generate heat.
@ -139,6 +151,7 @@ biogas_upgrading_cc,--,"{true, false}",Add option to capture CO2 from biomass up
 conventional_generation,,,Add a more detailed description of conventional carriers. Any power generation requires the consumption of fuel from nodes representing that fuel.
 biomass_to_liquid,--,"{true, false}",Add option for transforming solid biomass into liquid fuel with the same properties as oil
 biosng,--,"{true, false}",Add option for transforming solid biomass into synthesis gas with the same properties as natural gas
 municipal_solid_waste,--,"{true, false}",Add option for municipal solid waste
 limit_max_growth,,,
 -- enable,--,"{true, false}",Add option to limit the maximum growth of a carrier
 -- factor,p.u.,float,The maximum growth factor of a carrier (e.g. 1.3 allows  30% larger than max historic growth)
--- a/doc/preparation.rst
+++ b/doc/preparation.rst
@ -41,11 +41,6 @@ Rule ``build_cutout``
 .. automodule:: build_cutout
 Rule ``prepare_links_p_nom``
 ===============================
 .. automodule:: prepare_links_p_nom
 .. _base:
 Rule ``base_network``
--- a/doc/release_notes.rst
+++ b/doc/release_notes.rst
@ -10,9 +10,19 @@ Release Notes
 Upcoming Release
 ================
 * Added unsustainable biomass potentials for solid, gaseous, and liquid biomass. The potentials can be phased-out and/or
  substituted by the phase-in of sustainable biomass types using the config parameters
  ``biomass: share_unsustainable_use_retained`` and ``biomass: share_sustainable_potential_available``.
 * The rule ``prepare_links_p_nom`` was removed since it was outdated and not used.
 * Changed heat pump COP approximation for central heating to be based on `Jensen et al. (2018) <https://backend.orbit.dtu.dk/ws/portalfiles/portal/151965635/MAIN_Final.pdf>`__ and a default forward temperature of 90C. This is more realistic for district heating than the previously used approximation method.
 * split solid biomass potentials into solid biomass and municipal solid waste. Add option to use municipal solid waste. This option is only activated in combination with the flag ``waste_to_energy``
 * Add option to import solid biomass
-* Add option to produce electrobiofuels (flag ``electrobiofuels`) from solid biomass and hydrogen, as a combination of BtL and Fischer-Tropsch to make more use of the biogenic carbon
+* Add option to produce electrobiofuels (flag ``electrobiofuels``) from solid biomass and hydrogen, as a combination of BtL and Fischer-Tropsch to make more use of the biogenic carbon
 * Add flag ``sector: fossil_fuels`` in config to remove the option of importing fossil fuels
--- a/doc/requirements.txt
+++ b/doc/requirements.txt
@ -17,7 +17,6 @@ tabula-py
 # cartopy
 scikit-learn
 pycountry
 pyyaml
 seaborn
 memory_profiler
--- a/envs/environment.yaml
+++ b/envs/environment.yaml
@ -18,7 +18,6 @@ dependencies:
  # Dependencies of the workflow itself
 - xlrd
 - openpyxl!=3.1.1
 - pycountry
 - seaborn
 - snakemake-minimal>=8.14
 - memory_profiler
--- a/rules/build_electricity.smk
+++ b/rules/build_electricity.smk
@ -2,21 +2,6 @@
 #
 # SPDX-License-Identifier: MIT
 if config["enable"].get("prepare_links_p_nom", False):
    rule prepare_links_p_nom:
        output:
            "data/links_p_nom.csv",
        log:
            logs("prepare_links_p_nom.log"),
        threads: 1
        resources:
            mem_mb=1500,
        conda:
            "../envs/environment.yaml"
        script:
            "../scripts/prepare_links_p_nom.py"
 rule build_electricity_demand:
    params:
@ -106,8 +91,8 @@ rule build_shapes:
    params:
        countries=config_provider("countries"),
    input:
-        naturalearth=ancient("data/bundle/naturalearth/ne_10m_admin_0_countries.shp"),
+        naturalearth=ancient("data/naturalearth/ne_10m_admin_0_countries_deu.shp"),
-        eez=ancient("data/bundle/eez/World_EEZ_v8_2014.shp"),
+        eez=ancient("data/eez/World_EEZ_v12_20231025_gpkg/eez_v12.gpkg"),
        nuts3=ancient("data/bundle/NUTS_2013_60M_SH/data/NUTS_RG_60M_2013.shp"),
        nuts3pop=ancient("data/bundle/nama_10r_3popgdp.tsv.gz"),
        nuts3gdp=ancient("data/bundle/nama_10r_3gdp.tsv.gz"),
--- a/rules/build_sector.smk
+++ b/rules/build_sector.smk
@ -217,13 +217,27 @@ rule build_temperature_profiles:
 rule build_cop_profiles:
    params:
-        heat_pump_sink_T=config_provider("sector", "heat_pump_sink_T"),
+        heat_pump_sink_T_decentral_heating=config_provider(
            "sector", "heat_pump_sink_T_individual_heating"
        ),
        forward_temperature_central_heating=config_provider(
            "sector", "district_heating", "forward_temperature"
        ),
        return_temperature_central_heating=config_provider(
            "sector", "district_heating", "return_temperature"
        ),
        heat_source_cooling_central_heating=config_provider(
            "sector", "district_heating", "heat_source_cooling"
        ),
        heat_pump_cop_approximation_central_heating=config_provider(
            "sector", "district_heating", "heat_pump_cop_approximation"
        ),
        heat_pump_sources=config_provider("sector", "heat_pump_sources"),
    input:
        temp_soil_total=resources("temp_soil_total_elec_s{simpl}_{clusters}.nc"),
        temp_air_total=resources("temp_air_total_elec_s{simpl}_{clusters}.nc"),
    output:
-        cop_soil_total=resources("cop_soil_total_elec_s{simpl}_{clusters}.nc"),
+        cop_profiles=resources("cop_profiles_elec_s{simpl}_{clusters}.nc"),
        cop_air_total=resources("cop_air_total_elec_s{simpl}_{clusters}.nc"),
    resources:
        mem_mb=20000,
    log:
@ -233,7 +247,7 @@ rule build_cop_profiles:
    conda:
        "../envs/environment.yaml"
    script:
-        "../scripts/build_cop_profiles.py"
+        "../scripts/build_cop_profiles/run.py"
 def solar_thermal_cutout(wildcards):
@ -331,7 +345,8 @@ rule build_biomass_potentials:
            "https://zenodo.org/records/10356004/files/ENSPRESO_BIOMASS.xlsx",
            keep_local=True,
        ),
-        nuts2="data/bundle/nuts/NUTS_RG_10M_2013_4326_LEVL_2.geojson",  # https://gisco-services.ec.europa.eu/distribution/v2/nuts/download/#nuts21
+        eurostat="data/eurostat/Balances-April2023",
        nuts2="data/bundle/nuts/NUTS_RG_10M_2013_4326_LEVL_2.geojson",
        regions_onshore=resources("regions_onshore_elec_s{simpl}_{clusters}.geojson"),
        nuts3_population=ancient("data/bundle/nama_10r_3popgdp.tsv.gz"),
        swiss_cantons=ancient("data/ch_cantons.csv"),
@ -344,7 +359,7 @@ rule build_biomass_potentials:
        biomass_potentials=resources(
            "biomass_potentials_s{simpl}_{clusters}_{planning_horizons}.csv"
        ),
-    threads: 1
+    threads: 8
    resources:
        mem_mb=1000,
    log:
@ -940,7 +955,10 @@ rule prepare_sector_network:
        countries=config_provider("countries"),
        adjustments=config_provider("adjustments", "sector"),
        emissions_scope=config_provider("energy", "emissions"),
        biomass=config_provider("biomass"),
        RDIR=RDIR,
        heat_pump_sources=config_provider("sector", "heat_pump_sources"),
        heat_systems=config_provider("sector", "heat_systems"),
    input:
        unpack(input_profile_offwind),
        **rules.cluster_gas_network.output,
@ -1020,8 +1038,7 @@ rule prepare_sector_network:
        ),
        temp_soil_total=resources("temp_soil_total_elec_s{simpl}_{clusters}.nc"),
        temp_air_total=resources("temp_air_total_elec_s{simpl}_{clusters}.nc"),
-        cop_soil_total=resources("cop_soil_total_elec_s{simpl}_{clusters}.nc"),
+        cop_profiles=resources("cop_profiles_elec_s{simpl}_{clusters}.nc"),
        cop_air_total=resources("cop_air_total_elec_s{simpl}_{clusters}.nc"),
        solar_thermal_total=lambda w: (
            resources("solar_thermal_total_elec_s{simpl}_{clusters}.nc")
            if config_provider("sector", "solar_thermal")(w)
--- a/rules/retrieve.smk
+++ b/rules/retrieve.smk
@ -4,6 +4,7 @@
 import requests
 from datetime import datetime, timedelta
 from shutil import move, unpack_archive
 if config["enable"].get("retrieve", "auto") == "auto":
    config["enable"]["retrieve"] = has_internet_access()
@ -15,8 +16,6 @@ if config["enable"]["retrieve"] is False:
 if config["enable"]["retrieve"] and config["enable"].get("retrieve_databundle", True):
    datafiles = [
        "je-e-21.03.02.xls",
        "eez/World_EEZ_v8_2014.shp",
        "naturalearth/ne_10m_admin_0_countries.shp",
        "NUTS_2013_60M_SH/data/NUTS_RG_60M_2013.shp",
        "nama_10r_3popgdp.tsv.gz",
        "nama_10r_3gdp.tsv.gz",
@ -53,6 +52,8 @@ if config["enable"]["retrieve"] and config["enable"].get("retrieve_databundle",
        log:
            "logs/retrieve_eurostat_data.log",
        retries: 2
        conda:
            "../envs/retrieve.yaml"
        script:
            "../scripts/retrieve_eurostat_data.py"
@ -62,6 +63,8 @@ if config["enable"]["retrieve"] and config["enable"].get("retrieve_databundle",
        log:
            "logs/retrieve_eurostat_household_data.log",
        retries: 2
        conda:
            "../envs/retrieve.yaml"
        script:
            "../scripts/retrieve_eurostat_household_data.py"
@ -211,6 +214,64 @@ if config["enable"]["retrieve"]:
            move(input[0], output[0])
 if config["enable"]["retrieve"]:
    rule retrieve_eez:
        params:
            zip="data/eez/World_EEZ_v12_20231025_gpkg.zip",
        output:
            gpkg="data/eez/World_EEZ_v12_20231025_gpkg/eez_v12.gpkg",
        run:
            import os
            import requests
            from uuid import uuid4
            name = str(uuid4())[:8]
            org = str(uuid4())[:8]
            response = requests.post(
                "https://www.marineregions.org/download_file.php",
                params={"name": "World_EEZ_v12_20231025_gpkg.zip"},
                data={
                    "name": name,
                    "organisation": org,
                    "email": f"{name}@{org}.org",
                    "country": "Germany",
                    "user_category": "academia",
                    "purpose_category": "Research",
                    "agree": "1",
                },
            )
            with open(params["zip"], "wb") as f:
                f.write(response.content)
            output_folder = Path(params["zip"]).parent
            unpack_archive(params["zip"], output_folder)
            os.remove(params["zip"])
 if config["enable"]["retrieve"]:
    # Download directly from naciscdn.org which is a redirect from naturalearth.com
    # (https://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-admin-0-countries/)
    # Use point-of-view (POV) variant of Germany so that Crimea is included.
    rule retrieve_naturalearth_countries:
        input:
            storage(
                "https://naciscdn.org/naturalearth/10m/cultural/ne_10m_admin_0_countries_deu.zip"
            ),
        params:
            zip="data/naturalearth/ne_10m_admin_0_countries_deu.zip",
        output:
            countries="data/naturalearth/ne_10m_admin_0_countries_deu.shp",
        run:
            move(input[0], params["zip"])
            output_folder = Path(output["countries"]).parent
            unpack_archive(params["zip"], output_folder)
            os.remove(params["zip"])
 if config["enable"]["retrieve"]:
    # Some logic to find the correct file URL
    # Sometimes files are released delayed or ahead of schedule, check which file is currently available
--- a/rules/solve_myopic.smk
+++ b/rules/solve_myopic.smk
@ -9,6 +9,7 @@ rule add_existing_baseyear:
        sector=config_provider("sector"),
        existing_capacities=config_provider("existing_capacities"),
        costs=config_provider("costs"),
        heat_pump_sources=config_provider("sector", "heat_pump_sources"),
    input:
        network=RESULTS
        + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
@ -21,8 +22,7 @@ rule add_existing_baseyear:
                config_provider("scenario", "planning_horizons", 0)(w)
            )
        ),
-        cop_soil_total=resources("cop_soil_total_elec_s{simpl}_{clusters}.nc"),
+        cop_profiles=resources("cop_profiles_elec_s{simpl}_{clusters}.nc"),
        cop_air_total=resources("cop_air_total_elec_s{simpl}_{clusters}.nc"),
        existing_heating_distribution=resources(
            "existing_heating_distribution_elec_s{simpl}_{clusters}_{planning_horizons}.csv"
        ),
@ -69,6 +69,7 @@ rule add_brownfield:
        snapshots=config_provider("snapshots"),
        drop_leap_day=config_provider("enable", "drop_leap_day"),
        carriers=config_provider("electricity", "renewable_carriers"),
        heat_pump_sources=config_provider("sector", "heat_pump_sources"),
    input:
        unpack(input_profile_tech_brownfield),
        simplify_busmap=resources("busmap_elec_s{simpl}.csv"),
@ -77,8 +78,7 @@ rule add_brownfield:
        + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
        network_p=solved_previous_horizon,  #solved network at previous time step
        costs=resources("costs_{planning_horizons}.csv"),
-        cop_soil_total=resources("cop_soil_total_elec_s{simpl}_{clusters}.nc"),
+        cop_profiles=resources("cop_profiles_elec_s{simpl}_{clusters}.nc"),
        cop_air_total=resources("cop_air_total_elec_s{simpl}_{clusters}.nc"),
    output:
        RESULTS
        + "prenetworks-brownfield/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
--- a/rules/solve_perfect.smk
+++ b/rules/solve_perfect.smk
@ -7,6 +7,7 @@ rule add_existing_baseyear:
        sector=config_provider("sector"),
        existing_capacities=config_provider("existing_capacities"),
        costs=config_provider("costs"),
        heat_pump_sources=config_provider("sector", "heat_pump_sources"),
    input:
        network=RESULTS
        + "prenetworks/elec_s{simpl}_{clusters}_l{ll}_{opts}_{sector_opts}_{planning_horizons}.nc",
@ -19,8 +20,7 @@ rule add_existing_baseyear:
                config_provider("scenario", "planning_horizons", 0)(w)
            )
        ),
-        cop_soil_total=resources("cop_soil_total_elec_s{simpl}_{clusters}.nc"),
+        cop_profiles=resources("cop_profiles_elec_s{simpl}_{clusters}.nc"),
        cop_air_total=resources("cop_air_total_elec_s{simpl}_{clusters}.nc"),
        existing_heating_distribution=resources(
            "existing_heating_distribution_elec_s{simpl}_{clusters}_{planning_horizons}.csv"
        ),
--- a/scripts/add_existing_baseyear.py
+++ b/scripts/add_existing_baseyear.py
@ -24,6 +24,10 @@ from _helpers import (
 from add_electricity import sanitize_carriers
 from prepare_sector_network import cluster_heat_buses, define_spatial, prepare_costs
 from scripts.definitions.heat_sector import HeatSector
 from scripts.definitions.heat_system import HeatSystem
 from scripts.definitions.heat_system_type import HeatSystemType
 logger = logging.getLogger(__name__)
 cc = coco.CountryConverter()
 idx = pd.IndexSlice
@ -416,14 +420,14 @@ def add_power_capacities_installed_before_baseyear(n, grouping_years, costs, bas
 def add_heating_capacities_installed_before_baseyear(
-    n,
+    n: pypsa.Network,
-    baseyear,
+    baseyear: int,
-    grouping_years,
+    grouping_years: list,
-    ashp_cop,
+    cop: dict,
-    gshp_cop,
+    time_dep_hp_cop: bool,
-    time_dep_hp_cop,
+    costs: pd.DataFrame,
-    costs,
+    default_lifetime: int,
-    default_lifetime,
+    existing_heating: pd.DataFrame,
 ):
    """
    Parameters
@ -435,141 +439,158 @@ def add_heating_capacities_installed_before_baseyear(
        currently assumed heating capacities split between residential and
        services proportional to heating load in both 50% capacities
        in rural buses 50% in urban buses
    cop: xr.DataArray
        DataArray with time-dependent coefficients of performance (COPs) heat pumps. Coordinates are heat sources (see config), heat system types (see :file:`scripts/enums/HeatSystemType.py`), nodes and snapshots.
    time_dep_hp_cop: bool
        If True, time-dependent (dynamic) COPs are used for heat pumps
    """
    logger.debug(f"Adding heating capacities installed before {baseyear}")
-    existing_heating = pd.read_csv(
+    for heat_system in existing_heating.columns.get_level_values(0).unique():
-        snakemake.input.existing_heating_distribution, header=[0, 1], index_col=0
+        heat_system = HeatSystem(heat_system)
    )
-    for name in existing_heating.columns.get_level_values(0).unique():
+        nodes = pd.Index(
-        name_type = "central" if name == "urban central" else "decentral"
+            n.buses.location[n.buses.index.str.contains(f"{heat_system} heat")]
        )
-        nodes = pd.Index(n.buses.location[n.buses.index.str.contains(f"{name} heat")])
+        if (not heat_system == HeatSystem.URBAN_CENTRAL) and options[
-
+            "electricity_distribution_grid"
-        if (name_type != "central") and options["electricity_distribution_grid"]:
+        ]:
            nodes_elec = nodes + " low voltage"
        else:
            nodes_elec = nodes
-        heat_pump_type = "air" if "urban" in name else "ground"
+            too_large_grouping_years = [
-
+                gy for gy in grouping_years if gy >= int(baseyear)
        # Add heat pumps
        costs_name = f"decentral {heat_pump_type}-sourced heat pump"
        cop = {"air": ashp_cop, "ground": gshp_cop}
        if time_dep_hp_cop:
            efficiency = cop[heat_pump_type][nodes]
        else:
            efficiency = costs.at[costs_name, "efficiency"]
        too_large_grouping_years = [gy for gy in grouping_years if gy >= int(baseyear)]
        if too_large_grouping_years:
            logger.warning(
                f"Grouping years >= baseyear are ignored. Dropping {too_large_grouping_years}."
            )
        valid_grouping_years = pd.Series(
            [
                int(grouping_year)
                for grouping_year in grouping_years
                if int(grouping_year) + default_lifetime > int(baseyear)
                and int(grouping_year) < int(baseyear)
            ]
-        )
+            if too_large_grouping_years:
                logger.warning(
                    f"Grouping years >= baseyear are ignored. Dropping {too_large_grouping_years}."
                )
            valid_grouping_years = pd.Series(
                [
                    int(grouping_year)
                    for grouping_year in grouping_years
                    if int(grouping_year) + default_lifetime > int(baseyear)
                    and int(grouping_year) < int(baseyear)
                ]
            )
-        assert valid_grouping_years.is_monotonic_increasing
+            assert valid_grouping_years.is_monotonic_increasing
-        # get number of years of each interval
+            # get number of years of each interval
-        _years = valid_grouping_years.diff()
+            _years = valid_grouping_years.diff()
-        # Fill NA from .diff() with value for the first interval
+            # Fill NA from .diff() with value for the first interval
-        _years[0] = valid_grouping_years[0] - baseyear + default_lifetime
+            _years[0] = valid_grouping_years[0] - baseyear + default_lifetime
-        # Installation is assumed to be linear for the past
+            # Installation is assumed to be linear for the past
-        ratios = _years / _years.sum()
+            ratios = _years / _years.sum()
        for ratio, grouping_year in zip(ratios, valid_grouping_years):
            # Add heat pumps
            for heat_source in snakemake.params.heat_pump_sources[
                heat_system.system_type.value
            ]:
                costs_name = heat_system.heat_pump_costs_name(heat_source)
-            n.madd(
+                efficiency = (
-                "Link",
+                    cop.sel(
-                nodes,
+                        heat_system=heat_system.system_type.value,
-                suffix=f" {name} {heat_pump_type} heat pump-{grouping_year}",
+                        heat_source=heat_source,
-                bus0=nodes_elec,
+                        name=nodes,
-                bus1=nodes + " " + name + " heat",
+                    )
-                carrier=f"{name} {heat_pump_type} heat pump",
+                    .to_pandas()
-                efficiency=efficiency,
+                    .reindex(index=n.snapshots)
-                capital_cost=costs.at[costs_name, "efficiency"]
+                    if time_dep_hp_cop
-                * costs.at[costs_name, "fixed"],
+                    else costs.at[costs_name, "efficiency"]
-                p_nom=existing_heating.loc[nodes, (name, f"{heat_pump_type} heat pump")]
+                )
-                * ratio
+
-                / costs.at[costs_name, "efficiency"],
+                n.madd(
-                build_year=int(grouping_year),
+                    "Link",
-                lifetime=costs.at[costs_name, "lifetime"],
+                    nodes,
-            )
+                    suffix=f" {heat_system} {heat_source} heat pump-{grouping_year}",
                    bus0=nodes_elec,
                    bus1=nodes + " " + heat_system.value + " heat",
                    carrier=f"{heat_system} {heat_source} heat pump",
                    efficiency=efficiency,
                    capital_cost=costs.at[costs_name, "efficiency"]
                    * costs.at[costs_name, "fixed"],
                    p_nom=existing_heating.loc[
                        nodes, (heat_system.value, f"{heat_source} heat pump")
                    ]
                    * 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
            n.madd(
                "Link",
                nodes,
-                suffix=f" {name} resistive heater-{grouping_year}",
+                suffix=f" {heat_system} resistive heater-{grouping_year}",
                bus0=nodes_elec,
-                bus1=nodes + " " + name + " heat",
+                bus1=nodes + " " + heat_system.value + " heat",
-                carrier=name + " resistive heater",
+                carrier=heat_system.value + " resistive heater",
-                efficiency=costs.at[f"{name_type} resistive heater", "efficiency"],
+                efficiency=costs.at[
                    heat_system.resistive_heater_costs_name, "efficiency"
                ],
                capital_cost=(
-                    costs.at[f"{name_type} resistive heater", "efficiency"]
+                    costs.at[heat_system.resistive_heater_costs_name, "efficiency"]
-                    * costs.at[f"{name_type} resistive heater", "fixed"]
+                    * costs.at[heat_system.resistive_heater_costs_name, "fixed"]
                ),
                p_nom=(
-                    existing_heating.loc[nodes, (name, "resistive heater")]
+                    existing_heating.loc[nodes, (heat_system.value, "resistive heater")]
                    * ratio
-                    / costs.at[f"{name_type} resistive heater", "efficiency"]
+                    / costs.at[heat_system.resistive_heater_costs_name, "efficiency"]
                ),
                build_year=int(grouping_year),
-                lifetime=costs.at[f"{name_type} resistive heater", "lifetime"],
+                lifetime=costs.at[heat_system.resistive_heater_costs_name, "lifetime"],
            )
            n.madd(
                "Link",
                nodes,
-                suffix=f" {name} gas boiler-{grouping_year}",
+                suffix=f"{heat_system} gas boiler-{grouping_year}",
                bus0="EU gas" if "EU gas" in spatial.gas.nodes else nodes + " gas",
-                bus1=nodes + " " + name + " heat",
+                bus1=nodes + " " + heat_system.value + " heat",
                bus2="co2 atmosphere",
-                carrier=name + " gas boiler",
+                carrier=heat_system.value + " gas boiler",
-                efficiency=costs.at[f"{name_type} gas boiler", "efficiency"],
+                efficiency=costs.at[heat_system.gas_boiler_costs_name, "efficiency"],
                efficiency2=costs.at["gas", "CO2 intensity"],
                capital_cost=(
-                    costs.at[f"{name_type} gas boiler", "efficiency"]
+                    costs.at[heat_system.gas_boiler_costs_name, "efficiency"]
-                    * costs.at[f"{name_type} gas boiler", "fixed"]
+                    * costs.at[heat_system.gas_boiler_costs_name, "fixed"]
                ),
                p_nom=(
-                    existing_heating.loc[nodes, (name, "gas boiler")]
+                    existing_heating.loc[nodes, (heat_system.value, "gas boiler")]
                    * ratio
-                    / costs.at[f"{name_type} gas boiler", "efficiency"]
+                    / costs.at[heat_system.gas_boiler_costs_name, "efficiency"]
                ),
                build_year=int(grouping_year),
-                lifetime=costs.at[f"{name_type} gas boiler", "lifetime"],
+                lifetime=costs.at[heat_system.gas_boiler_costs_name, "lifetime"],
            )
            n.madd(
                "Link",
                nodes,
-                suffix=f" {name} oil boiler-{grouping_year}",
+                suffix=f" {heat_system} oil boiler-{grouping_year}",
                bus0=spatial.oil.nodes,
-                bus1=nodes + " " + name + " heat",
+                bus1=nodes + " " + heat_system.value + " heat",
                bus2="co2 atmosphere",
-                carrier=name + " oil boiler",
+                carrier=heat_system.value + " oil boiler",
-                efficiency=costs.at["decentral oil boiler", "efficiency"],
+                efficiency=costs.at[heat_system.oil_boiler_costs_name, "efficiency"],
                efficiency2=costs.at["oil", "CO2 intensity"],
-                capital_cost=costs.at["decentral oil boiler", "efficiency"]
+                capital_cost=costs.at[heat_system.oil_boiler_costs_name, "efficiency"]
-                * costs.at["decentral oil boiler", "fixed"],
+                * costs.at[heat_system.oil_boiler_costs_name, "fixed"],
                p_nom=(
-                    existing_heating.loc[nodes, (name, "oil boiler")]
+                    existing_heating.loc[nodes, (heat_system.value, "oil boiler")]
                    * ratio
-                    / costs.at["decentral oil boiler", "efficiency"]
+                    / costs.at[heat_system.oil_boiler_costs_name, "efficiency"]
                ),
                build_year=int(grouping_year),
-                lifetime=costs.at[f"{name_type} gas boiler", "lifetime"],
+                lifetime=costs.at[
                    f"{heat_system.central_or_decentral} gas boiler", "lifetime"
                ],
            )
            # delete links with p_nom=nan corresponding to extra nodes in country
@ -639,29 +660,22 @@ if __name__ == "__main__":
    )
    if options["heating"]:
-        time_dep_hp_cop = options["time_dep_hp_cop"]
+
        ashp_cop = (
            xr.open_dataarray(snakemake.input.cop_air_total)
            .to_pandas()
            .reindex(index=n.snapshots)
        )
        gshp_cop = (
            xr.open_dataarray(snakemake.input.cop_soil_total)
            .to_pandas()
            .reindex(index=n.snapshots)
        )
        default_lifetime = snakemake.params.existing_capacities[
            "default_heating_lifetime"
        ]
        add_heating_capacities_installed_before_baseyear(
-            n,
+            n=n,
-            baseyear,
+            baseyear=baseyear,
-            grouping_years_heat,
+            grouping_years=grouping_years_heat,
-            ashp_cop,
+            cop=xr.open_dataarray(snakemake.input.cop_profiles),
-            gshp_cop,
+            time_dep_hp_cop=options["time_dep_hp_cop"],
-            time_dep_hp_cop,
+            costs=costs,
-            costs,
+            default_lifetime=snakemake.params.existing_capacities[
-            default_lifetime,
+                "default_heating_lifetime"
            ],
            existing_heating=pd.read_csv(
                snakemake.input.existing_heating_distribution,
                header=[0, 1],
                index_col=0,
            ),
        )
    if options.get("cluster_heat_buses", False):
--- a/scripts/build_biomass_potentials.py
+++ b/scripts/build_biomass_potentials.py
@ -13,11 +13,51 @@ import geopandas as gpd
 import numpy as np
 import pandas as pd
 from _helpers import configure_logging, set_scenario_config
 from build_energy_totals import build_eurostat
 logger = logging.getLogger(__name__)
 AVAILABLE_BIOMASS_YEARS = [2010, 2020, 2030, 2040, 2050]
 def _calc_unsustainable_potential(df, df_unsustainable, share_unsus, resource_type):
    """
    Calculate the unsustainable biomass potential for a given resource type or
    regex.
    Parameters
    ----------
    df : pd.DataFrame
        The dataframe with sustainable biomass potentials.
    df_unsustainable : pd.DataFrame
        The dataframe with unsustainable biomass potentials.
    share_unsus : float
        The share of unsustainable biomass potential retained.
    resource_type : str or regex
        The resource type to calculate the unsustainable potential for.
    Returns
    -------
    pd.Series
        The unsustainable biomass potential for the given resource type or regex.
    """
    if "|" in resource_type:
        resource_potential = df_unsustainable.filter(regex=resource_type).sum(axis=1)
    else:
        resource_potential = df_unsustainable[resource_type]
    return (
        df.apply(
            lambda c: c.sum()
            / df.loc[df.index.str[:2] == c.name[:2]].sum().sum()
            * resource_potential.loc[c.name[:2]],
            axis=1,
        )
        .mul(share_unsus)
        .clip(lower=0)
    )
 def build_nuts_population_data(year=2013):
    pop = pd.read_csv(
        snakemake.input.nuts3_population,
@ -211,15 +251,104 @@ def convert_nuts2_to_regions(bio_nuts2, regions):
    return bio_regions
 def add_unsustainable_potentials(df):
    """
    Add unsustainable biomass potentials to the given dataframe. The difference
    between the data of JRC and Eurostat is assumed to be unsustainable
    biomass.
    Parameters
    ----------
    df : pd.DataFrame
        The dataframe with sustainable biomass potentials.
    unsustainable_biomass : str
        Path to the file with unsustainable biomass potentials.
    Returns
    -------
    pd.DataFrame
        The dataframe with added unsustainable biomass potentials.
    """
    if "GB" in snakemake.config["countries"]:
        latest_year = 2019
    else:
        latest_year = 2021
    idees_rename = {"GR": "EL", "GB": "UK"}
    df_unsustainable = (
        build_eurostat(
            countries=snakemake.config["countries"],
            input_eurostat=snakemake.input.eurostat,
            nprocesses=int(snakemake.threads),
        )
        .xs(
            max(min(latest_year, int(snakemake.wildcards.planning_horizons)), 1990),
            level=1,
        )
        .xs("Primary production", level=2)
        .droplevel([1, 2, 3])
    )
    df_unsustainable.index = df_unsustainable.index.str.strip()
    df_unsustainable = df_unsustainable.rename(
        {v: k for k, v in idees_rename.items()}, axis=0
    )
    bio_carriers = [
        "Primary solid biofuels",
        "Biogases",
        "Renewable municipal waste",
        "Pure biogasoline",
        "Blended biogasoline",
        "Pure biodiesels",
        "Blended biodiesels",
        "Pure bio jet kerosene",
        "Blended bio jet kerosene",
        "Other liquid biofuels",
    ]
    df_unsustainable = df_unsustainable[bio_carriers]
    # Phase out unsustainable biomass potentials linearly from 2020 to 2035 while phasing in sustainable potentials
    share_unsus = params.get("share_unsustainable_use_retained").get(investment_year)
    df_wo_ch = df.drop(df.filter(regex="CH\d", axis=0).index)
    # Calculate unsustainable solid biomass
    df_wo_ch["unsustainable solid biomass"] = _calc_unsustainable_potential(
        df_wo_ch, df_unsustainable, share_unsus, "Primary solid biofuels"
    )
    # Calculate unsustainable biogas
    df_wo_ch["unsustainable biogas"] = _calc_unsustainable_potential(
        df_wo_ch, df_unsustainable, share_unsus, "Biogases"
    )
    # Calculate unsustainable bioliquids
    df_wo_ch["unsustainable bioliquids"] = _calc_unsustainable_potential(
        df_wo_ch,
        df_unsustainable,
        share_unsus,
        resource_type="gasoline|diesel|kerosene|liquid",
    )
    share_sus = params.get("share_sustainable_potential_available").get(investment_year)
    df *= share_sus
    df = df.join(df_wo_ch.filter(like="unsustainable")).fillna(0)
    return df
 if __name__ == "__main__":
    if "snakemake" not in globals():
        from _helpers import mock_snakemake
        snakemake = mock_snakemake(
            "build_biomass_potentials",
            simpl="",
-            clusters="5",
+            clusters="37",
-            planning_horizons=2050,
+            planning_horizons=2020,
        )
    configure_logging(snakemake)
@ -269,6 +398,8 @@ if __name__ == "__main__":
    grouper = {v: k for k, vv in params["classes"].items() for v in vv}
    df = df.T.groupby(grouper).sum().T
    df = add_unsustainable_potentials(df)
    df *= 1e6  # TWh/a to MWh/a
    df.index.name = "MWh/a"
--- a/scripts/build_cop_profiles.py
+++ b/scripts/build_cop_profiles.py
@ -1,69 +0,0 @@
 # -*- coding: utf-8 -*-
 # SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
 #
 # SPDX-License-Identifier: MIT
 """
 Build coefficient of performance (COP) time series for air- or ground-sourced
 heat pumps.
 The COP is approximated as a quatratic function of the temperature difference between source and
 sink, based on Staffell et al. 2012.
 This rule is executed in ``build_sector.smk``.
 Relevant Settings
 -----------------
 .. code:: yaml
    heat_pump_sink_T:
 Inputs:
 -------
 - ``resources/<run_name>/temp_soil_total_elec_s<simpl>_<clusters>.nc``: Soil temperature (total) time series.
 - ``resources/<run_name>/temp_air_total_elec_s<simpl>_<clusters>.nc``: Ambient air temperature (total) time series.
 Outputs:
 --------
 - ``resources/cop_soil_total_elec_s<simpl>_<clusters>.nc``: COP (ground-sourced) time series (total).
 - ``resources/cop_air_total_elec_s<simpl>_<clusters>.nc``: COP (air-sourced) time series (total).
 References
 ----------
 [1] Staffell et al., Energy & Environmental Science 11 (2012): A review of domestic heat pumps, https://doi.org/10.1039/C2EE22653G.
 """
 import xarray as xr
 from _helpers import set_scenario_config
 def coefficient_of_performance(delta_T, source="air"):
    if source == "air":
        return 6.81 - 0.121 * delta_T + 0.000630 * delta_T**2
    elif source == "soil":
        return 8.77 - 0.150 * delta_T + 0.000734 * delta_T**2
    else:
        raise NotImplementedError("'source' must be one of  ['air', 'soil']")
 if __name__ == "__main__":
    if "snakemake" not in globals():
        from _helpers import mock_snakemake
        snakemake = mock_snakemake(
            "build_cop_profiles",
            simpl="",
            clusters=48,
        )
    set_scenario_config(snakemake)
    for source in ["air", "soil"]:
        source_T = xr.open_dataarray(snakemake.input[f"temp_{source}_total"])
        delta_T = snakemake.params.heat_pump_sink_T - source_T
        cop = coefficient_of_performance(delta_T, source)
        cop.to_netcdf(snakemake.output[f"cop_{source}_total"])
--- a/scripts/build_cop_profiles/BaseCopApproximator.py
+++ b/scripts/build_cop_profiles/BaseCopApproximator.py
@ -0,0 +1,111 @@
 # -*- coding: utf-8 -*-
 # SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
 #
 # SPDX-License-Identifier: MIT
 from abc import ABC, abstractmethod
 from typing import Union
 import numpy as np
 import xarray as xr
 class BaseCopApproximator(ABC):
    """
    Abstract class for approximating the coefficient of performance (COP) of a
    heat pump.
    Attributes:
    ----------
    forward_temperature_celsius : Union[xr.DataArray, np.array]
        The forward temperature in Celsius.
    source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
        The source inlet temperature in Celsius.
    Methods:
    -------
    __init__(self, forward_temperature_celsius, source_inlet_temperature_celsius)
        Initialize CopApproximator.
    approximate_cop(self)
        Approximate heat pump coefficient of performance (COP).
    celsius_to_kelvin(t_celsius)
        Convert temperature from Celsius to Kelvin.
    logarithmic_mean(t_hot, t_cold)
        Calculate the logarithmic mean temperature difference.
    """
    def __init__(
        self,
        forward_temperature_celsius: Union[xr.DataArray, np.array],
        source_inlet_temperature_celsius: Union[xr.DataArray, np.array],
    ):
        """
        Initialize CopApproximator.
        Parameters:
        ----------
        forward_temperature_celsius : Union[xr.DataArray, np.array]
            The forward temperature in Celsius.
        source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
            The source inlet temperature in Celsius.
        """
        pass
    @abstractmethod
    def approximate_cop(self) -> Union[xr.DataArray, np.array]:
        """
        Approximate heat pump coefficient of performance (COP).
        Returns:
        -------
        Union[xr.DataArray, np.array]
            The calculated COP values.
        """
        pass
    @staticmethod
    def celsius_to_kelvin(
        t_celsius: Union[float, xr.DataArray, np.array]
    ) -> Union[float, xr.DataArray, np.array]:
        """
        Convert temperature from Celsius to Kelvin.
        Parameters:
        ----------
        t_celsius : Union[float, xr.DataArray, np.array]
            Temperature in Celsius.
        Returns:
        -------
        Union[float, xr.DataArray, np.array]
            Temperature in Kelvin.
        """
        if (np.asarray(t_celsius) > 200).any():
            raise ValueError(
                "t_celsius > 200. Are you sure you are using the right units?"
            )
        return t_celsius + 273.15
    @staticmethod
    def logarithmic_mean(
        t_hot: Union[float, xr.DataArray, np.ndarray],
        t_cold: Union[float, xr.DataArray, np.ndarray],
    ) -> Union[float, xr.DataArray, np.ndarray]:
        """
        Calculate the logarithmic mean temperature difference.
        Parameters:
        ----------
        t_hot : Union[float, xr.DataArray, np.ndarray]
            Hot temperature.
        t_cold : Union[float, xr.DataArray, np.ndarray]
            Cold temperature.
        Returns:
        -------
        Union[float, xr.DataArray, np.ndarray]
            Logarithmic mean temperature difference.
        """
        if (np.asarray(t_hot <= t_cold)).any():
            raise ValueError("t_hot must be greater than t_cold")
        return (t_hot - t_cold) / np.log(t_hot / t_cold)
--- a/scripts/build_cop_profiles/CentralHeatingCopApproximator.py
+++ b/scripts/build_cop_profiles/CentralHeatingCopApproximator.py
@ -0,0 +1,392 @@
 # -*- coding: utf-8 -*-
 # SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
 #
 # SPDX-License-Identifier: MIT
 from typing import Union
 import numpy as np
 import xarray as xr
 from BaseCopApproximator import BaseCopApproximator
 class CentralHeatingCopApproximator(BaseCopApproximator):
    """
    Approximate the coefficient of performance (COP) for a heat pump in a
    central heating system (district heating).
    Uses an approximation method proposed by Jensen et al. (2018) and
    default parameters from Pieper et al. (2020). The method is based on
    a thermodynamic heat pump model with some hard-to-know parameters
    being approximated.
    Attributes:
    ----------
    forward_temperature_celsius : Union[xr.DataArray, np.array]
        The forward temperature in Celsius.
    return_temperature_celsius : Union[xr.DataArray, np.array]
        The return temperature in Celsius.
    source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
        The source inlet temperature in Celsius.
    source_outlet_temperature_celsius : Union[xr.DataArray, np.array]
        The source outlet temperature in Celsius.
    delta_t_pinch_point : float, optional
        The pinch point temperature difference, by default 5.
    isentropic_compressor_efficiency : float, optional
        The isentropic compressor efficiency, by default 0.8.
    heat_loss : float, optional
        The heat loss, by default 0.0.
    Methods:
    -------
    __init__(
        forward_temperature_celsius: Union[xr.DataArray, np.array],
        source_inlet_temperature_celsius: Union[xr.DataArray, np.array],
        return_temperature_celsius: Union[xr.DataArray, np.array],
        source_outlet_temperature_celsius: Union[xr.DataArray, np.array],
        delta_t_pinch_point: float = 5,
        isentropic_compressor_efficiency: float = 0.8,
        heat_loss: float = 0.0,
    ) -> None:
        Initializes the CentralHeatingCopApproximator object.
    approximate_cop(self) -> Union[xr.DataArray, np.array]:
        Calculate the coefficient of performance (COP) for the system.
    _approximate_delta_t_refrigerant_source(
        self, delta_t_source: Union[xr.DataArray, np.array]
    ) -> Union[xr.DataArray, np.array]:
        Approximates the temperature difference between the refrigerant and the source.
    _approximate_delta_t_refrigerant_sink(
        self,
        refrigerant: str = "ammonia",
        a: float = {"ammonia": 0.2, "isobutane": -0.0011},
        b: float = {"ammonia": 0.2, "isobutane": 0.3},
        c: float = {"ammonia": 0.016, "isobutane": 2.4},
    ) -> Union[xr.DataArray, np.array]:
        Approximates the temperature difference between the refrigerant and heat sink.
    _ratio_evaporation_compression_work_approximation(
        self,
        refrigerant: str = "ammonia",
        a: float = {"ammonia": 0.0014, "isobutane": 0.0035},
    ) -> Union[xr.DataArray, np.array]:
        Calculate the ratio of evaporation to compression work based on approximation.
    _approximate_delta_t_refrigerant_sink(
        self,
        refrigerant: str = "ammonia",
        a: float = {"ammonia": 0.2, "isobutane": -0.0011},
        b: float = {"ammonia": 0.2, "isobutane": 0.3},
        c: float = {"ammonia": 0.016, "isobutane": 2.4},
    ) -> Union[xr.DataArray, np.array]:
        Approximates the temperature difference between the refrigerant and heat sink.
    _ratio_evaporation_compression_work_approximation(
        self,
        refrigerant: str = "ammonia",
        a: float = {"ammonia": 0.0014, "isobutane": 0.0035},
    ) -> Union[xr.DataArray, np.array]:
        Calculate the ratio of evaporation to compression work based on approximation.
    """
    def __init__(
        self,
        forward_temperature_celsius: Union[xr.DataArray, np.array],
        source_inlet_temperature_celsius: Union[xr.DataArray, np.array],
        return_temperature_celsius: Union[xr.DataArray, np.array],
        source_outlet_temperature_celsius: Union[xr.DataArray, np.array],
        delta_t_pinch_point: float = 5,
        isentropic_compressor_efficiency: float = 0.8,
        heat_loss: float = 0.0,
    ) -> None:
        """
        Initializes the CentralHeatingCopApproximator object.
        Parameters:
        ----------
        forward_temperature_celsius : Union[xr.DataArray, np.array]
            The forward temperature in Celsius.
        return_temperature_celsius : Union[xr.DataArray, np.array]
            The return temperature in Celsius.
        source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
            The source inlet temperature in Celsius.
        source_outlet_temperature_celsius : Union[xr.DataArray, np.array]
            The source outlet temperature in Celsius.
        delta_t_pinch_point : float, optional
            The pinch point temperature difference, by default 5.
        isentropic_compressor_efficiency : float, optional
            The isentropic compressor efficiency, by default 0.8.
        heat_loss : float, optional
            The heat loss, by default 0.0.
        """
        self.t_source_in_kelvin = BaseCopApproximator.celsius_to_kelvin(
            source_inlet_temperature_celsius
        )
        self.t_sink_out_kelvin = BaseCopApproximator.celsius_to_kelvin(
            forward_temperature_celsius
        )
        self.t_sink_in_kelvin = BaseCopApproximator.celsius_to_kelvin(
            return_temperature_celsius
        )
        self.t_source_out = BaseCopApproximator.celsius_to_kelvin(
            source_outlet_temperature_celsius
        )
        self.isentropic_efficiency_compressor_kelvin = isentropic_compressor_efficiency
        self.heat_loss = heat_loss
        self.delta_t_pinch = delta_t_pinch_point
    def approximate_cop(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the coefficient of performance (COP) for the system.
        Returns:
        --------
            Union[xr.DataArray, np.array]: The calculated COP values.
        """
        return (
            self.ideal_lorenz_cop
            * (
                (
                    1
                    + (self.delta_t_refrigerant_sink + self.delta_t_pinch)
                    / self.t_sink_mean_kelvin
                )
                / (
                    1
                    + (
                        self.delta_t_refrigerant_sink
                        + self.delta_t_refrigerant_source
                        + 2 * self.delta_t_pinch
                    )
                    / self.delta_t_lift
                )
            )
            * self.isentropic_efficiency_compressor_kelvin
            * (1 - self.ratio_evaporation_compression_work)
            + 1
            - self.isentropic_efficiency_compressor_kelvin
            - self.heat_loss
        )
    @property
    def t_sink_mean_kelvin(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the logarithmic mean temperature difference between the cold
        and hot sinks.
        Returns
        -------
        Union[xr.DataArray, np.array]
            The mean temperature difference.
        """
        return BaseCopApproximator.logarithmic_mean(
            t_cold=self.t_sink_in_kelvin, t_hot=self.t_sink_out_kelvin
        )
    @property
    def t_source_mean_kelvin(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the logarithmic mean temperature of the heat source.
        Returns
        -------
        Union[xr.DataArray, np.array]
            The mean temperature of the heat source.
        """
        return BaseCopApproximator.logarithmic_mean(
            t_hot=self.t_source_in_kelvin, t_cold=self.t_source_out
        )
    @property
    def delta_t_lift(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the temperature lift as the difference between the
        logarithmic sink and source temperatures.
        Returns
        -------
        Union[xr.DataArray, np.array]
            The temperature difference between the sink and source.
        """
        return self.t_sink_mean_kelvin - self.t_source_mean_kelvin
    @property
    def ideal_lorenz_cop(self) -> Union[xr.DataArray, np.array]:
        """
        Ideal Lorenz coefficient of performance (COP).
        The ideal Lorenz COP is calculated as the ratio of the mean sink temperature
        to the lift temperature difference.
        Returns
        -------
        np.array
            The ideal Lorenz COP.
        """
        return self.t_sink_mean_kelvin / self.delta_t_lift
    @property
    def delta_t_refrigerant_source(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the temperature difference between the refrigerant source
        inlet and outlet.
        Returns
        -------
        Union[xr.DataArray, np.array]
            The temperature difference between the refrigerant source inlet and outlet.
        """
        return self._approximate_delta_t_refrigerant_source(
            delta_t_source=self.t_source_in_kelvin - self.t_source_out
        )
    @property
    def delta_t_refrigerant_sink(self) -> Union[xr.DataArray, np.array]:
        """
        Temperature difference between the refrigerant and the sink based on
        approximation.
        Returns
        -------
        Union[xr.DataArray, np.array]
            The temperature difference between the refrigerant and the sink.
        """
        return self._approximate_delta_t_refrigerant_sink()
    @property
    def ratio_evaporation_compression_work(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the ratio of evaporation to compression work based on
        approximation.
        Returns
        -------
        Union[xr.DataArray, np.array]
            The calculated ratio of evaporation to compression work.
        """
        return self._ratio_evaporation_compression_work_approximation()
    @property
    def delta_t_sink(self) -> Union[xr.DataArray, np.array]:
        """
        Calculate the temperature difference at the sink.
        Returns
        -------
        Union[xr.DataArray, np.array]
            The temperature difference at the sink.
        """
        return self.t_sink_out_kelvin - self.t_sink_in_kelvin
    def _approximate_delta_t_refrigerant_source(
        self, delta_t_source: Union[xr.DataArray, np.array]
    ) -> Union[xr.DataArray, np.array]:
        """
        Approximates the temperature difference between the refrigerant and the
        source.
        Parameters
        ----------
        delta_t_source : Union[xr.DataArray, np.array]
            The temperature difference for the refrigerant source.
        Returns
        -------
        Union[xr.DataArray, np.array]
            The approximate temperature difference between the refrigerant and heat source.
        """
        return delta_t_source / 2
    def _approximate_delta_t_refrigerant_sink(
        self,
        refrigerant: str = "ammonia",
        a: float = {"ammonia": 0.2, "isobutane": -0.0011},
        b: float = {"ammonia": 0.2, "isobutane": 0.3},
        c: float = {"ammonia": 0.016, "isobutane": 2.4},
    ) -> Union[xr.DataArray, np.array]:
        """
        Approximates the temperature difference between the refrigerant and
        heat sink.
        Parameters:
        ----------
        refrigerant : str, optional
            The refrigerant used in the system. Either 'isobutane' or 'ammonia. Default is 'ammonia'.
        a : float, optional
            Coefficient for the temperature difference between the sink and source, default is 0.2.
        b : float, optional
            Coefficient for the temperature difference at the sink, default is 0.2.
        c : float, optional
            Constant term, default is 0.016.
        Returns:
        -------
        Union[xr.DataArray, np.array]
            The approximate temperature difference between the refrigerant and heat sink.
        Notes:
        ------
        This function assumes ammonia as the refrigerant.
        The approximate temperature difference at the refrigerant sink is calculated using the following formula:
        a * (t_sink_out - t_source_out + 2 * delta_t_pinch) + b * delta_t_sink + c
        """
        if refrigerant not in a.keys():
            raise ValueError(
                f"Invalid refrigerant '{refrigerant}'. Must be one of {a.keys()}"
            )
        return (
            a[refrigerant]
            * (self.t_sink_out_kelvin - self.t_source_out + 2 * self.delta_t_pinch)
            + b[refrigerant] * self.delta_t_sink
            + c[refrigerant]
        )
    def _ratio_evaporation_compression_work_approximation(
        self,
        refrigerant: str = "ammonia",
        a: float = {"ammonia": 0.0014, "isobutane": 0.0035},
        b: float = {"ammonia": -0.0015, "isobutane": -0.0033},
        c: float = {"ammonia": 0.039, "isobutane": 0.053},
    ) -> Union[xr.DataArray, np.array]:
        """
        Calculate the ratio of evaporation to compression work approximation.
        Parameters:
        ----------
        refrigerant : str, optional
            The refrigerant used in the system. Either 'isobutane' or 'ammonia. Default is 'ammonia'.
        a : float, optional
            Coefficient 'a' in the approximation equation. Default is 0.0014.
        b : float, optional
            Coefficient 'b' in the approximation equation. Default is -0.0015.
        c : float, optional
            Coefficient 'c' in the approximation equation. Default is 0.039.
        Returns:
        -------
        Union[xr.DataArray, np.array]
            The approximated ratio of evaporation to compression work.
        Notes:
        ------
        This function assumes ammonia as the refrigerant.
        The approximation equation used is:
        ratio = a * (t_sink_out - t_source_out + 2 * delta_t_pinch) + b * delta_t_sink + c
        """
        if refrigerant not in a.keys():
            raise ValueError(
                f"Invalid refrigerant '{refrigerant}'. Must be one of {a.keys()}"
            )
        return (
            a[refrigerant]
            * (self.t_sink_out_kelvin - self.t_source_out + 2 * self.delta_t_pinch)
            + b[refrigerant] * self.delta_t_sink
            + c[refrigerant]
        )
--- a/scripts/build_cop_profiles/DecentralHeatingCopApproximator.py
+++ b/scripts/build_cop_profiles/DecentralHeatingCopApproximator.py
@ -0,0 +1,110 @@
 # -*- coding: utf-8 -*-
 # SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
 #
 # SPDX-License-Identifier: MIT
 from typing import Union
 import numpy as np
 import xarray as xr
 from BaseCopApproximator import BaseCopApproximator
 class DecentralHeatingCopApproximator(BaseCopApproximator):
    """
    Approximate the coefficient of performance (COP) for a heat pump in a
    decentral heating system (individual/household heating).
    Uses a quadratic regression on the temperature difference between the source and sink based on empirical data proposed by Staffell et al. 2012.
    Attributes
    ----------
    forward_temperature_celsius : Union[xr.DataArray, np.array]
        The forward temperature in Celsius.
    source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
        The source inlet temperature in Celsius.
    source_type : str
        The source of the heat pump. Must be either 'air' or 'ground'.
    Methods
    -------
    __init__(forward_temperature_celsius, source_inlet_temperature_celsius, source_type)
        Initialize the DecentralHeatingCopApproximator object.
    approximate_cop()
        Compute the COP values using quadratic regression for air-/ground-source heat pumps.
    _approximate_cop_air_source()
        Evaluate quadratic regression for an air-sourced heat pump.
    _approximate_cop_ground_source()
        Evaluate quadratic regression for a ground-sourced heat pump.
    References
    ----------
    [1] Staffell et al., Energy & Environmental Science 11 (2012): A review of domestic heat pumps, https://doi.org/10.1039/C2EE22653G.
    """
    def __init__(
        self,
        forward_temperature_celsius: Union[xr.DataArray, np.array],
        source_inlet_temperature_celsius: Union[xr.DataArray, np.array],
        source_type: str,
    ):
        """
        Initialize the DecentralHeatingCopApproximator object.
        Parameters
        ----------
        forward_temperature_celsius : Union[xr.DataArray, np.array]
            The forward temperature in Celsius.
        source_inlet_temperature_celsius : Union[xr.DataArray, np.array]
            The source inlet temperature in Celsius.
        source_type : str
            The source of the heat pump. Must be either 'air' or 'ground'.
        """
        self.delta_t = forward_temperature_celsius - source_inlet_temperature_celsius
        if source_type not in ["air", "ground"]:
            raise ValueError("'source_type' must be one of ['air', 'ground']")
        else:
            self.source_type = source_type
    def approximate_cop(self) -> Union[xr.DataArray, np.array]:
        """
        Compute the COP values using quadratic regression for air-/ground-
        source heat pumps.
        Returns
        -------
        Union[xr.DataArray, np.array]
            The calculated COP values.
        """
        if self.source_type == "air":
            return self._approximate_cop_air_source()
        elif self.source_type == "ground":
            return self._approximate_cop_ground_source()
    def _approximate_cop_air_source(self) -> Union[xr.DataArray, np.array]:
        """
        Evaluate quadratic regression for an air-sourced heat pump.
        COP = 6.81 - 0.121 * delta_T + 0.000630 * delta_T^2
        Returns
        -------
        Union[xr.DataArray, np.array]
            The calculated COP values.
        """
        return 6.81 - 0.121 * self.delta_t + 0.000630 * self.delta_t**2
    def _approximate_cop_ground_source(self) -> Union[xr.DataArray, np.array]:
        """
        Evaluate quadratic regression for a ground-sourced heat pump.
        COP = 8.77 - 0.150 * delta_T + 0.000734 * delta_T^2
        Returns
        -------
        Union[xr.DataArray, np.array]
            The calculated COP values.
        """
        return 8.77 - 0.150 * self.delta_t + 0.000734 * self.delta_t**2
--- a/scripts/build_cop_profiles/run.py
+++ b/scripts/build_cop_profiles/run.py
@ -0,0 +1,95 @@
 # -*- coding: utf-8 -*-
 # SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
 #
 # SPDX-License-Identifier: MIT
 import sys
 import numpy as np
 import pandas as pd
 import xarray as xr
 from _helpers import set_scenario_config
 from CentralHeatingCopApproximator import CentralHeatingCopApproximator
 from DecentralHeatingCopApproximator import DecentralHeatingCopApproximator
 from scripts.definitions.heat_system_type import HeatSystemType
 sys.path.append("..")
 def get_cop(
    heat_system_type: str,
    heat_source: str,
    source_inlet_temperature_celsius: xr.DataArray,
 ) -> xr.DataArray:
    """
    Calculate the coefficient of performance (COP) for a heating system.
    Parameters
    ----------
    heat_system_type : str
        The type of heating system.
    heat_source : str
        The heat source used in the heating system.
    source_inlet_temperature_celsius : xr.DataArray
        The inlet temperature of the heat source in Celsius.
    Returns
    -------
    xr.DataArray
        The calculated coefficient of performance (COP) for the heating system.
    """
    if HeatSystemType(heat_system_type).is_central:
        return CentralHeatingCopApproximator(
            forward_temperature_celsius=snakemake.params.forward_temperature_central_heating,
            return_temperature_celsius=snakemake.params.return_temperature_central_heating,
            source_inlet_temperature_celsius=source_inlet_temperature_celsius,
            source_outlet_temperature_celsius=source_inlet_temperature_celsius
            - snakemake.params.heat_source_cooling_central_heating,
        ).approximate_cop()
    else:
        return DecentralHeatingCopApproximator(
            forward_temperature_celsius=snakemake.params.heat_pump_sink_T_decentral_heating,
            source_inlet_temperature_celsius=source_inlet_temperature_celsius,
            source_type=heat_source,
        ).approximate_cop()
 if __name__ == "__main__":
    if "snakemake" not in globals():
        from _helpers import mock_snakemake
        snakemake = mock_snakemake(
            "build_cop_profiles",
            simpl="",
            clusters=48,
        )
    set_scenario_config(snakemake)
    cop_all_system_types = []
    for heat_system_type, heat_sources in snakemake.params.heat_pump_sources.items():
        cop_this_system_type = []
        for heat_source in heat_sources:
            source_inlet_temperature_celsius = xr.open_dataarray(
                snakemake.input[f"temp_{heat_source.replace('ground', 'soil')}_total"]
            )
            cop_da = get_cop(
                heat_system_type=heat_system_type,
                heat_source=heat_source,
                source_inlet_temperature_celsius=source_inlet_temperature_celsius,
            )
            cop_this_system_type.append(cop_da)
        cop_all_system_types.append(
            xr.concat(
                cop_this_system_type, dim=pd.Index(heat_sources, name="heat_source")
            )
        )
    cop_dataarray = xr.concat(
        cop_all_system_types,
        dim=pd.Index(snakemake.params.heat_pump_sources.keys(), name="heat_system"),
    )
    cop_dataarray.to_netcdf(snakemake.output.cop_profiles)
--- a/scripts/build_shapes.py
+++ b/scripts/build_shapes.py
@ -26,7 +26,7 @@ Inputs
    .. image:: img/countries.png
        :scale: 33 %
- ``data/bundle/eez/World_EEZ_v8_2014.shp``: World `exclusive economic zones <https://en.wikipedia.org/wiki/Exclusive_economic_zone>`_ (EEZ)
+- ``data/eez/World_EEZ_v12_20231025_gpkg/eez_v12.gpkg   ``: World `exclusive economic zones <https://en.wikipedia.org/wiki/Exclusive_economic_zone>`_ (EEZ)
    .. image:: img/eez.png
        :scale: 33 %
@ -73,22 +73,16 @@ from functools import reduce
 from itertools import takewhile
 from operator import attrgetter
 import country_converter as coco
 import geopandas as gpd
 import numpy as np
 import pandas as pd
 import pycountry as pyc
 from _helpers import configure_logging, set_scenario_config
 from shapely.geometry import MultiPolygon, Polygon
 logger = logging.getLogger(__name__)
-
+cc = coco.CountryConverter()
 def _get_country(target, **keys):
    assert len(keys) == 1
    try:
        return getattr(pyc.countries.get(**keys), target)
    except (KeyError, AttributeError):
        return np.nan
 def _simplify_polys(polys, minarea=0.1, tolerance=None, filterremote=True):
@ -135,22 +129,15 @@ def countries(naturalearth, country_list):
    return s
-def eez(country_shapes, eez, country_list):
+def eez(eez, country_list):
    df = gpd.read_file(eez)
-    df = df.loc[
+    iso3_list = cc.convert(country_list, src="ISO2", to="ISO3")
-        df["ISO_3digit"].isin(
+    df = df.query("ISO_TER1 in @iso3_list and POL_TYPE == '200NM'").copy()
-            [_get_country("alpha_3", alpha_2=c) for c in country_list]
+    df["name"] = cc.convert(df.ISO_TER1, src="ISO3", to="ISO2")
        )
    ]
    df["name"] = df["ISO_3digit"].map(lambda c: _get_country("alpha_2", alpha_3=c))
    s = df.set_index("name").geometry.map(
        lambda s: _simplify_polys(s, filterremote=False)
    )
-    s = gpd.GeoSeries(
+    s = s.to_frame("geometry").set_crs(df.crs)
        {k: v for k, v in s.items() if v.distance(country_shapes[k]) < 1e-3},
        crs=df.crs,
    )
    s = s.to_frame("geometry")
    s.index.name = "name"
    return s
@ -262,9 +249,7 @@ if __name__ == "__main__":
    country_shapes = countries(snakemake.input.naturalearth, snakemake.params.countries)
    country_shapes.reset_index().to_file(snakemake.output.country_shapes)
-    offshore_shapes = eez(
+    offshore_shapes = eez(snakemake.input.eez, snakemake.params.countries)
        country_shapes, snakemake.input.eez, snakemake.params.countries
    )
    offshore_shapes.reset_index().to_file(snakemake.output.offshore_shapes)
    europe_shape = gpd.GeoDataFrame(
--- a/scripts/definitions/heat_sector.py
+++ b/scripts/definitions/heat_sector.py
@ -0,0 +1,28 @@
 # -*- coding: utf-8 -*-
 # SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
 #
 # SPDX-License-Identifier: MIT
 from enum import Enum
 class HeatSector(Enum):
    """
    Enumeration class representing different heat sectors.
    Attributes:
        RESIDENTIAL (str): Represents the residential heat sector.
        SERVICES (str): Represents the services heat sector.
    """
    RESIDENTIAL = "residential"
    SERVICES = "services"
    def __str__(self) -> str:
        """
        Returns the string representation of the heat sector.
        Returns:
            str: The string representation of the heat sector.
        """
        return self.value
--- a/scripts/definitions/heat_system.py
+++ b/scripts/definitions/heat_system.py
@ -0,0 +1,267 @@
 # -*- coding: utf-8 -*-
 # SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
 #
 # SPDX-License-Identifier: MIT
 from enum import Enum
 from scripts.definitions.heat_sector import HeatSector
 from scripts.definitions.heat_system_type import HeatSystemType
 class HeatSystem(Enum):
    """
    Enumeration representing different heat systems.
    Attributes
    ----------
    RESIDENTIAL_RURAL : str
        Heat system for residential areas in rural locations.
    SERVICES_RURAL : str
        Heat system for service areas in rural locations.
    RESIDENTIAL_URBAN_DECENTRAL : str
        Heat system for residential areas in urban decentralized locations.
    SERVICES_URBAN_DECENTRAL : str
        Heat system for service areas in urban decentralized locations.
    URBAN_CENTRAL : str
        Heat system for urban central areas.
    Methods
    -------
    __str__()
        Returns the string representation of the heat system.
    central_or_decentral()
        Returns whether the heat system is central or decentralized.
    system_type()
        Returns the type of the heat system.
    sector()
        Returns the sector of the heat system.
    rural()
        Returns whether the heat system is for rural areas.
    urban_decentral()
        Returns whether the heat system is for urban decentralized areas.
    urban()
        Returns whether the heat system is for urban areas.
    heat_demand_weighting(urban_fraction=None, dist_fraction=None)
        Calculates the heat demand weighting based on urban fraction and distribution fraction.
    heat_pump_costs_name(heat_source)
        Generates the name for the heat pump costs based on the heat source.
    """
    RESIDENTIAL_RURAL = "residential rural"
    SERVICES_RURAL = "services rural"
    RESIDENTIAL_URBAN_DECENTRAL = "residential urban decentral"
    SERVICES_URBAN_DECENTRAL = "services urban decentral"
    URBAN_CENTRAL = "urban central"
    def __init__(self, *args):
        super().__init__(*args)
    def __str__(self) -> str:
        """
        Returns the string representation of the heat system.
        Returns
        -------
        str
            The string representation of the heat system.
        """
        return self.value
    @property
    def central_or_decentral(self) -> str:
        """
        Returns whether the heat system is central or decentralized.
        Returns
        -------
        str
            "central" if the heat system is central, "decentral" otherwise.
        """
        if self == HeatSystem.URBAN_CENTRAL:
            return "central"
        else:
            return "decentral"
    @property
    def system_type(self) -> HeatSystemType:
        """
        Returns the type of the heat system.
        Returns
        -------
        str
            The type of the heat system.
        Raises
        ------
        RuntimeError
            If the heat system is invalid.
        """
        if self == HeatSystem.URBAN_CENTRAL:
            return HeatSystemType.URBAN_CENTRAL
        elif (
            self == HeatSystem.RESIDENTIAL_URBAN_DECENTRAL
            or self == HeatSystem.SERVICES_URBAN_DECENTRAL
        ):
            return HeatSystemType.URBAN_DECENTRAL
        elif self == HeatSystem.RESIDENTIAL_RURAL or self == HeatSystem.SERVICES_RURAL:
            return HeatSystemType.RURAL
        else:
            raise RuntimeError(f"Invalid heat system: {self}")
    @property
    def sector(self) -> HeatSector:
        """
        Returns the sector of the heat system.
        Returns
        -------
        HeatSector
            The sector of the heat system.
        """
        if (
            self == HeatSystem.RESIDENTIAL_RURAL
            or self == HeatSystem.RESIDENTIAL_URBAN_DECENTRAL
        ):
            return HeatSector.RESIDENTIAL
        elif (
            self == HeatSystem.SERVICES_RURAL
            or self == HeatSystem.SERVICES_URBAN_DECENTRAL
        ):
            return HeatSector.SERVICES
        else:
            "tot"
    @property
    def is_rural(self) -> bool:
        """
        Returns whether the heat system is for rural areas.
        Returns
        -------
        bool
            True if the heat system is for rural areas, False otherwise.
        """
        if self == HeatSystem.RESIDENTIAL_RURAL or self == HeatSystem.SERVICES_RURAL:
            return True
        else:
            return False
    @property
    def is_urban_decentral(self) -> bool:
        """
        Returns whether the heat system is for urban decentralized areas.
        Returns
        -------
        bool
            True if the heat system is for urban decentralized areas, False otherwise.
        """
        if (
            self == HeatSystem.RESIDENTIAL_URBAN_DECENTRAL
            or self == HeatSystem.SERVICES_URBAN_DECENTRAL
        ):
            return True
        else:
            return False
    @property
    def is_urban(self) -> bool:
        """
        Returns whether the heat system is for urban areas.
        Returns
        -------
        bool True if the heat system is for urban areas, False otherwise.
        """
        return not self.is_rural
    def heat_demand_weighting(self, urban_fraction=None, dist_fraction=None) -> float:
        """
        Calculates the heat demand weighting based on urban fraction and
        distribution fraction.
        Parameters
        ----------
        urban_fraction : float, optional
            The fraction of urban heat demand.
        dist_fraction : float, optional
            The fraction of distributed heat demand.
        Returns
        -------
        float
            The heat demand weighting.
        Raises
        ------
        RuntimeError
            If the heat system is invalid.
        """
        if "rural" in self.value:
            return 1 - urban_fraction
        elif "urban central" in self.value:
            return dist_fraction
        elif "urban decentral" in self.value:
            return urban_fraction - dist_fraction
        else:
            raise RuntimeError(f"Invalid heat system: {self}")
    def heat_pump_costs_name(self, heat_source: str) -> str:
        """
        Generates the name for the heat pump costs based on the heat source and
        system.
        Used to retrieve data from `technology-data <https://github.com/PyPSA/technology-data>`.
        Parameters
        ----------
        heat_source : str
            The heat source.
        Returns
        -------
        str
            The name for the heat pump costs.
        """
        return f"{self.central_or_decentral} {heat_source}-sourced heat pump"
    @property
    def resistive_heater_costs_name(self) -> str:
        """
        Generates the name for the resistive heater costs based on the heat
        system.
        Used to retrieve data from `technology-data <https://github.com/PyPSA/technology-data>`.
        Returns
        -------
        str
            The name for the heater costs.
        """
        return f"{self.central_or_decentral} resistive heater"
    @property
    def gas_boiler_costs_name(self) -> str:
        """
        Generates the name for the gas boiler costs based on the heat system.
        Used to retrieve data from `technology-data <https://github.com/PyPSA/technology-data>`.
        Returns
        -------
        str
            The name for the gas boiler costs.
        """
        return f"{self.central_or_decentral} gas boiler"
    @property
    def oil_boiler_costs_name(self) -> str:
        """
        Generates the name for the oil boiler costs based on the heat system.
        Used to retrieve data from `technology-data <https://github.com/PyPSA/technology-data>`.
        Returns
        -------
        str
            The name for the oil boiler costs.
        """
        return "decentral oil boiler"
--- a/scripts/definitions/heat_system_type.py
+++ b/scripts/definitions/heat_system_type.py
@ -0,0 +1,35 @@
 # -*- coding: utf-8 -*-
 # SPDX-FileCopyrightText: : 2020-2024 The PyPSA-Eur Authors
 #
 # SPDX-License-Identifier: MIT
 from enum import Enum
 class HeatSystemType(Enum):
    """
    Enumeration representing different types of heat systems.
    """
    URBAN_CENTRAL = "urban central"
    URBAN_DECENTRAL = "urban decentral"
    RURAL = "rural"
    def __str__(self) -> str:
        """
        Returns the string representation of the heat system type.
        Returns:
            str: The string representation of the heat system type.
        """
        return self.value
    @property
    def is_central(self) -> bool:
        """
        Returns whether the heat system type is central.
        Returns:
            bool: True if the heat system type is central, False otherwise.
        """
        return self == HeatSystemType.URBAN_CENTRAL
--- a/scripts/prepare_links_p_nom.py
+++ b/scripts/prepare_links_p_nom.py
@ -1,95 +0,0 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 # SPDX-FileCopyrightText: : 2017-2024 The PyPSA-Eur Authors
 #
 # SPDX-License-Identifier: MIT
 """
 Extracts capacities of HVDC links from `Wikipedia.
 <https://en.wikipedia.org/wiki/List_of_HVDC_projects>`_.
 Relevant Settings
 -----------------
 .. code:: yaml
    enable:
        prepare_links_p_nom:
 .. seealso::
    Documentation of the configuration file ``config/config.yaml`` at
    :ref:`toplevel_cf`
 Inputs
 ------
 *None*
 Outputs
 -------
 - ``data/links_p_nom.csv``: A plain download of https://en.wikipedia.org/wiki/List_of_HVDC_projects#Europe plus extracted coordinates.
 Description
 -----------
 *None*
 """
 import logging
 import pandas as pd
 from _helpers import configure_logging, set_scenario_config
 logger = logging.getLogger(__name__)
 def multiply(s):
    return s.str[0].astype(float) * s.str[1].astype(float)
 def extract_coordinates(s):
    regex = (
        r"(\d{1,2})°(\d{1,2})′(\d{1,2})″(N|S) " r"(\d{1,2})°(\d{1,2})′(\d{1,2})″(E|W)"
    )
    e = s.str.extract(regex, expand=True)
    lat = (
        e[0].astype(float) + (e[1].astype(float) + e[2].astype(float) / 60.0) / 60.0
    ) * e[3].map({"N": +1.0, "S": -1.0})
    lon = (
        e[4].astype(float) + (e[5].astype(float) + e[6].astype(float) / 60.0) / 60.0
    ) * e[7].map({"E": +1.0, "W": -1.0})
    return lon, lat
 if __name__ == "__main__":
    if "snakemake" not in globals():
        from _helpers import mock_snakemake  # rule must be enabled in config
        snakemake = mock_snakemake("prepare_links_p_nom", simpl="")
    configure_logging(snakemake)
    set_scenario_config(snakemake)
    links_p_nom = pd.read_html(
        "https://en.wikipedia.org/wiki/List_of_HVDC_projects", header=0, match="SwePol"
    )[0]
    mw = "Power (MW)"
    m_b = links_p_nom[mw].str.contains("x").fillna(False)
    links_p_nom.loc[m_b, mw] = links_p_nom.loc[m_b, mw].str.split("x").pipe(multiply)
    links_p_nom[mw] = (
        links_p_nom[mw].str.extract("[-/]?([\d.]+)", expand=False).astype(float)
    )
    links_p_nom["x1"], links_p_nom["y1"] = extract_coordinates(
        links_p_nom["Converterstation 1"]
    )
    links_p_nom["x2"], links_p_nom["y2"] = extract_coordinates(
        links_p_nom["Converterstation 2"]
    )
    links_p_nom.dropna(subset=["x1", "y1", "x2", "y2"]).to_csv(
        snakemake.output[0], index=False
    )
--- a/scripts/prepare_sector_network.py
+++ b/scripts/prepare_sector_network.py
@ -37,6 +37,10 @@ from pypsa.geo import haversine_pts
 from pypsa.io import import_components_from_dataframe
 from scipy.stats import beta
 from scripts.definitions.heat_sector import HeatSector
 from scripts.definitions.heat_system import HeatSystem
 from scripts.definitions.heat_system_type import HeatSystemType
 spatial = SimpleNamespace()
 logger = logging.getLogger(__name__)
@ -56,19 +60,27 @@ def define_spatial(nodes, options):
    # biomass
    spatial.biomass = SimpleNamespace()
    spatial.msw = SimpleNamespace()
    if options.get("biomass_spatial", options["biomass_transport"]):
        spatial.biomass.nodes = nodes + " solid biomass"
        spatial.biomass.bioliquids = nodes + " bioliquids"
        spatial.biomass.locations = nodes
        spatial.biomass.industry = nodes + " solid biomass for industry"
        spatial.biomass.industry_cc = nodes + " solid biomass for industry CC"
        spatial.msw.nodes = nodes + " municipal solid waste"
        spatial.msw.locations = nodes
    else:
        spatial.biomass.nodes = ["EU solid biomass"]
        spatial.biomass.bioliquids = ["EU unsustainable bioliquids"]
        spatial.biomass.locations = ["EU"]
        spatial.biomass.industry = ["solid biomass for industry"]
        spatial.biomass.industry_cc = ["solid biomass for industry CC"]
        spatial.msw.nodes = ["EU municipal solid waste"]
        spatial.msw.locations = ["EU"]
    spatial.biomass.df = pd.DataFrame(vars(spatial.biomass), index=nodes)
    spatial.msw.df = pd.DataFrame(vars(spatial.msw), index=nodes)
    # co2
@ -2102,7 +2114,7 @@ def build_heat_demand(n):
        .unstack(level=1)
    )
-    sectors = ["residential", "services"]
+    sectors = [sector.value for sector in HeatSector]
    uses = ["water", "space"]
    heat_demand = {}
@ -2130,10 +2142,21 @@ def build_heat_demand(n):
    return heat_demand
-def add_heat(n, costs):
+def add_heat(n: pypsa.Network, costs: pd.DataFrame, cop: xr.DataArray):
    """
    Add heat sector to the network.
    Parameters:
        n (pypsa.Network): The PyPSA network object.
        costs (pd.DataFrame): DataFrame containing cost information.
        cop (xr.DataArray): DataArray containing coefficient of performance (COP) values.
    Returns:
        None
    """
    logger.info("Add heat sector")
-    sectors = ["residential", "services"]
+    sectors = [sector.value for sector in HeatSector]
    heat_demand = build_heat_demand(n)
@ -2152,23 +2175,6 @@ def add_heat(n, costs):
        for sector in sectors:
            heat_demand[sector + " space"] = (1 - dE) * heat_demand[sector + " space"]
    heat_systems = [
        "residential rural",
        "services rural",
        "residential urban decentral",
        "services urban decentral",
        "urban central",
    ]
    cop = {
        "air": xr.open_dataarray(snakemake.input.cop_air_total)
        .to_pandas()
        .reindex(index=n.snapshots),
        "ground": xr.open_dataarray(snakemake.input.cop_soil_total)
        .to_pandas()
        .reindex(index=n.snapshots),
    }
    if options["solar_thermal"]:
        solar_thermal = (
            xr.open_dataarray(snakemake.input.solar_thermal_total)
@ -2178,31 +2184,34 @@ def add_heat(n, costs):
        # 1e3 converts from W/m^2 to MW/(1000m^2) = kW/m^2
        solar_thermal = options["solar_cf_correction"] * solar_thermal / 1e3
-    for name in heat_systems:
+    for (
-        name_type = "central" if name == "urban central" else "decentral"
+        heat_system
    ) in (
        HeatSystem
    ):  # this loops through all heat systems defined in _entities.HeatSystem
-        if name == "urban central":
+        if heat_system == HeatSystem.URBAN_CENTRAL:
            nodes = dist_fraction.index[dist_fraction > 0]
        else:
            nodes = pop_layout.index
-        n.add("Carrier", name + " heat")
+        n.add("Carrier", f"{heat_system} heat")
        n.madd(
            "Bus",
-            nodes + f" {name} heat",
+            nodes + f" {heat_system.value} heat",
            location=nodes,
-            carrier=name + " heat",
+            carrier=f"{heat_system.value} heat",
            unit="MWh_th",
        )
-        if name == "urban central" and options.get("central_heat_vent"):
+        if heat_system == HeatSystem.URBAN_CENTRAL and options.get("central_heat_vent"):
            n.madd(
                "Generator",
-                nodes + f" {name} heat vent",
+                nodes + f" {heat_system} heat vent",
-                bus=nodes + f" {name} heat",
+                bus=nodes + f" {heat_system} heat",
                location=nodes,
-                carrier=name + " heat vent",
+                carrier=f"{heat_system} heat vent",
                p_nom_extendable=True,
                p_max_pu=0,
                p_min_pu=-1,
@ -2210,30 +2219,24 @@ def add_heat(n, costs):
            )
        ## Add heat load
        factor = heat_system.heat_demand_weighting(
            urban_fraction=urban_fraction[nodes], dist_fraction=dist_fraction[nodes]
        )
        if not heat_system == HeatSystem.URBAN_CENTRAL:
            heat_load = (
                heat_demand[
                    [
                        heat_system.sector.value + " water",
                        heat_system.sector.value + " space",
                    ]
                ]
                .T.groupby(level=1)
                .sum()
                .T[nodes]
                .multiply(factor)
            )
-        for sector in sectors:
+        if heat_system == HeatSystem.URBAN_CENTRAL:
            # heat demand weighting
            if "rural" in name:
                factor = 1 - urban_fraction[nodes]
            elif "urban central" in name:
                factor = dist_fraction[nodes]
            elif "urban decentral" in name:
                factor = urban_fraction[nodes] - dist_fraction[nodes]
            else:
                raise NotImplementedError(
                    f" {name} not in " f"heat systems: {heat_systems}"
                )
            if sector in name:
                heat_load = (
                    heat_demand[[sector + " water", sector + " space"]]
                    .T.groupby(level=1)
                    .sum()
                    .T[nodes]
                    .multiply(factor)
                )
        if name == "urban central":
            heat_load = (
                heat_demand.T.groupby(level=1)
                .sum()
@ -2246,20 +2249,25 @@ def add_heat(n, costs):
        n.madd(
            "Load",
            nodes,
-            suffix=f" {name} heat",
+            suffix=f" {heat_system} heat",
-            bus=nodes + f" {name} heat",
+            bus=nodes + f" {heat_system} heat",
-            carrier=name + " heat",
+            carrier=f"{heat_system} heat",
            p_set=heat_load,
        )
        ## Add heat pumps
-
+        for heat_source in snakemake.params.heat_pump_sources[
-        heat_pump_types = ["air"] if "urban" in name else ["ground", "air"]
+            heat_system.system_type.value
-
+        ]:
-        for heat_pump_type in heat_pump_types:
+            costs_name = heat_system.heat_pump_costs_name(heat_source)
            costs_name = f"{name_type} {heat_pump_type}-sourced heat pump"
            efficiency = (
-                cop[heat_pump_type][nodes]
+                cop.sel(
                    heat_system=heat_system.system_type.value,
                    heat_source=heat_source,
                    name=nodes,
                )
                .to_pandas()
                .reindex(index=n.snapshots)
                if options["time_dep_hp_cop"]
                else costs.at[costs_name, "efficiency"]
            )
@ -2267,10 +2275,10 @@ def add_heat(n, costs):
            n.madd(
                "Link",
                nodes,
-                suffix=f" {name} {heat_pump_type} heat pump",
+                suffix=f" {heat_system} {heat_source} heat pump",
                bus0=nodes,
-                bus1=nodes + f" {name} heat",
+                bus1=nodes + f" {heat_system} heat",
-                carrier=f"{name} {heat_pump_type} heat pump",
+                carrier=f"{heat_system} {heat_source} heat pump",
                efficiency=efficiency,
                capital_cost=costs.at[costs_name, "efficiency"]
                * costs.at[costs_name, "fixed"]
@ -2280,59 +2288,65 @@ def add_heat(n, costs):
            )
        if options["tes"]:
-            n.add("Carrier", name + " water tanks")
+            n.add("Carrier", f"{heat_system} water tanks")
            n.madd(
                "Bus",
-                nodes + f" {name} water tanks",
+                nodes + f" {heat_system} water tanks",
                location=nodes,
-                carrier=name + " water tanks",
+                carrier=f"{heat_system} water tanks",
                unit="MWh_th",
            )
            n.madd(
                "Link",
-                nodes + f" {name} water tanks charger",
+                nodes + f" {heat_system} water tanks charger",
-                bus0=nodes + f" {name} heat",
+                bus0=nodes + f" {heat_system} heat",
-                bus1=nodes + f" {name} water tanks",
+                bus1=nodes + f" {heat_system} water tanks",
                efficiency=costs.at["water tank charger", "efficiency"],
-                carrier=name + " water tanks charger",
+                carrier=f"{heat_system} water tanks charger",
                p_nom_extendable=True,
            )
            n.madd(
                "Link",
-                nodes + f" {name} water tanks discharger",
+                nodes + f" {heat_system} water tanks discharger",
-                bus0=nodes + f" {name} water tanks",
+                bus0=nodes + f" {heat_system} water tanks",
-                bus1=nodes + f" {name} heat",
+                bus1=nodes + f" {heat_system} heat",
-                carrier=name + " water tanks discharger",
+                carrier=f"{heat_system} water tanks discharger",
                efficiency=costs.at["water tank discharger", "efficiency"],
                p_nom_extendable=True,
            )
-            tes_time_constant_days = options["tes_tau"][name_type]
+            tes_time_constant_days = options["tes_tau"][
                heat_system.central_or_decentral
            ]
            n.madd(
                "Store",
-                nodes + f" {name} water tanks",
+                nodes + f" {heat_system} water tanks",
-                bus=nodes + f" {name} water tanks",
+                bus=nodes + f" {heat_system} water tanks",
                e_cyclic=True,
                e_nom_extendable=True,
-                carrier=name + " water tanks",
+                carrier=f"{heat_system} water tanks",
                standing_loss=1 - np.exp(-1 / 24 / tes_time_constant_days),
-                capital_cost=costs.at[name_type + " water tank storage", "fixed"],
+                capital_cost=costs.at[
-                lifetime=costs.at[name_type + " water tank storage", "lifetime"],
+                    heat_system.central_or_decentral + " water tank storage", "fixed"
                ],
                lifetime=costs.at[
                    heat_system.central_or_decentral + " water tank storage", "lifetime"
                ],
            )
        if options["resistive_heaters"]:
-            key = f"{name_type} resistive heater"
+            key = f"{heat_system.central_or_decentral} resistive heater"
            n.madd(
                "Link",
-                nodes + f" {name} resistive heater",
+                nodes + f" {heat_system} resistive heater",
                bus0=nodes,
-                bus1=nodes + f" {name} heat",
+                bus1=nodes + f" {heat_system} heat",
-                carrier=name + " resistive heater",
+                carrier=f"{heat_system} resistive heater",
                efficiency=costs.at[key, "efficiency"],
                capital_cost=costs.at[key, "efficiency"]
                * costs.at[key, "fixed"]
@ -2342,16 +2356,16 @@ def add_heat(n, costs):
            )
        if options["boilers"]:
-            key = f"{name_type} gas boiler"
+            key = f"{heat_system.central_or_decentral} gas boiler"
            n.madd(
                "Link",
-                nodes + f" {name} gas boiler",
+                nodes + f" {heat_system} gas boiler",
                p_nom_extendable=True,
                bus0=spatial.gas.df.loc[nodes, "nodes"].values,
-                bus1=nodes + f" {name} heat",
+                bus1=nodes + f" {heat_system} heat",
                bus2="co2 atmosphere",
-                carrier=name + " gas boiler",
+                carrier=f"{heat_system} gas boiler",
                efficiency=costs.at[key, "efficiency"],
                efficiency2=costs.at["gas", "CO2 intensity"],
                capital_cost=costs.at[key, "efficiency"]
@ -2361,22 +2375,26 @@ def add_heat(n, costs):
            )
        if options["solar_thermal"]:
-            n.add("Carrier", name + " solar thermal")
+            n.add("Carrier", f"{heat_system} solar thermal")
            n.madd(
                "Generator",
                nodes,
-                suffix=f" {name} solar thermal collector",
+                suffix=f" {heat_system} solar thermal collector",
-                bus=nodes + f" {name} heat",
+                bus=nodes + f" {heat_system} heat",
-                carrier=name + " solar thermal",
+                carrier=f"{heat_system} solar thermal",
                p_nom_extendable=True,
-                capital_cost=costs.at[name_type + " solar thermal", "fixed"]
+                capital_cost=costs.at[
                    heat_system.central_or_decentral + " solar thermal", "fixed"
                ]
                * overdim_factor,
                p_max_pu=solar_thermal[nodes],
-                lifetime=costs.at[name_type + " solar thermal", "lifetime"],
+                lifetime=costs.at[
                    heat_system.central_or_decentral + " solar thermal", "lifetime"
                ],
            )
-        if options["chp"] and name == "urban central":
+        if options["chp"] and heat_system == HeatSystem.URBAN_CENTRAL:
            # add gas CHP; biomass CHP is added in biomass section
            n.madd(
                "Link",
@ -2433,16 +2451,20 @@ def add_heat(n, costs):
                lifetime=costs.at["central gas CHP", "lifetime"],
            )
-        if options["chp"] and options["micro_chp"] and name != "urban central":
+        if (
            options["chp"]
            and options["micro_chp"]
            and heat_system.value != "urban central"
        ):
            n.madd(
                "Link",
-                nodes + f" {name} micro gas CHP",
+                nodes + f" {heat_system} micro gas CHP",
                p_nom_extendable=True,
                bus0=spatial.gas.df.loc[nodes, "nodes"].values,
                bus1=nodes,
-                bus2=nodes + f" {name} heat",
+                bus2=nodes + f" {heat_system} heat",
                bus3="co2 atmosphere",
-                carrier=name + " micro gas CHP",
+                carrier=heat_system.value + " micro gas CHP",
                efficiency=costs.at["micro CHP", "efficiency"],
                efficiency2=costs.at["micro CHP", "efficiency-heat"],
                efficiency3=costs.at["gas", "CO2 intensity"],
@ -2478,7 +2500,7 @@ def add_heat(n, costs):
        ) / heat_demand.T.groupby(level=[1]).sum().T
        for name in n.loads[
-            n.loads.carrier.isin([x + " heat" for x in heat_systems])
+            n.loads.carrier.isin([x + " heat" for x in HeatSystem])
        ].index:
            node = n.buses.loc[name, "location"]
            ct = pop_layout.loc[node, "ct"]
@ -2632,19 +2654,83 @@ def add_biomass(n, costs):
        biogas_potentials_spatial = biomass_potentials["biogas"].rename(
            index=lambda x: x + " biogas"
        )
        unsustainable_biogas_potentials_spatial = biomass_potentials[
            "unsustainable biogas"
        ].rename(index=lambda x: x + " biogas")
    else:
        biogas_potentials_spatial = biomass_potentials["biogas"].sum()
        unsustainable_biogas_potentials_spatial = biomass_potentials[
            "unsustainable biogas"
        ].sum()
    if options.get("biomass_spatial", options["biomass_transport"]):
        solid_biomass_potentials_spatial = biomass_potentials["solid biomass"].rename(
            index=lambda x: x + " solid biomass"
        )
        msw_biomass_potentials_spatial = biomass_potentials[
            "municipal solid waste"
        ].rename(index=lambda x: x + " municipal solid waste")
        unsustainable_solid_biomass_potentials_spatial = biomass_potentials[
            "unsustainable solid biomass"
        ].rename(index=lambda x: x + " solid biomass")
    else:
        solid_biomass_potentials_spatial = biomass_potentials["solid biomass"].sum()
        msw_biomass_potentials_spatial = biomass_potentials[
            "municipal solid waste"
        ].sum()
        unsustainable_solid_biomass_potentials_spatial = biomass_potentials[
            "unsustainable solid biomass"
        ].sum()
    if options["regional_oil_demand"]:
        unsustainable_liquid_biofuel_potentials_spatial = biomass_potentials[
            "unsustainable bioliquids"
        ].rename(index=lambda x: x + " bioliquids")
    else:
        unsustainable_liquid_biofuel_potentials_spatial = biomass_potentials[
            "unsustainable bioliquids"
        ].sum()
    n.add("Carrier", "biogas")
    n.add("Carrier", "solid biomass")
    if (
        options["municipal_solid_waste"]
        and not options["industry"]
        and cf_industry["waste_to_energy"]
        or cf_industry["waste_to_energy_cc"]
    ):
        logger.warning(
            "Flag municipal_solid_waste can be only used with industry "
            "sector waste to energy."
            "Setting municipal_solid_waste=False."
        )
        options["municipal_solid_waste"] = False
    if options["municipal_solid_waste"]:
        n.add("Carrier", "municipal solid waste")
        n.madd(
            "Bus",
            spatial.msw.nodes,
            location=spatial.msw.locations,
            carrier="municipal solid waste",
        )
        e_max_pu = pd.Series([1] * (len(n.snapshots) - 1) + [0], index=n.snapshots)
        n.madd(
            "Store",
            spatial.msw.nodes,
            bus=spatial.msw.nodes,
            carrier="municipal solid waste",
            e_nom=msw_biomass_potentials_spatial,
            marginal_cost=0,  # costs.at["municipal solid waste", "fuel"],
            e_max_pu=e_max_pu,
            e_initial=msw_biomass_potentials_spatial,
        )
    n.madd(
        "Bus",
        spatial.gas.biogas,
@ -2729,6 +2815,81 @@ def add_biomass(n, costs):
            p_nom_extendable=True,
        )
    if biomass_potentials.filter(like="unsustainable").sum().sum() > 0:
        # Create timeseries to force usage of unsustainable potentials
        e_max_pu = pd.DataFrame(1, index=n.snapshots, columns=spatial.gas.biogas)
        e_max_pu.iloc[-1] = 0
        n.madd(
            "Store",
            spatial.gas.biogas,
            suffix=" unsustainable",
            bus=spatial.gas.biogas,
            carrier="unsustainable biogas",
            e_nom=unsustainable_biogas_potentials_spatial,
            marginal_cost=costs.at["biogas", "fuel"],
            e_initial=unsustainable_biogas_potentials_spatial,
            e_nom_extendable=False,
            e_max_pu=e_max_pu,
        )
        e_max_pu = pd.DataFrame(1, index=n.snapshots, columns=spatial.biomass.nodes)
        e_max_pu.iloc[-1] = 0
        n.madd(
            "Store",
            spatial.biomass.nodes,
            suffix=" unsustainable",
            bus=spatial.biomass.nodes,
            carrier="unsustainable solid biomass",
            e_nom=unsustainable_solid_biomass_potentials_spatial,
            marginal_cost=costs.at["fuelwood", "fuel"],
            e_initial=unsustainable_solid_biomass_potentials_spatial,
            e_nom_extendable=False,
            e_max_pu=e_max_pu,
        )
        n.madd(
            "Bus",
            spatial.biomass.bioliquids,
            location=spatial.biomass.locations,
            carrier="unsustainable bioliquids",
            unit="MWh_LHV",
        )
        e_max_pu = pd.DataFrame(
            1, index=n.snapshots, columns=spatial.biomass.bioliquids
        )
        e_max_pu.iloc[-1] = 0
        n.madd(
            "Store",
            spatial.biomass.bioliquids,
            suffix=" unsustainable",
            bus=spatial.biomass.bioliquids,
            carrier="unsustainable bioliquids",
            e_nom=unsustainable_liquid_biofuel_potentials_spatial,
            marginal_cost=costs.at["biodiesel crops", "fuel"],
            e_initial=unsustainable_liquid_biofuel_potentials_spatial,
            e_nom_extendable=False,
            e_max_pu=e_max_pu,
        )
        n.madd(
            "Link",
            spatial.biomass.bioliquids,
            bus0=spatial.biomass.bioliquids,
            bus1=spatial.oil.nodes,
            bus2="co2 atmosphere",
            carrier="unsustainable bioliquids",
            efficiency=1,
            efficiency2=-costs.at["solid biomass", "CO2 intensity"]
            + costs.at["BtL", "CO2 stored"],
            p_nom=unsustainable_liquid_biofuel_potentials_spatial,
            marginal_cost=costs.at["BtL", "VOM"],
        )
    n.madd(
        "Link",
        spatial.gas.biogas_to_gas,
@ -2800,6 +2961,19 @@ def add_biomass(n, costs):
            carrier="solid biomass transport",
        )
        if options["municipal_solid_waste"]:
            n.madd(
                "Link",
                biomass_transport.index,
                bus0=biomass_transport.bus0 + " municipal solid waste",
                bus1=biomass_transport.bus1 + " municipal solid waste",
                p_nom_extendable=False,
                p_nom=5e4,
                length=biomass_transport.length.values,
                marginal_cost=biomass_transport.costs * biomass_transport.length.values,
                carrier="municipal solid waste transport",
            )
    elif options["biomass_spatial"]:
        # add artificial biomass generators at nodes which include transport costs
        transport_costs = pd.read_csv(
@ -2829,6 +3003,26 @@ def add_biomass(n, costs):
            type="operational_limit",
        )
        if options["municipal_solid_waste"]:
            # Add municipal solid waste
            n.madd(
                "Generator",
                spatial.msw.nodes,
                bus=spatial.msw.nodes,
                carrier="municipal solid waste",
                p_nom=10000,
                marginal_cost=0  # costs.at["municipal solid waste", "fuel"]
                + bus_transport_costs * average_distance,
            )
            n.add(
                "GlobalConstraint",
                "msw limit",
                carrier_attribute="municipal solid waste",
                sense="<=",
                constant=biomass_potentials["municipal solid waste"].sum(),
                type="operational_limit",
            )
    # AC buses with district heating
    urban_central = n.buses.index[n.buses.carrier == "urban central heat"]
    if not urban_central.empty and options["chp"]:
@ -3420,27 +3614,23 @@ def add_industry(n, costs):
    if options["oil_boilers"]:
        nodes = pop_layout.index
-        for name in [
+        for heat_system in HeatSystem:
-            "residential rural",
+            if not heat_system == HeatSystem.URBAN_CENTRAL:
-            "services rural",
+                n.madd(
-            "residential urban decentral",
+                    "Link",
-            "services urban decentral",
+                    nodes + f" {heat_system} oil boiler",
-        ]:
+                    p_nom_extendable=True,
-            n.madd(
+                    bus0=spatial.oil.nodes,
-                "Link",
+                    bus1=nodes + f" {heat_system} heat",
-                nodes + f" {name} oil boiler",
+                    bus2="co2 atmosphere",
-                p_nom_extendable=True,
+                    carrier=f"{heat_system} oil boiler",
-                bus0=spatial.oil.nodes,
+                    efficiency=costs.at["decentral oil boiler", "efficiency"],
-                bus1=nodes + f" {name} heat",
+                    efficiency2=costs.at["oil", "CO2 intensity"],
-                bus2="co2 atmosphere",
+                    capital_cost=costs.at["decentral oil boiler", "efficiency"]
-                carrier=f"{name} oil boiler",
+                    * costs.at["decentral oil boiler", "fixed"]
-                efficiency=costs.at["decentral oil boiler", "efficiency"],
+                    * options["overdimension_individual_heating"],
-                efficiency2=costs.at["oil", "CO2 intensity"],
+                    lifetime=costs.at["decentral oil boiler", "lifetime"],
-                capital_cost=costs.at["decentral oil boiler", "efficiency"]
+                )
                * costs.at["decentral oil boiler", "fixed"]
                * options["overdimension_individual_heating"],
                lifetime=costs.at["decentral oil boiler", "lifetime"],
            )
    n.madd(
        "Link",
@ -3535,6 +3725,17 @@ def add_industry(n, costs):
            efficiency3=process_co2_per_naphtha,
        )
        if options.get("biomass", True) and options["municipal_solid_waste"]:
            n.madd(
                "Link",
                spatial.msw.locations,
                bus0=spatial.msw.nodes,
                bus1=non_sequestered_hvc_locations,
                carrier="municipal solid waste",
                p_nom_extendable=True,
                efficiency=1.0,
            )
        n.madd(
            "Link",
            spatial.oil.demand_locations,
@ -4427,13 +4628,14 @@ def add_enhanced_geothermal(n, egs_potentials, egs_overlap, costs):
 # %%
 if __name__ == "__main__":
    if "snakemake" not in globals():
        from _helpers import mock_snakemake
        snakemake = mock_snakemake(
            "prepare_sector_network",
            simpl="",
            opts="",
-            clusters="1",
+            clusters="37",
            ll="vopt",
            sector_opts="",
            planning_horizons="2050",
@ -4492,7 +4694,7 @@ if __name__ == "__main__":
        add_land_transport(n, costs)
    if options["heating"]:
-        add_heat(n, costs)
+        add_heat(n=n, costs=costs, cop=xr.open_dataarray(snakemake.input.cop_profiles))
    if options["biomass"]:
        add_biomass(n, costs)