Merge pull request #151 from PyPSA/enspreso-biomass

spatially-explicit biomass potentials from ENSPRESO

.gitignore

@@ -28,7 +28,8 @@ gurobi.log
 /data/.nfs*
 /data/Industrial_Database.csv
 /data/retro/tabula-calculator-calcsetbuilding.csv
-/data
+/data/nuts*
+
 *.org
 
 *.nc

Snakefile

@@ -24,7 +24,6 @@ subworkflow pypsaeur:
     snakefile: "../pypsa-eur/Snakefile"
     configfile: "../pypsa-eur/config.yaml"
 
-
 rule all:
     input: SDIR + '/graphs/costs.pdf'
 
@@ -173,13 +172,19 @@ rule build_energy_totals:
 
 rule build_biomass_potentials:
     input:
-        jrc_potentials="data/biomass/JRC Biomass Potentials.xlsx"
+        enspreso_biomass=HTTP.remote("https://cidportal.jrc.ec.europa.eu/ftp/jrc-opendata/ENSPRESO/ENSPRESO_BIOMASS.xlsx", keep_local=True),
+        nuts2="data/nuts/NUTS_RG_10M_2013_4326_LEVL_2.geojson", # https://gisco-services.ec.europa.eu/distribution/v2/nuts/download/#nuts21
+        regions_onshore=pypsaeur("resources/regions_onshore_elec_s{simpl}_{clusters}.geojson"),
+        nuts3_population="../pypsa-eur/data/bundle/nama_10r_3popgdp.tsv.gz",
+        swiss_cantons="../pypsa-eur/data/bundle/ch_cantons.csv",
+        swiss_population="../pypsa-eur/data/bundle/je-e-21.03.02.xls",
+        country_shapes=pypsaeur('resources/country_shapes.geojson')
     output:
-        biomass_potentials_all='resources/biomass_potentials_all.csv',
+        biomass_potentials_all='resources/biomass_potentials_all_s{simpl}_{clusters}.csv',
-        biomass_potentials='resources/biomass_potentials.csv'
+        biomass_potentials='resources/biomass_potentials_s{simpl}_{clusters}.csv'
     threads: 1
     resources: mem_mb=1000
-    benchmark: "benchmarks/build_biomass_potentials"
+    benchmark: "benchmarks/build_biomass_potentials_s{simpl}_{clusters}"
     script: 'scripts/build_biomass_potentials.py'
 
 
@@ -339,9 +344,9 @@ rule prepare_sector_network:
         energy_totals_name='resources/energy_totals.csv',
         co2_totals_name='resources/co2_totals.csv',
         transport_name='resources/transport_data.csv',
-        traffic_data_KFZ="data/emobility/KFZ__count",
+        traffic_data_KFZ = "data/emobility/KFZ__count",
-        traffic_data_Pkw="data/emobility/Pkw__count",
+        traffic_data_Pkw = "data/emobility/Pkw__count",
-        biomass_potentials='resources/biomass_potentials.csv',
+        biomass_potentials='resources/biomass_potentials_s{simpl}_{clusters}.csv',
         heat_profile="data/heat_load_profile_BDEW.csv",
         costs=CDIR + "costs_{planning_horizons}.csv",
         profile_offwind_ac=pypsaeur("resources/profile_offwind-ac.nc"),

config.default.yaml

@@ -99,28 +99,28 @@ energy:
 
 biomass:
   year: 2030
-  scenario: Med
+  scenario: ENS_Med
   classes:
     solid biomass:
-      - Primary agricultural residues
+      - Agricultural waste
-      - Forestry energy residue
+      - Fuelwood residues
-      - Secondary forestry residues
+      - Secondary Forestry residues - woodchips
-      - Secondary Forestry residues sawdust
+      - Sawdust
-      - Forestry residues from landscape care biomass
+      - Residues from landscape care
       - Municipal waste
     not included:
