fix sector coupeled model for single, non EU28 countries - e.g. UA (#1)

* fix sector coupeled model for single, non EU28 countries - e.g. UA

* Modification to industrial energy demand to use 0s if sector data unavailable

---------

Co-authored-by: Sermisha <sermisha.narayana@openenergytransition.org>

doc/release_notes.rst

@@ -19,6 +19,9 @@ Upcoming Release
   <https://www.sciencedirect.com/science/article/pii/S0306261920312551>`__.
   See configuration ``sector: enhanced_geothermal`` for details; by default switched off.
 
+* Allow running the sector model for isolated non-EU28 countries, by filling missing sectoral
+  data with defaults, average EU values or zeros, if not available.
+
 
 PyPSA-Eur 0.11.0 (25th May 2024)
 =====================================

scripts/build_district_heat_share.py

@@ -33,9 +33,11 @@
     pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0)
 
     year = str(snakemake.params.energy_totals_year)
-    district_heat_share = pd.read_csv(snakemake.input.district_heat_share, index_col=0)[
-        year
-    ]
+    district_heat_share = pd.read_csv(snakemake.input.district_heat_share, index_col=0)
+    if not district_heat_share.empty:
+        district_heat_share = district_heat_share[year]
+    else:
+        district_heat_share = pd.Series(index=pop_layout.index, data=0)
 
     # make ct-based share nodal
     district_heat_share = district_heat_share.reindex(pop_layout.ct).fillna(0)
@@ -48,15 +50,16 @@
     pop_layout["urban_ct_fraction"] = pop_layout.urban / pop_layout.ct.map(ct_urban.get)
 
     # fraction of node that is urban
-    urban_fraction = pop_layout.urban / pop_layout[["rural", "urban"]].sum(axis=1)
+    urban_fraction = (pop_layout.urban / pop_layout[["rural", "urban"]].sum(axis=1)).fillna(0)
+
 
     # maximum potential of urban demand covered by district heating
     central_fraction = snakemake.config["sector"]["district_heating"]["potential"]
 
     # district heating share at each node
     dist_fraction_node = (
         district_heat_share * pop_layout["urban_ct_fraction"] / pop_layout["fraction"]
-    )
+    ).fillna(0)
 
     # if district heating share larger than urban fraction -> set urban
     # fraction to district heating share

scripts/build_energy_totals.py

@@ -1026,7 +1026,40 @@ def update_residential_from_eurostat(energy):
         disable_progressbar=snakemake.config["run"].get("disable_progressbar", False),
     )
     swiss = build_swiss()
-    idees = build_idees(idees_countries)
+    if len(idees_countries) > 0:
+        idees = build_idees(idees_countries)
+    else:
+        # e.g. UA and MD
+        logger.info(f"No IDEES data available for {countries} and years 2000-2015. Filling with zeros.")
+        years = range(2000, 2016)
+        idees = pd.DataFrame(
+            index=pd.MultiIndex.from_tuples([(country, year) for country in countries for year in years]),
+            columns=[
+                "passenger cars",
+                "passenger car efficiency",
+                "total passenger cars",
+                "total other road passenger",
+                "total light duty road freight",
+                "total two-wheel",
+                "total heavy duty road freight",
+                "electricity passenger cars",
+                "electricity other road passenger",
+                "electricity light duty road freight",
+                "total rail passenger",
+                "total rail freight",
+                "electricity rail passenger",
+                "electricity rail freight",
+                "total domestic aviation passenger",
+                "total domestic aviation freight",
+                "total international aviation passenger",
+                "total international aviation freight",
+                "derived heat residential",
+                "derived heat services",
+                "thermal uses residential",
+                "thermal uses services",
+            ],
+            data=0
+        )
 
     energy = build_energy_totals(countries, eurostat, swiss, idees)
 

scripts/build_industrial_energy_demand_per_country_today.py

@@ -65,6 +65,7 @@
 
 import country_converter as coco
 import pandas as pd
+import numpy as np
 from _helpers import set_scenario_config
 from tqdm import tqdm
 
@@ -195,12 +196,15 @@ def add_non_eu28_industrial_energy_demand(countries, demand, production):
         return demand
 
     eu28_production = production.loc[countries.intersection(eu28)].sum()
+    if eu28_production.sum() == 0:
+        logger.info("EU production is zero. Fallback: Filling non EU28 countries with zeros.")
     eu28_energy = demand.groupby(level=1).sum()
     eu28_averages = eu28_energy / eu28_production
 
     demand_non_eu28 = pd.concat(
         {k: v * eu28_averages for k, v in production.loc[non_eu28].iterrows()}
-    )
+    ).fillna(0)
+    demand_non_eu28.replace([np.inf, -np.inf], 0, inplace=True)
 
     return pd.concat([demand, demand_non_eu28])
 
@@ -235,7 +239,14 @@ def industrial_energy_demand(countries, year):
     year = params.get("reference_year", 2015)
     countries = pd.Index(snakemake.params.countries)
 
-    demand = industrial_energy_demand(countries.intersection(eu28), year)
+    if len(countries.intersection(eu28)) > 0:
+        demand = industrial_energy_demand(countries.intersection(eu28), year)
+    else:
+        # e.g. only UA or MD
+        logger.info(
+            f"No industrial energy demand available for {countries}. Filling with average values of EU."
+        )
+        demand = industrial_energy_demand(eu28, year)
 
     # output in MtMaterial/a
     production = (

scripts/build_industrial_energy_demand_per_node_today.py

@@ -72,7 +72,15 @@ def build_nodal_industrial_energy_demand():
         buses = keys.index[keys.country == country]
         mapping = sector_mapping.get(sector, "population")
 
-        key = keys.loc[buses, mapping]
+        try:
+            key = keys.loc[buses, mapping].fillna(0)
+        except:
+            logger.info(
+                f"No industrial demand available for {mapping}. Filling with zeros."
+            )
+            keys[mapping] = 0
+            key = keys.loc[buses, mapping].fillna(0)
+
         demand = industrial_demand[country, sector]
 
         outer = pd.DataFrame(

scripts/build_industrial_production_per_node.py

@@ -75,6 +75,15 @@ def build_nodal_industrial_production():
         buses = keys.index[keys.country == country]
         mapping = sector_mapping.get(sector, "population")
 
+        try:
+            key = keys.loc[buses, mapping].fillna(0)
+        except:
+            logger.info(
+                f"No industrial production available for {mapping}. Filling with zeros."
+            )
+            keys[mapping] = 0
+            key = keys.loc[buses, mapping].fillna(0)
+
         key = keys.loc[buses, mapping]
         nodal_production.loc[buses, sector] = (
             industrial_production.at[country, sector] * key