-      - Bioethanol sugar beet biomass
+      - Sugar from sugar beet
-      - Rapeseeds for biodiesel
+      - Rape seed
-      - sunflower and soya for Biodiesel
+      - "Sunflower, soya seed "
-      - Starchy crops biomass
+      - Bioethanol barley, wheat, grain maize, oats, other cereals and rye
-      - Grassy crops biomass
+      - Miscanthus, switchgrass, RCG
-      - Willow biomass
+      - Willow
-      - Poplar biomass potential
+      - Poplar
-      - Roundwood fuelwood
+      - FuelwoodRW
-      - Roundwood Chips & Pellets
+      - C&P_RW
     biogas:
-      - Manure biomass potential
+      - Manure solid, liquid
-      - Sludge biomass
+      - Sludge
 
 
 solar_thermal:

doc/data.csv

@@ -2,6 +2,7 @@ description,file/folder,licence,source
 JRC IDEES database,jrc-idees-2015/,CC BY 4.0,https://ec.europa.eu/jrc/en/potencia/jrc-idees
 urban/rural fraction,urban_percent.csv,unknown,unknown
 JRC biomass potentials,biomass/,unknown,https://doi.org/10.2790/39014
+JRC ENSPRESO biomass potentials,remote,CC BY 4.0,https://data.jrc.ec.europa.eu/dataset/74ed5a04-7d74-4807-9eab-b94774309d9f
 EEA emission statistics,eea/UNFCCC_v23.csv,EEA standard re-use policy,https://www.eea.europa.eu/data-and-maps/data/national-emissions-reported-to-the-unfccc-and-to-the-eu-greenhouse-gas-monitoring-mechanism-16
 Eurostat Energy Balances,eurostat-energy_balances-*/,Eurostat,https://ec.europa.eu/eurostat/web/energy/data/energy-balances
 Swiss energy statistics from Swiss Federal Office of Energy,switzerland-sfoe/,unknown,http://www.bfe.admin.ch/themen/00526/00541/00542/02167/index.html?dossier_id=02169

doc/release_notes.rst

@@ -82,6 +82,9 @@ Future release
   The transport costs are determined based on the `JRC-EU-Times Bioenergy report <http://dx.doi.org/10.2790/01017>`_
   in the new optional rule ``build_biomass_transport_costs``.
   Biomass transport can be activated with the setting ``sector: biomass_transport: true``.
+* Use `JRC ENSPRESO database <https://data.jrc.ec.europa.eu/dataset/74ed5a04-7d74-4807-9eab-b94774309d9f>`_ to
+  spatially disaggregate biomass potentials to PyPSA-Eur regions based on overlaps with NUTS2 regions from ENSPRESO
+  (proportional to area) (`#151 <https://github.com/PyPSA/pypsa-eur-sec/pull/151>`_).
 * Compatibility with ``xarray`` version 0.19.
 * Separate basic chemicals into HVC, chlorine, methanol and ammonia [`#166 <https://github.com/PyPSA/PyPSA-Eur-Sec/pull/166>`_].
 * Add option to specify reuse, primary production, and mechanical and chemical recycling fraction of platics [`#166 <https://github.com/PyPSA/PyPSA-Eur-Sec/pull/166>`_].

doc/supply_demand.rst

@@ -183,7 +183,7 @@ Solid biomass provides process heat up to 500 Celsius in industry, as well as fe
 Solid biomass supply
 =====================
 
-Only wastes and residues from the JRC biomass dataset.
+Only wastes and residues from the JRC ENSPRESO biomass dataset.
 
 
 Oil product demand

scripts/build_biomass_potentials.py

@@ -1,55 +1,194 @@
 import pandas as pd
-
+import geopandas as gpd
-rename = {"UK" : "GB", "BH" : "BA"}
 
 
-def build_biomass_potentials():
+def build_nuts_population_data(year=2013):
 
-    config = snakemake.config['biomass']
+    pop = pd.read_csv(
-    year = config["year"]
+        snakemake.input.nuts3_population,
-    scenario = config["scenario"]
+        sep=r'\,| \t|\t',
+        engine='python',
+        na_values=[":"],
+        index_col=1
+    )[str(year)]
+    
+    # only countries
+    pop.drop("EU28", inplace=True)
 
-    df = pd.read_excel(snakemake.input.jrc_potentials,
+    # mapping from Cantons to NUTS3
-                    "Potentials (PJ)",
+    cantons = pd.read_csv(snakemake.input.swiss_cantons)
-                    index_col=[0,1])
+    cantons = cantons.set_index(cantons.HASC.str[3:]).NUTS
+    cantons = cantons.str.pad(5, side='right', fillchar='0')
 
-    df.rename(columns={"Unnamed: 18": "Municipal waste"}, inplace=True)
+    # get population by NUTS3
-    df.drop(columns="Total", inplace=True)
+    swiss = pd.read_excel(snakemake.input.swiss_population, skiprows=3, index_col=0).loc["Residents in 1000"]
-    df.replace("-", 0., inplace=True)
+    swiss = swiss.rename(cantons).filter(like="CH")
 
-    column = df.iloc[:,0]
+    # aggregate also to higher order NUTS levels
-    countries = column.where(column.str.isalpha()).pad()
+    swiss = [swiss.groupby(swiss.index.str[:i]).sum() for i in range(2, 6)]
-    countries = [rename.get(ct, ct) for ct in countries]
-    countries_i = pd.Index(countries, name='country')
-    df.set_index(countries_i, append=True, inplace=True)
 
-    df.drop(index='MS', level=0, inplace=True)
+    # merge Europe + Switzerland
+    pop = pd.DataFrame(pop.append(swiss), columns=["total"])
+    
+    # add missing manually
+    pop["AL"] = 2893
+    pop["BA"] = 3871
+    pop["RS"] = 7210
+    
+    pop["ct"] = pop.index.str[:2]
+    
+    return pop
 
-    # convert from PJ to MWh
-    df = df / 3.6 * 1e6
 
-    df.to_csv(snakemake.output.biomass_potentials_all)
+def enspreso_biomass_potentials(year=2020, scenario="ENS_Low"):
+    """
+    Loads the JRC ENSPRESO biomass potentials.
+    
+    Parameters
+    ----------
+    year : int
+        The year for which potentials are to be taken.
+        Can be {2010, 2020, 2030, 2040, 2050}.
+    scenario : str
+        The scenario. Can be {"ENS_Low", "ENS_Med", "ENS_High"}.
+    
+    Returns
+    -------
+    pd.DataFrame
+        Biomass potentials for given year and scenario
+        in TWh/a by commodity and NUTS2 region.
+    """
 
-    # solid biomass includes:
+    glossary = pd.read_excel(
-    # Primary agricultural residues (MINBIOAGRW1),
+        str(snakemake.input.enspreso_biomass),
-    # Forestry energy residue (MINBIOFRSF1),
+        sheet_name="Glossary",
-    # Secondary forestry residues (MINBIOWOOW1),
+        usecols="B:D",
-    # Secondary Forestry residues – sawdust (MINBIOWOO1a)',
+        skiprows=1,
-    # Forestry residues from landscape care biomass (MINBIOFRSF1a),
+        index_col=0
-    # Municipal waste (MINBIOMUN1)',
+    )
+    
+    df = pd.read_excel(
+        str(snakemake.input.enspreso_biomass),
+        sheet_name="ENER - NUTS2 BioCom E",
+        usecols="A:H"
+    )
 
-    # biogas includes:
+    df["group"] = df["E-Comm"].map(glossary.group)
-    # Manure biomass potential (MINBIOGAS1),
+    df["commodity"] = df["E-Comm"].map(glossary.description)
-    # Sludge biomass (MINBIOSLU1),
 
-    df = df.loc[year, scenario, :]
+    to_rename = {
+        "NUTS2 Potential available by Bio Commodity": "potential",
+        "NUST2": "NUTS2",
+    }
+    df.rename(columns=to_rename, inplace=True)
+    
+    # fill up with NUTS0 if NUTS2 is not given
+    df.NUTS2 = df.apply(lambda x: x.NUTS0 if x.NUTS2 == '-' else x.NUTS2, axis=1)
 
-    grouper = {v: k for k, vv in config["classes"].items() for v in vv}
+    # convert PJ to TWh
-    df = df.groupby(grouper, axis=1).sum()
+    df.potential /= 3.6
+    df.Unit = "TWh/a"
 
-    df.index.name = "MWh/a"
+    dff = df.query("Year == @year and Scenario == @scenario")
 
-    df.to_csv(snakemake.output.biomass_potentials)
+    bio = dff.groupby(["NUTS2", "commodity"]).potential.sum().unstack()
+    
+    # currently Serbia and Kosovo not split, so aggregate
+    bio.loc["RS"] += bio.loc["XK"]
+    bio.drop("XK", inplace=True)
+    
+    return bio
+
+
+def disaggregate_nuts0(bio):    
+    """
+    Some commodities are only given on NUTS0 level.
+    These are disaggregated here using the NUTS2
+    population as distribution key.
+    
+    Parameters
+    ----------
+    bio : pd.DataFrame
+        from enspreso_biomass_potentials()
+    
+    Returns
+    -------
+    pd.DataFrame
+    """
+    
+    pop = build_nuts_population_data()
+    
+    # get population in nuts2
+    pop_nuts2 = pop.loc[pop.index.str.len() == 4]
+    by_country = pop_nuts2.total.groupby(pop_nuts2.ct).sum()
+    pop_nuts2["fraction"] = pop_nuts2.total / pop_nuts2.ct.map(by_country)
+
+    # distribute nuts0 data to nuts2 by population
+    bio_nodal = bio.loc[pop_nuts2.ct]
+    bio_nodal.index = pop_nuts2.index
+    bio_nodal = bio_nodal.mul(pop_nuts2.fraction, axis=0)
+
+    # update inplace
+    bio.update(bio_nodal)
+    
+    return bio
+
+
+def build_nuts2_shapes():
+    """
+    - load NUTS2 geometries
+    - add RS, AL, BA country shapes (not covered in NUTS 2013)
+    - consistently name ME, MK
+    """
+
+    nuts2 = gpd.GeoDataFrame(gpd.read_file(snakemake.input.nuts2).set_index('id').geometry)
+
+    countries = gpd.read_file(snakemake.input.country_shapes).set_index('name')
+    missing = countries.loc[["AL", "RS", "BA"]]
+    nuts2.rename(index={"ME00": "ME", "MK00": "MK"}, inplace=True)
+
+    return nuts2.append(missing)
+
+
+def area(gdf):
+    """Returns area of GeoDataFrame geometries in square kilometers."""
+    return gdf.to_crs(epsg=3035).area.div(1e6)
+
+
+def convert_nuts2_to_regions(bio_nuts2, regions):
+    """
+    Converts biomass potentials given in NUTS2 to PyPSA-Eur regions based on the
+    overlay of both GeoDataFrames in proportion to the area.
+
+    Parameters
+    ----------
+    bio_nuts2 : gpd.GeoDataFrame
+        JRC ENSPRESO biomass potentials indexed by NUTS2 shapes.
+    regions : gpd.GeoDataFrame
+        PyPSA-Eur clustered onshore regions
+
+    Returns
+    -------
+    gpd.GeoDataFrame
+    """
+    
+    # calculate area of nuts2 regions
+    bio_nuts2["area_nuts2"] = area(bio_nuts2)
+
+    overlay = gpd.overlay(regions, bio_nuts2)
+
+    # calculate share of nuts2 area inside region
+    overlay["share"] = area(overlay) / overlay["area_nuts2"]
+
+    # multiply all nuts2-level values with share of nuts2 inside region
+    adjust_cols = overlay.columns.difference({"name", "area_nuts2", "geometry", "share"})
+    overlay[adjust_cols] = overlay[adjust_cols].multiply(overlay["share"], axis=0)
+
+    bio_regions = overlay.groupby("name").sum()
+
+    bio_regions.drop(["area_nuts2", "share"], axis=1, inplace=True)
+    
+    return bio_regions
 
 
 if __name__ == "__main__":
@@ -57,12 +196,28 @@ if __name__ == "__main__":
         from helper import mock_snakemake
         snakemake = mock_snakemake('build_biomass_potentials')
 
+    config = snakemake.config['biomass']
+    year = config["year"]
+    scenario = config["scenario"]
 
-    # This is a hack, to be replaced once snakemake is unicode-conform
+    enspreso = enspreso_biomass_potentials(year, scenario)
 
-    solid_biomass = snakemake.config['biomass']['classes']['solid biomass']
+    enspreso = disaggregate_nuts0(enspreso)
-    if 'Secondary Forestry residues sawdust' in solid_biomass:
-        solid_biomass.remove('Secondary Forestry residues sawdust')
-        solid_biomass.append('Secondary Forestry residues – sawdust')
 
-    build_biomass_potentials()
+    nuts2 = build_nuts2_shapes()
+
+    df_nuts2 = gpd.GeoDataFrame(nuts2.geometry).join(enspreso)
+
+    regions = gpd.read_file(snakemake.input.regions_onshore)
+
+    df = convert_nuts2_to_regions(df_nuts2, regions)
+
+    df.to_csv(snakemake.output.biomass_potentials_all)
+
+    grouper = {v: k for k, vv in config["classes"].items() for v in vv}
+    df = df.groupby(grouper, axis=1).sum()
+
+    df *= 1e6 # TWh/a to MWh/a
+    df.index.name = "MWh/a"
+
+    df.to_csv(snakemake.output.biomass_potentials)

scripts/prepare_sector_network.py

@@ -26,7 +26,7 @@ spatial = SimpleNamespace()
 def define_spatial(nodes):
     """
     Namespace for spatial
-    
+
     Parameters
     ----------
     nodes : list-like
@@ -36,7 +36,7 @@ def define_spatial(nodes):
     global options
 
     spatial.nodes = nodes
-    
+
     # biomass
 
     spatial.biomass = SimpleNamespace()
@@ -53,11 +53,11 @@ def define_spatial(nodes):
         spatial.biomass.industry_cc = ["solid biomass for industry CC"]
 
     spatial.biomass.df = pd.DataFrame(vars(spatial.biomass), index=nodes)
-    
+
     # co2
 
     spatial.co2 = SimpleNamespace()
-    
+
     if options["co2_network"]:
         spatial.co2.nodes = nodes + " co2 stored"
         spatial.co2.locations = nodes
@@ -72,7 +72,7 @@ def define_spatial(nodes):
 
 def emission_sectors_from_opts(opts):
 
-    sectors = ["electricity"]   
+    sectors = ["electricity"]
     if "T" in opts:
         sectors += [
             "rail non-elec",
@@ -107,13 +107,13 @@ def get(item, investment_year=None):
 def create_network_topology(n, prefix, connector=" -> "):
     """
     Create a network topology like the power transmission network.
-    
+
     Parameters
     ----------
     n : pypsa.Network
     prefix : str
     connector : str
-        
+
     Returns
     -------
     pd.DataFrame with columns bus0, bus1 and length
@@ -228,7 +228,7 @@ def add_lifetime_wind_solar(n, costs):
 def create_network_topology(n, prefix, connector=" -> ", bidirectional=True):
     """
     Create a network topology like the power transmission network.
-    
+
     Parameters
     ----------
     n : pypsa.Network
@@ -237,7 +237,7 @@ def create_network_topology(n, prefix, connector=" -> ", bidirectional=True):
     bidirectional : bool, default True
         True: one link for each connection
         False: one link for each connection and direction (back and forth)
-        
+
     Returns
     -------
     pd.DataFrame with columns bus0, bus1 and length
@@ -256,19 +256,19 @@ def create_network_topology(n, prefix, connector=" -> ", bidirectional=True):
     swap_buses = {"bus0": "bus1", "bus1": "bus0"}
     candidates_n = candidates[~positive_order].rename(columns=swap_buses)
     candidates = pd.concat([candidates_p, candidates_n])
-    
+
     def make_index(c):
         return prefix + c.bus0 + connector + c.bus1
 
     topo = candidates.groupby(["bus0", "bus1"], as_index=False).mean()
     topo.index = topo.apply(make_index, axis=1)
-    
+
     if not bidirectional:
         topo_reverse = topo.copy()
         topo_reverse.rename(columns=swap_buses, inplace=True)
         topo_reverse.index = topo_reverse.apply(make_index, axis=1)
         topo = topo.append(topo_reverse)
-    
+
     return topo
 
 
@@ -439,7 +439,7 @@ def add_co2_tracking(n, options):
     n.madd("Store",
         spatial.co2.nodes,
         e_nom_extendable=True,
-        e_nom_max=np.inf, 
+        e_nom_max=np.inf,
         capital_cost=options['co2_sequestration_cost'],
         carrier="co2 stored",
         bus=spatial.co2.nodes
@@ -458,14 +458,14 @@ def add_co2_tracking(n, options):
 
 
 def add_co2_network(n, costs):
-    
+
     logger.info("Adding CO2 network.")
     co2_links = create_network_topology(n, "CO2 pipeline ")
 
     cost_onshore = (1 - co2_links.underwater_fraction) * costs.at['CO2 pipeline', 'fixed'] * co2_links.length
     cost_submarine = co2_links.underwater_fraction * costs.at['CO2 submarine pipeline', 'fixed'] * co2_links.length
     capital_cost = cost_onshore + cost_submarine
-    
+
     n.madd("Link",
         co2_links.index,
         bus0=co2_links.bus0.values + " co2 stored",
@@ -1684,17 +1684,10 @@ def add_biomass(n, costs):
 
     print("adding biomass")
 
-    # biomass distributed at country level - i.e. transport within country allowed
-    countries = n.buses.country.dropna().unique()
-
     biomass_potentials = pd.read_csv(snakemake.input.biomass_potentials, index_col=0)
 
     if options["biomass_transport"]:
-        # potential per node distributed within country by population
+        biomass_potentials_spatial = biomass_potentials.rename(index=lambda x: x + " solid biomass")
-        biomass_potentials_spatial = (biomass_potentials.loc[pop_layout.ct]
-                            .set_index(pop_layout.index)
-                            .mul(pop_layout.fraction, axis="index")
-                            .rename(index=lambda x: x + " solid biomass"))
     else:
         biomass_potentials_spatial = biomass_potentials.sum()
 
@@ -1717,9 +1710,9 @@ def add_biomass(n, costs):
         "EU biogas",
         bus="EU biogas",
         carrier="biogas",
-        e_nom=biomass_potentials.loc[countries, "biogas"].sum(),
+        e_nom=biomass_potentials["biogas"].sum(),
         marginal_cost=costs.at['biogas', 'fuel'],
-        e_initial=biomass_potentials.loc[countries, "biogas"].sum()
+        e_initial=biomass_potentials["biogas"].sum()
     )
 
     n.madd("Store",
@@ -1750,7 +1743,7 @@ def add_biomass(n, costs):
             index_col=0,
             squeeze=True
         )
-        
+
         # add biomass transport
         biomass_transport = create_network_topology(n, "biomass transport ", bidirectional=False)
 
@@ -2205,11 +2198,11 @@ if __name__ == "__main__":
         snakemake = mock_snakemake(
             'prepare_sector_network',
             simpl='',
-            clusters=45,
+            clusters="45",
             lv=1.0,
             opts='',
             sector_opts='Co2L0-168H-T-H-B-I-solar3-dist1',
-            planning_horizons=2030,
+            planning_horizons="2030",
         )
 
     logging.basicConfig(level=snakemake.config['logging_level'])