From 0a0248e0c435c2d07d4425045ca9b1a4d5724f7b Mon Sep 17 00:00:00 2001 From: BertoGBG Date: Tue, 23 Jun 2026 09:28:13 +0200 Subject: [PATCH 01/10] Add perennialisation as a carbon dioxide removal technology Ported from the a_CDRs development branch. Adds perennialisation (conversion of 1st-generation biofuel cropland to perennial grasses) as an optional CDR technology: land area displaced from 1G biofuel production is backed out from existing biomass potentials and NUTS2 crop yields, converted to a CO2 sequestration potential, and dispatched via an April-October harvest profile. Retrieves NUTS2 Eurostat crop-yield data from the same Zenodo-archived CO2-removal data package afforestation will also depend on (rules/retrieve.smk: retrieve_co2_removal_data, kept byte-identical to avoid conflicts whichever PR merges first). Off by default (sector.perennials: false). Co-Authored-By: Claude Sonnet 4.6 --- config/config.default.yaml | 12 + config/plotting.default.yaml | 1 + data/versions.csv | 1 + rules/build_sector.smk | 38 ++ rules/retrieve.smk | 34 ++ scripts/_check_utils.py | 218 +++++++++++ scripts/_helpers.py | 45 +++ scripts/build_biomass_potentials.py | 7 +- scripts/build_perennials_crop_yields_nuts2.py | 357 ++++++++++++++++++ scripts/build_perennials_potentials.py | 203 ++++++++++ scripts/check_perennials_pipeline.py | 207 ++++++++++ scripts/prepare_sector_network.py | 108 ++++++ 12 files changed, 1229 insertions(+), 2 deletions(-) create mode 100644 scripts/_check_utils.py create mode 100644 scripts/build_perennials_crop_yields_nuts2.py create mode 100644 scripts/build_perennials_potentials.py create mode 100644 scripts/check_perennials_pipeline.py diff --git a/config/config.default.yaml b/config/config.default.yaml index eee27b2fb3..05bcb35da7 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -419,6 +419,14 @@ renewable: eia_correct_by_capacity: false eia_approximate_missing: false +# alternative carbon dioxide removal technologies +perennials: + potential_co2: 2 # tCO2e/(ha y) sequestered + biofuel_conversion: # Biofuel conversion factors from JRC Technical Report doi:10.2760/69179 + MINBIOCRP11: 0.295 # t_ethanol / t_wheat grain (13.5% moisture), Table 93 + MINBIOCRP21: 0.07777 # t_ethanol / t_sugar beet (16% sugar content), Table 133 + MINBIORPS1: 0.4176 # t_crude_oil / t_rapeseed (9% moisture) x crude-to-FAME (1/1.0063), Tables 155+159 + # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#conventional conventional: unit_commitment: false @@ -836,6 +844,7 @@ sector: methanation: true coal_cc: false dac: true + perennials: false co2_vent: false heat_vent: urban central: true @@ -1309,6 +1318,9 @@ data: instrat_co2_prices: source: primary version: latest + co2_removal_data: + source: primary + version: latest co2stop: source: archive version: latest diff --git a/config/plotting.default.yaml b/config/plotting.default.yaml index a2d7405123..35bb9046a0 100644 --- a/config/plotting.default.yaml +++ b/config/plotting.default.yaml @@ -737,3 +737,4 @@ plotting: import NH3: '#e2ed74' import oil: '#93eda2' import methanol: '#87d0e6' + co2 perennials: '#008ffc' diff --git a/data/versions.csv b/data/versions.csv index c08fa55057..2e55b6e978 100644 --- a/data/versions.csv +++ b/data/versions.csv @@ -11,6 +11,7 @@ bidding_zones_electricitymaps,v1.238.0,primary,latest supported,2026-01-22,,http bidding_zones_electricitymaps,v1.238.0,archive,latest supported,2026-01-22,,https://data.pypsa.org/workflows/eur/bidding_zones_electricitymaps/v1.238.0/world.geojson bidding_zones_entsoepy,v0.6.18,primary,latest supported,2026-01-22,,https://raw.githubusercontent.com/EnergieID/entsoe-py/refs/tags/V0.6.18/entsoe/geo/geojson bidding_zones_entsoepy,v0.6.18,archive,latest supported,2026-01-22,,https://data.pypsa.org/workflows/eur/bidding_zones_entsoepy/v0.6.18 +co2_removal_data,v1.0.1,primary,latest supported,2026-06-22,"Afforestation and perennialisation input data for the CDR sector technologies, archived from https://github.com/BertoGBG/CO2-removal",https://zenodo.org/records/20799337/files/BertoGBG/CO2-removal-v1.0.1.zip co2stop,26-august-2020,primary,latest supported,2025-12-02,,https://setis.ec.europa.eu/document/download/786a884f-0b33-4789-b744-28004b16bd1a_en?filename=co2jrc_openformats.zip co2stop,26-august-2020,archive,latest supported,2026-01-13,,https://data.pypsa.org/workflows/eur/co2stop/26-august-2020/co2jrc_openformats.zip copernicus_land_cover,v3.0.1,primary,latest supported,2025-12-02,"The primary is already from Zenodo, documentation in https://zenodo.org/records/4723921",https://zenodo.org/records/3939050/files/PROBAV_LC100_global_v3.0.1_2019-nrt_Discrete-Classification-map_EPSG-4326.tif diff --git a/rules/build_sector.smk b/rules/build_sector.smk index 59b6540fea..c876b3f3f6 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -920,6 +920,39 @@ rule build_biomass_potentials: scripts("build_biomass_potentials.py") +rule build_perennials_yields_nuts_file: + params: + biofuel_conversion=config_provider("perennials", "biofuel_conversion"), + input: + nuts2021=rules.retrieve_eu_nuts_2021.output.shapes_level_2, + crops_nuts2=rules.retrieve_co2_removal_data.output.eurostat_crops_nuts2, + crops_nuts0=rules.retrieve_co2_removal_data.output.eurostat_crops_nuts0, + output: + yields_all=resources("perennials_yields_1G_biofuels.csv"), + log: + logs("build_perennials_yields_nuts_file.log"), + script: + scripts("build_perennials_crop_yields_nuts2.py") + + +rule build_perennial_potentials: + params: + biomass=config_provider("biomass"), + input: + nuts2=rules.retrieve_eu_nuts_2021.output.shapes_level_2, + country_shapes=resources("country_shapes.geojson"), + perennials_yields_1G_biofuels=resources("perennials_yields_1G_biofuels.csv"), + regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"), + output: + csv_file=resources("perennials_yields_1G_biofuels_s_{clusters}.csv"), + log: + logs("build_perennial_potentials_s_{clusters}.log"), + resources: + mem_mb=8000, + script: + scripts("build_perennials_potentials.py") + + rule build_biomass_transport_costs: input: sc1="data/biomass_transport_costs_supplychain1.csv", @@ -1689,6 +1722,11 @@ rule prepare_sector_network: if config_provider("sector", "district_heating", "ates", "enable")(w) else [] ), + perennials_yields_1G_biofuels=lambda w: ( + resources("perennials_yields_1G_biofuels_s_{clusters}.csv") + if config_provider("sector", "perennials")(w) + else [] + ), output: resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc" diff --git a/rules/retrieve.smk b/rules/retrieve.smk index cde7f32151..d95cfc7a71 100755 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -1689,3 +1689,37 @@ if (MOBILITY_PROFILES_DATASET := dataset_version("mobility_profiles"))["source"] run: copy2(input["kfz"], output["kfz"]) copy2(input["pkw"], output["pkw"]) + + +if (CO2_REMOVAL_DATASET := dataset_version("co2_removal_data"))["source"] in [ + "primary", + "archive", +]: + + rule retrieve_co2_removal_data: + message: + "Downloading carbon dioxide removal data (afforestation, perennialisation inputs)" + input: + zip=storage(CO2_REMOVAL_DATASET["url"]), + output: + afforestation_nuts_biomass_densities=resources("afforestation_nuts_biomass_densities.xlsx"), + afforestation_nuts2_afforestation_rates=resources("afforestation_rates_nuts2_full.csv"), + afforestation_nuts2_monthly_weights=resources("afforestation_nuts2_monthly_weights.csv"), + eurostat_crops_nuts2=resources("eurostat_apro_cpshr_nuts2_raw.csv"), + eurostat_crops_nuts0=resources("eurostat_apro_cpshr_nuts0_raw.csv"), + retries: 2 + run: + with ZipFile(input.zip) as z: + # GitHub's release archive nests everything under a single + # top-level "--/" folder whose name + # changes with every release, so resolve it at runtime. + top_dir = z.namelist()[0].split("/")[0] + for src_path, dest in [ + ("outputs/afforestation/afforestation_nuts_biomass_densities.xlsx", output.afforestation_nuts_biomass_densities), + ("outputs/afforestation/afforestation_rates_nuts2_full.csv", output.afforestation_nuts2_afforestation_rates), + ("outputs/afforestation/afforestation_nuts2_monthly_weights.csv", output.afforestation_nuts2_monthly_weights), + ("outputs/perennialisation/eurostat_apro_cpshr_nuts2_raw.csv", output.eurostat_crops_nuts2), + ("outputs/perennialisation/eurostat_apro_cpshr_nuts0_raw.csv", output.eurostat_crops_nuts0), + ]: + with z.open(f"{top_dir}/{src_path}") as src, open(dest, "wb") as dst: + dst.write(src.read()) diff --git a/scripts/_check_utils.py b/scripts/_check_utils.py new file mode 100644 index 0000000000..2dc9a6f1a0 --- /dev/null +++ b/scripts/_check_utils.py @@ -0,0 +1,218 @@ +""" +Shared utilities for carbon dioxide removal pipeline check scripts. + +Loads run parameters from the saved pypsa-eur config so check scripts need +only a single input: the run name. + +Usage in any check script: + + from _check_utils import parse_check_args, load_check_params + + args = parse_check_args() # positional run_name, optional overrides + p = load_check_params(args) # returns dict with RDIR, WC, paths, cfg, … + +Run from the pypsa-eur root: + + python scripts/check_EW_pipeline.py EW_2050 + python scripts/check_EW_pipeline.py --config results/EW_2050/EW_2050/configs/config.base_s_90__168h_2050.yaml +""" + +import argparse +import sys +from pathlib import Path + +import yaml + +BASE_DIR = Path(".") # check scripts are run from the pypsa-eur root + + +# ── Saved-config discovery ───────────────────────────────────────────────────── + +def find_saved_config(run_name: str, base_dir: Path) -> Path: + """Find the config saved by pypsa-eur under results//**/configs/.""" + pattern = f"results/{run_name}/**/configs/config.*.yaml" + hits = sorted(base_dir.glob(pattern)) + if not hits: + sys.exit( + f"ERROR: no saved config found for run '{run_name}'.\n" + f" Looked for: {base_dir / pattern}\n" + f" Use --config to point directly at the config file." + ) + if len(hits) > 1: + print( + "WARNING: multiple saved configs found; using first:\n " + + "\n ".join(str(h) for h in hits) + ) + return hits[0] + + +def _load_yaml(path: Path) -> dict: + with open(path) as f: + return yaml.safe_load(f) or {} + + +def load_config(args: argparse.Namespace, base_dir: Path = BASE_DIR) -> tuple: + """ + Load the saved run config. Returns (cfg_dict, config_path). + + Priority: + 1. --config → load that file directly + 2. run_name (positional) → auto-discover under results//**/configs/ + """ + if getattr(args, "config", None): + config_path = Path(args.config) + if not config_path.exists(): + sys.exit(f"ERROR: config file not found: {config_path}") + elif getattr(args, "run_name", None): + config_path = find_saved_config(args.run_name, base_dir) + else: + sys.exit("ERROR: provide a run_name or --config .") + + cfg = _load_yaml(config_path) + print(f"Config: {config_path}") + return cfg, config_path + + +# ── CLI ──────────────────────────────────────────────────────────────────────── + +def parse_check_args(extra_args=None) -> argparse.Namespace: + """ + Common CLI parser for all carbon dioxide removal check scripts. + + Positional: + run_name Run directory name (e.g. rock_weathering_2050). Script finds the saved + config under results//**/configs/ automatically. + + Optional overrides (take precedence over config values): + --config Direct path to a saved config YAML (skips auto-discovery). + --base-dir pypsa-eur root directory (default: current directory). + --run-name Override run name from config. + --clusters Override cluster count. + --opts Override opts wildcard. + --sector-opts Override sector opts wildcard. + --horizon Override planning horizon year. + + extra_args: list of ([flags], kwargs) for script-specific arguments. + """ + p = argparse.ArgumentParser( + description="Carbon dioxide removal pipeline check — reads wildcards from the saved run config.", + formatter_class=argparse.ArgumentDefaultsHelpFormatter, + ) + p.add_argument( + "run_name", nargs="?", + help="Run name (e.g. EW_2050). Auto-discovers saved config under results//.", + ) + p.add_argument( + "--config", "-c", default=None, metavar="PATH", + help="Direct path to a saved config YAML (overrides run_name auto-discovery).", + ) + p.add_argument( + "--base-dir", default=".", metavar="DIR", + help="pypsa-eur root directory (default: current directory).", + ) + p.add_argument("--run-name", default=None, metavar="NAME", + help="Override run name (= results/resources sub-directory).") + p.add_argument("--clusters", default=None, metavar="N", + help="Override cluster count (e.g. 90).") + p.add_argument("--opts", default=None, metavar="OPTS", + help="Override opts wildcard.") + p.add_argument("--sector-opts", default=None, metavar="OPTS", + help="Override sector opts wildcard (e.g. 168h).") + p.add_argument("--horizon", default=None, metavar="YEAR", + help="Override planning horizon year (e.g. 2050).") + + if extra_args: + for flags, kwargs in extra_args: + p.add_argument(*flags, **kwargs) + + args = p.parse_args() + if not args.run_name and not args.config: + p.print_help() + sys.exit(1) + return args + + +# ── Parameter extraction ─────────────────────────────────────────────────────── + +def load_check_params(args: argparse.Namespace) -> dict: + """ + Load the saved config and return a dict with all derived parameters: + + RDIR – run directory name + CLUSTERS – cluster count string (e.g. "90") + OPTS – opts wildcard (e.g. "") + SECTOR_OPTS – sector opts wildcard (e.g. "168h") + PLANNING_HORIZON – planning horizon string (e.g. "2050") + WC – full wildcard stem (e.g. "base_s_90__168h_2050") + SHARED_RES – shared-resources Path (= RES_RUN if no shared policy) + BASE_DIR – Path to pypsa-eur root + RES – Path to resources/ + RES_RUN – Path to resources// + RESULTS – Path to results// + cfg – full config dict (for script-specific lookups) + """ + base_dir = Path(getattr(args, "base_dir", ".")) + cfg, _ = load_config(args, base_dir) + + run_cfg = cfg.get("run", {}) + scenario_cfg = cfg.get("scenario", {}) + + # CLI --run-name > config run.name > positional run_name + RDIR = ( + getattr(args, "run_name_override", None) # --run-name flag (argparse stores as run_name_override below) + or run_cfg.get("name") + or getattr(args, "run_name", None) + or "run" + ) + # Note: argparse stores --run-name as args.run_name which clashes with the + # positional. We use dest="run_name_cli" to separate them. + # In practice, parse_check_args stores --run-name in args.run_name (the flag) + # and the positional in args.run_name too — last one wins in argparse. + # Simpler: just use config value as primary, CLI flags as overrides. + RDIR = run_cfg.get("name") or getattr(args, "run_name", None) or "run" + if getattr(args, "run_name", None) and not run_cfg.get("name"): + RDIR = args.run_name + + CLUSTERS = str( + args.clusters + or (scenario_cfg.get("clusters") or [90])[0] + ) + OPTS = ( + args.opts if args.opts is not None + else str((scenario_cfg.get("opts") or [""])[0]) + ) + SECTOR_OPTS = str( + args.sector_opts + or (scenario_cfg.get("sector_opts") or ["168h"])[0] + ) + PLANNING_HORIZON = str( + args.horizon + or (scenario_cfg.get("planning_horizons") or [2050])[-1] + ) + + WC = f"base_s_{CLUSTERS}_{OPTS}_{SECTOR_OPTS}_{PLANNING_HORIZON}" + + RES = base_dir / "resources" + RES_RUN = base_dir / "resources" / RDIR + RESULTS = base_dir / "results" / RDIR + + # shared resources: if run.shared_resources.policy is a string, resources live there + shared_policy = run_cfg.get("shared_resources", {}).get("policy", False) + SHARED_RES = (base_dir / "resources" / shared_policy) if shared_policy else RES_RUN + + print(f"Run: {RDIR} | WC: {WC}") + + return dict( + RDIR=RDIR, + CLUSTERS=CLUSTERS, + OPTS=OPTS, + SECTOR_OPTS=SECTOR_OPTS, + PLANNING_HORIZON=PLANNING_HORIZON, + WC=WC, + BASE_DIR=base_dir, + RES=RES, + RES_RUN=RES_RUN, + RESULTS=RESULTS, + SHARED_RES=SHARED_RES, + cfg=cfg, + ) diff --git a/scripts/_helpers.py b/scripts/_helpers.py index 4c74e4b718..5c670bfdde 100644 --- a/scripts/_helpers.py +++ b/scripts/_helpers.py @@ -1088,3 +1088,48 @@ def load_costs(cost_file: str) -> pd.DataFrame: """ return pd.read_csv(cost_file, index_col=0) + + +# 1G biofuel crop groups and their target biomass class names. +# In the upstream default config these sit in "not included". +# When sector: perennials=True, resolve_biomass_classes() moves them automatically. +ONE_G_BIOFUEL_CLASSES = { + "Bioethanol barley, wheat, grain maize, oats, other cereals and rye": "biofuels_1G_bioethanol_cereals", + "Sugar from sugar beet": "biofuels_1G_bioethanol_sugar", + "Rape seed": "biofuels_1G_biodiesel", + "Sunflower, soya seed ": "biofuels_1G_biodiesel", +} + + +def resolve_biomass_classes(classes, perennials_enabled): + """ + Auto-reallocate 1G biofuel crop groups from 'not included' into their + biofuels_1G_* target classes when perennials are enabled. + + Cases: + perennials=False → return classes unchanged. + perennials=True, all groups in 'not included' → move them (normal case). + perennials=True, no groups in 'not included' → raise AssertionError. + perennials=True, some groups elsewhere → move available ones, warn about rest. + """ + if not perennials_enabled: + return classes + classes = copy.deepcopy(classes) + not_incl = classes.get("not included", []) + in_ni = [g for g in ONE_G_BIOFUEL_CLASSES if g in not_incl] + elsewhere = [g for g in ONE_G_BIOFUEL_CLASSES if g not in not_incl] + if not in_ni: + raise AssertionError( + "sector: perennials=true but no 1G-biofuel groups are in " + "biomass: classes: 'not included'. Restore upstream defaults " + "so the groups can be reallocated automatically." + ) + if elsewhere: + logger.warning( + "perennials: 1G groups already allocated outside 'not included' " + f"— skipped (running perennials only for available groups): {elsewhere}" + ) + for g in in_ni: + classes["not included"].remove(g) + classes.setdefault(ONE_G_BIOFUEL_CLASSES[g], []).append(g) + return classes diff --git a/scripts/build_biomass_potentials.py b/scripts/build_biomass_potentials.py index be0d9bf3ff..70ab27d80c 100755 --- a/scripts/build_biomass_potentials.py +++ b/scripts/build_biomass_potentials.py @@ -12,7 +12,7 @@ import numpy as np import pandas as pd -from scripts._helpers import configure_logging, set_scenario_config +from scripts._helpers import configure_logging, resolve_biomass_classes, set_scenario_config logger = logging.getLogger(__name__) AVAILABLE_BIOMASS_YEARS = [2010, 2020, 2030, 2040, 2050] @@ -388,7 +388,10 @@ def add_unsustainable_potentials(df, input_eurostat): df.to_csv(snakemake.output.biomass_potentials_all) - grouper = {v: k for k, vv in params["classes"].items() for v in vv} + classes = resolve_biomass_classes( + params["classes"], snakemake.config["sector"].get("perennials", False) + ) + grouper = {v: k for k, vv in classes.items() for v in vv} df = df.T.groupby(grouper).sum().T input_eurostat = snakemake.input.eurostat diff --git a/scripts/build_perennials_crop_yields_nuts2.py b/scripts/build_perennials_crop_yields_nuts2.py new file mode 100644 index 0000000000..43d47f82ee --- /dev/null +++ b/scripts/build_perennials_crop_yields_nuts2.py @@ -0,0 +1,357 @@ +# SPDX-FileCopyrightText: Contributors to PyPSA-Eur +# +# SPDX-License-Identifier: MIT +""" +Retrieve crop harvest data from the Eurostat API (dataset ``apro_cpshr``) at +NUTS2 and NUTS0 resolution, compute area-weighted yields for 1st-generation +biofuel feedstocks and perennial grasses, and harmonize the results to the +NUTS2021 region definitions used by PyPSA-Eur. + +Outputs a single CSV with columns for each crop class (cereals, sugar beet, +rapeseed, perennials) indexed by NUTS2 region. + +Biofuel conversion efficiencies (t_biofuel / t_feedstock) are read from +``config["perennials"]["biofuel_conversion"]`` and sourced from: + + Banja et al. (2013), "Biofuels in the European Union - A general overview", + JRC Technical Report, doi:10.2760/69179, Tables 93, 133, 155, 159. +""" + +import logging +import os +from pathlib import Path + +import geopandas as gpd +import numpy as np +import pandas as pd +import requests + +logger = logging.getLogger(__name__) + +def harmonize_to_nuts2021(df, keep_col, nuts2021_n2): + """ + df : DataFrame indexed by ['geo', 'TIME_PERIOD', 'mapping'] + contains both NUTS2 and NUTS0 rows + keep_col : column to harmonize (e.g. 'weighted_YL_(t/ha)') + nuts2021_n2 : GeoDataFrame with index=NUTS2_ID and geometry + """ + + # NUTS2 target regions + target = nuts2021_n2.index.sort_values() + + # Split df + df = df.copy() + + # Identify NUTS2 and NUTS0 by index length + df["is_nuts2"] = df.index.str.len() == 4 + + # Split into NUTS2 and NUTS0 + df_nuts2 = df[df["is_nuts2"]] + df_nuts0 = df[~df["is_nuts2"]] + + # Prepare working layer for only NUTS2 + df_work = df_nuts2[[keep_col]].reindex(target) + + # Country code extraction (first 2 chars) + df_work["country"] = df_work.index.str[:2] + + # Fallback 2 — Use NUTS-0 values + df_work = df_work.join( + df_nuts0[[keep_col]].rename(columns={keep_col: "fallback"}), on="country" + ) + df_work[keep_col] = df_work[keep_col].fillna(df_work["fallback"]) + df_work.drop(columns=["fallback"], inplace=True) + + # Join geometry + nuts_proj = nuts2021_n2.to_crs(epsg=3035) + gdf = nuts_proj.join(df_work)[[keep_col, "country", "geometry"]] + + # Fallback 3 — Spatial neighbors mean or nearest region if island + mask_missing = gdf[keep_col].isna() | (gdf[keep_col] <= 0) + if mask_missing.any(): + logger.warning("Spatial fallback required for %d regions", mask_missing.sum()) + + # Pre-calc distance matrix only once + valid = gdf[gdf[keep_col].notna() & (gdf[keep_col] > 0)] + + for idx in gdf[mask_missing].index: + region = gdf.loc[idx, "geometry"] + + # Touching neighbors + neigh_idxs = gdf[gdf.geometry.touches(region)].index.tolist() + neigh_vals = gdf.loc[neigh_idxs, keep_col].dropna() + neigh_vals = neigh_vals[neigh_vals > 0] + + if len(neigh_vals) > 0: + gdf.at[idx, keep_col] = neigh_vals.mean() + else: + # Island fallback: nearest valid NUTS2 region + nearest_idx = valid.distance(region).idxmin() + gdf.at[idx, keep_col] = valid.at[nearest_idx, keep_col] + logger.debug("Spatial fallback for %s: nearest = %s", idx, nearest_idx) + + # Final output tidy + result = gdf[[keep_col]] + result.index.name = "NUTS2" + result.sort_index() + return result + +def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, biofuel_yields): + # filter columns - keep only relevant + df_crops_raw_nuts2 = pd.read_csv(filepath_nuts2) + df_crops_raw_nuts2["TIME_PERIOD"] = df_crops_raw_nuts2["TIME_PERIOD"].astype(int) + + df_crops_raw_nuts0 = pd.read_csv(filepath_nuts0) + df_crops_raw_nuts0["TIME_PERIOD"] = df_crops_raw_nuts0["TIME_PERIOD"].astype(int) + + df_crops_raw = pd.concat([df_crops_raw_nuts0, df_crops_raw_nuts2], ignore_index=True) + + # drop empty and irrelevant columns + columns_to_drop = [ + "Observation value", + "OBS_FLAG", + "Observation status (Flag) V2 structure", + "CONF_STATUS", + "Confidentiality status (flag)", + "Time", + "STRUCTURE_ID", + "STRUCTURE", + "STRUCTURE_NAME", + "Geopolitical entity (reporting)", + "Time frequency", + ] + + # useful columns: + # Index(['freq', 'crops', 'Crops', 'strucpro', 'Structure of production', 'geo','TIME_PERIOD', 'OBS_VALUE'], + # freq : 'A' meaning annual + # 'crops' : code of the crops e.g. 'G0000', 'G1000', 'G2000', 'G2100', 'G2900' + # 'Crops' : name of the corp e.g. Sugar beet (excluding seed + # 'strucpro' : type of data ['AR', 'MA', 'PR_HU_EU'], where AR is cultivated area, MA is and PR_HU_EU is production at standard EU humidity + # 'OBS_VALUE' : numerical value + + df_crops = df_crops_raw.drop(columns=columns_to_drop, errors="ignore") + df_crops["OBS_VALUE"] = df_crops["OBS_VALUE"].fillna(0) + + # Step 1: Filter to relevant rows for 2023 and strucpro of interest + df_sub = df_crops[ + (df_crops["strucpro"].isin(["AR", "PR_HU_EU"])) & (df_crops["crops"].isin(crops_sel)) + ][["crops", "geo", "TIME_PERIOD", "strucpro", "OBS_VALUE"]] + + # Pivot so AR and PR_HU_EU are columns for each (crop, geo, year) + df_pivot = ( + df_sub.pivot_table( + index=["crops", "geo", "TIME_PERIOD"], + columns="strucpro", + values="OBS_VALUE", + ) + .dropna(subset=["AR", "PR_HU_EU"]) + .reset_index() + ) + + # Compute yield per year (PR_HU_EU / AR) + df_pivot["YL_(t/ha)"] = np.divide( + df_pivot["PR_HU_EU"], + df_pivot["AR"], + out=np.zeros_like(df_pivot["PR_HU_EU"], dtype=float), + where=df_pivot["AR"] != 0, + ) + + # Average data across years + df_avg_yield = ( + df_pivot.groupby(["crops", "geo"], as_index=False)[["AR", "PR_HU_EU", "YL_(t/ha)"]] + .mean() + ) + + min_year = df_pivot["TIME_PERIOD"].min() + max_year = df_pivot["TIME_PERIOD"].max() + df_avg_yield["TIME_PERIOD"] = f"{min_year}-{max_year}" + + # map crops to categories unsustainable biofuels in + rev_map = { + code: key + for key, val in crops_mapping.items() + for code in (val if isinstance(val, list) else [val]) + } + df_avg_yield["mapping"] = df_avg_yield["crops"].map(rev_map) + + # calculate weighted production per crop within mapping classes + df_avg_yield["PR_share"] = df_avg_yield["PR_HU_EU"] / df_avg_yield.groupby( + ["geo", "TIME_PERIOD", "mapping"] + )["PR_HU_EU"].transform("sum") + df_avg_yield["PR_share"] = df_avg_yield["PR_share"].fillna(0) + + # calculated average weighted yield + df_avg_yield["weighted_YL_(t/ha)"] = df_avg_yield["YL_(t/ha)"] * df_avg_yield["PR_share"] + + # sanity check for very low yields due to small productions + thresholds = { + "MINBIOCRP11": 2.0, # cereals + "MINBIOCRP21": 50.0, # sugar beet + "MINBIORPS1": 1.5, # rapeseed + "PERENNIALS": 5.0, # perennial grasses + } + + # Apply crop-specific minimum threshold + df_avg_yield["weighted_YL_(t/ha)"] = df_avg_yield.apply( + lambda row: row["weighted_YL_(t/ha)"] + if row["weighted_YL_(t/ha)"] >= thresholds.get(row["mapping"], 0) + else 0, + axis=1, + ) + + # weighted yields from current production : applies to unsustainable biofuels + weighted_yields = df_avg_yield.groupby(["geo", "TIME_PERIOD", "mapping"])[ + "weighted_YL_(t/ha)" + ].sum() + + unsustainable_biofuels_yields = pd.DataFrame(weighted_yields) + + # unsustainable biomass yield units from t/ha to MWh/ha + unsustainable_biofuels_yields = unsustainable_biofuels_yields[ + unsustainable_biofuels_yields.index.get_level_values("mapping") != "PERENNIALS" + ] + + unsustainable_biofuels_yields["energy_yields_(MWh/ha)"] = ( + unsustainable_biofuels_yields["weighted_YL_(t/ha)"] + * unsustainable_biofuels_yields.index.get_level_values("mapping").map(biofuel_yields) + ) + + # yields of perennials per hectar in ton/ha + # standard humidity for perennials = 0.65 (tH2O/t_fresh) -> note production is for fresh until 2025 + std_moist_perennials = 0.65 + + perennial_yields = pd.DataFrame(weighted_yields) + perennial_yields = perennial_yields[ + perennial_yields.index.get_level_values("mapping") == "PERENNIALS" + ] * (1 - std_moist_perennials) + + # max yields from current production : applies to perennials for green biorefining + max_yields = df_avg_yield.groupby(["geo", "TIME_PERIOD", "mapping"])["YL_(t/ha)"].max() + + perennial_yields_max = pd.DataFrame(max_yields) + perennial_yields_max = perennial_yields_max[ + perennial_yields_max.index.get_level_values("mapping") == "PERENNIALS" + ] * (1 - std_moist_perennials) + + return unsustainable_biofuels_yields, perennial_yields, perennial_yields_max + + +if __name__ == "__main__": + if "snakemake" not in globals(): + from scripts._helpers import mock_snakemake + + snakemake = mock_snakemake("build_perennials_yields_nuts_file") + + from scripts._helpers import configure_logging, set_scenario_config + + configure_logging(snakemake) + set_scenario_config(snakemake) + + CROPS_CSV_NUTS2 = Path(snakemake.input["crops_nuts2"]) + CROPS_CSV_NUTS0 = Path(snakemake.input["crops_nuts0"]) + NUTS2_2021_GEOJSON = Path(snakemake.input["nuts2021"]) + OUT_CSV_YIELDS_ALL = Path(snakemake.output["yields_all"]) + + CROPS_CSV_NUTS2.parent.mkdir(parents=True, exist_ok=True) + OUT_CSV_YIELDS_ALL.parent.mkdir(parents=True, exist_ok=True) + + perennial_codes = ["G0000", "G1000", "G2000", "G2100", "G2900"] + + crops_mapping = dict( + MINBIOCRP11=["C0000", "C1000", "C1210", "C1300", "C1310", "C1320"], + MINBIOCRP21="R2000", + MINBIORPS1=["I1110", "I1120", "I1130", "I1110-1130", "I0000"], + PERENNIALS=perennial_codes, + ) + + conv = snakemake.params.biofuel_conversion + # LHV values are fixed physical constants, not parameters: JRC Technical Report doi:10.2760/69179 + LHV_fuels = {"ethanol": 7.447, "biodiesel": 10.194} # MWh/t (26.81 MJ/kg, 36.7 MJ/kg) + + biofuel_yields = { + "MINBIOCRP11": conv["MINBIOCRP11"] * LHV_fuels["ethanol"], + "MINBIOCRP21": conv["MINBIOCRP21"] * LHV_fuels["ethanol"], + "MINBIORPS1": conv["MINBIORPS1"] * LHV_fuels["biodiesel"], + } + + other_crops_codes = [ + item + for v in crops_mapping.values() + for item in (v if isinstance(v, list) else [v]) + ] + crops_sel = perennial_codes + other_crops_codes + + logger.info("Computing crop yields...") + unsustainable_biofuels_yields, perennial_yields, perennial_yields_max = calculate_yields( + filepath_nuts0=CROPS_CSV_NUTS0, + filepath_nuts2=CROPS_CSV_NUTS2, + crops_sel=crops_sel, + crops_mapping=crops_mapping, + biofuel_yields=biofuel_yields, + ) + + yield_MINBIOCRP11 = unsustainable_biofuels_yields[ + unsustainable_biofuels_yields.index.get_level_values("mapping") == "MINBIOCRP11" + ] + yield_MINBIOCRP21 = unsustainable_biofuels_yields[ + unsustainable_biofuels_yields.index.get_level_values("mapping") == "MINBIOCRP21" + ] + yield_MINBIORPS1 = unsustainable_biofuels_yields[ + unsustainable_biofuels_yields.index.get_level_values("mapping") == "MINBIORPS1" + ] + + yield_MINBIOCRP11 = yield_MINBIOCRP11.droplevel(["TIME_PERIOD", "mapping"]) + yield_MINBIOCRP11.index.name = "NUTS2" + yield_MINBIOCRP21 = yield_MINBIOCRP21.droplevel(["TIME_PERIOD", "mapping"]) + yield_MINBIOCRP21.index.name = "NUTS2" + yield_MINBIORPS1 = yield_MINBIORPS1.droplevel(["TIME_PERIOD", "mapping"]) + yield_MINBIORPS1.index.name = "NUTS2" + perennial_yields = perennial_yields.droplevel(["TIME_PERIOD", "mapping"]) + perennial_yields.index.name = "NUTS2" + perennial_yields_max = perennial_yields_max.droplevel(["TIME_PERIOD", "mapping"]) + perennial_yields_max.index.name = "NUTS2" + + logger.info("Harmonizing to NUTS2021 regions...") + nuts2021_n2 = ( + gpd.read_file(NUTS2_2021_GEOJSON) + .loc[:, ["NUTS_ID", "NUTS_NAME", "CNTR_CODE", "geometry"]] + .set_index("NUTS_ID") + ) + + yield_MINBIOCRP11_full = harmonize_to_nuts2021( + yield_MINBIOCRP11, "energy_yields_(MWh/ha)", nuts2021_n2 + ) + yield_MINBIOCRP21_full = harmonize_to_nuts2021( + yield_MINBIOCRP21, "energy_yields_(MWh/ha)", nuts2021_n2 + ) + yield_MINBIORPS1_full = harmonize_to_nuts2021( + yield_MINBIORPS1, "energy_yields_(MWh/ha)", nuts2021_n2 + ) + yields_perennials_max_full = harmonize_to_nuts2021( + perennial_yields_max, "YL_(t/ha)", nuts2021_n2 + ) + yields_perennials_full = harmonize_to_nuts2021( + perennial_yields, "weighted_YL_(t/ha)", nuts2021_n2 + ) + + df_yields_all = pd.concat( + { + "MINBIOCRP11": yield_MINBIOCRP11_full["energy_yields_(MWh/ha)"], + "MINBIOCRP21": yield_MINBIOCRP21_full["energy_yields_(MWh/ha)"], + "MINBIORPS1": yield_MINBIORPS1_full["energy_yields_(MWh/ha)"], + "PERENNIALS_MAX": yields_perennials_max_full["YL_(t/ha)"], + }, + axis=1, + ) + df_yields_all.columns = [ + "Bioethanol barley, wheat, grain maize, oats, other cereals and rye", + "Sugar from sugar beet", + "Rape seed", + "perennials", + ] + df_yields_all = df_yields_all.sort_index() + + logger.info("Saving output CSV files...") + df_yields_all.to_csv(OUT_CSV_YIELDS_ALL, index=True) + + logger.info("Done.") diff --git a/scripts/build_perennials_potentials.py b/scripts/build_perennials_potentials.py new file mode 100644 index 0000000000..747e50ac57 --- /dev/null +++ b/scripts/build_perennials_potentials.py @@ -0,0 +1,203 @@ +# SPDX-FileCopyrightText: Contributors to PyPSA-Eur +# +# SPDX-License-Identifier: MIT +""" +Reproject NUTS2-level 1st-generation (1G) biofuel and perennial crop yields +(from ``build_perennials_crop_yields_nuts2.py``) onto the clustered network +regions, for use by ``add_perennials()`` in ``prepare_sector_network.py``. + +NUTS2 yields (MWh/ha/y for 1G biofuels, t/ha/y for perennials) are mapped to +clustered regions via an area-weighted overlay (NUTS2 geometries intersected +with cluster region geometries, weighted by intersection area); NUTS2 +regions not covered by Eurostat data (non-EU countries, small islands, +city-states) are first filled from the nearest valid NUTS2 centroid. +Resulting columns are then grouped into biomass classes +(``resolve_biomass_classes()``) to match the class structure used by +``build_biomass_potentials.py``. + +``add_perennials()`` later combines this script's clustered-region 1G yields +with the ENSPRESO 1G biomass potential (MWh/y) to back out the land area +available for conversion to perennial grasses, and converts that area into a +CO2 sequestration potential via ``perennials.potential_co2`` (tCO2/ha/y). + +Outputs a single CSV with one column per crop class (cereals, sugar beet, +rapeseed, perennials) indexed by clustered region name. +""" + +import logging + +import geopandas as gpd +import pandas as pd +from _helpers import configure_logging, resolve_biomass_classes, set_scenario_config + +logger = logging.getLogger(__name__) + + +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("NUTS_ID").geometry + ) + + countries = gpd.read_file(snakemake.input.country_shapes).set_index("name") + missing_iso2 = countries.index.intersection(["AL", "RS", "XK", "BA"]) + missing = countries.loc[missing_iso2] + + nuts2.rename(index={"ME00": "ME", "MK00": "MK"}, inplace=True) + + return pd.concat([nuts2, missing]) + + +def area(gdf): + return gdf.to_crs(epsg=3035).area.div(1e6) + + +def impute_missing_values(df_nuts2, missing_shapes, yield_cols): + + # Keep only rows that have valid yields (drop NaN rows!) + df_valid = df_nuts2.dropna(subset=yield_cols).copy() + + # Project to 3035 + valid_proj = df_valid.to_crs(3035) + missing_proj = missing_shapes.to_crs(3035) + + valid_centroids = valid_proj.centroid + missing_centroids = missing_proj.centroid + + imputed_rows = [] + + for missing_id, c_geom in missing_centroids.items(): + + # Distance to only VALID NUTS2 rows + dists = valid_centroids.distance(c_geom) + nearest_nuts2 = dists.idxmin() + + # Copy yields + attrs = df_valid.loc[nearest_nuts2, yield_cols] + + # Create new row + new_row = missing_shapes.copy().loc[[missing_id]] + for col in yield_cols: + new_row[col] = attrs[col] + + imputed_rows.append(new_row) + + return pd.concat(imputed_rows) + + +def convert_nuts2_to_regions_yields(df_nuts2, regions, yield_cols=None): + """ + Convert NUTS2-level yields (intensive) to PyPSA regions using: + + y_n = Σ_i (y_i * A_i∩n) / Σ_i A_i∩n + + Only NUTS2 rows with non-NaN yields are used. Regions with no + overlapping valid NUTS2 get NaN. + """ + + nuts = df_nuts2.copy() + regs = regions.copy() + + # Ensure ID column for NUTS2 + nuts["NUTS_ID"] = nuts.index + + # Identify yield columns if not given + if yield_cols is None: + yield_cols = nuts.columns.difference(["geometry", "NUTS_ID"]) + + # 1) Keep only NUTS2 rows that have at least one non-NaN yield + nuts_valid = nuts.dropna(subset=yield_cols, how="all") + + # 2) Reproject to equal-area CRS for areas + nuts_valid = nuts_valid.to_crs(3035) + regs = regs.to_crs(3035) + + # 3) Overlay: intersection between regions and valid NUTS2 + overlay = gpd.overlay(regs, nuts_valid, keep_geom_type=True) + + # Return empty if nothing overlaps + if overlay.empty: + return pd.DataFrame(index=regs["name"].values, columns=yield_cols, dtype=float) + + # 4) Area of intersections (m² → km²) + overlay["area_intersection"] = overlay.geometry.area + + # Optional: drop absurdly tiny slivers + overlay = overlay[overlay["area_intersection"] > 0] + + # 5) Numerators: y_i * A_i∩n + numerators = overlay[yield_cols].multiply(overlay["area_intersection"], axis=0) + + # 6) Sum per region + numer_by_region = numerators.groupby(overlay["name"]).sum(min_count=1) + denom_by_region = overlay.groupby("name")["area_intersection"].sum() + + # 7) Weighted average + yields_regions = numer_by_region.div(denom_by_region, axis=0) + + # 8) Ensure one row per region (regions with no valid overlaps → NaN) + yields_regions = yields_regions.reindex(regs["name"].values) + + return yields_regions + + +if __name__ == "__main__": + if "snakemake" not in globals(): + from _helpers import mock_snakemake + + snakemake = mock_snakemake( + "build_perennials_potentials", + clusters="39", + planning_horizons=2050, + ) + + configure_logging(snakemake) + set_scenario_config(snakemake) + + regions = gpd.read_file(snakemake.input.regions_onshore) + nuts2 = build_nuts2_shapes() + + yields = pd.read_csv(snakemake.input.perennials_yields_1G_biofuels, index_col=0) + + df_nuts2 = gpd.GeoDataFrame(nuts2.geometry).join(yields) + + # Address the missign countries (df_nuts2 contains all NUTS2 + missing shapes (with NaNs)) + missing_countries = ["AL", "RS", "BA", "XK"] + missing_shapes = df_nuts2.loc[missing_countries, ["geometry"]] + + # Impute yields for missing shapes basest of nearest valid Nuts2 + yield_cols = [ + "Bioethanol barley, wheat, grain maize, oats, other cereals and rye", + "Sugar from sugar beet", + "Rape seed", + "perennials", + ] + imputed_missing = impute_missing_values(df_nuts2, missing_shapes, yield_cols) + df_nuts2 = pd.concat([ + df_nuts2.drop(index=missing_countries), + imputed_missing + ]) + + # convert nuts2 yields to regions + df = convert_nuts2_to_regions_yields(df_nuts2, regions) + + params = snakemake.params.biomass + + classes = resolve_biomass_classes( + params["classes"], snakemake.config["sector"].get("perennials", False) + ) + grouper = {v: k for k, vv in classes.items() for v in vv} + + # keep the column 'perennials' which is otherwise dropped + for col in df.columns: + if col not in grouper: + grouper[col] = col + + df = df.T.groupby(grouper).sum().T + + df.index.name = "name" + df.to_csv(snakemake.output.csv_file) diff --git a/scripts/check_perennials_pipeline.py b/scripts/check_perennials_pipeline.py new file mode 100644 index 0000000000..05336eca46 --- /dev/null +++ b/scripts/check_perennials_pipeline.py @@ -0,0 +1,207 @@ +""" +Check script for the perennials pipeline in pypsa-eur. +Run from the pypsa-eur root directory: + + python scripts/check_perennials_pipeline.py peren_2050 + +Wildcards are read from the saved run config automatically. +Override any value on the CLI: + --run-name NAME --clusters N --horizon YEAR --sector-opts OPTS + --shared-resources NAME + --config PATH (direct path to a saved config YAML) +""" + +import sys +from pathlib import Path + +import pandas as pd + +from _check_utils import parse_check_args, load_check_params + +# ── Configuration (from config file + CLI overrides) ────────────────────────── +_args = parse_check_args(extra_args=[ + (["--shared-resources"], + {"default": None, + "metavar": "NAME", + "help": "Override shared resources directory name. " + "Reads run.shared_resources.policy from config if not set, " + "falls back to run name."}), +]) +_p = load_check_params(_args) + +BASE_DIR = _p["BASE_DIR"] +RDIR = _p["RDIR"] +CLUSTERS = _p["CLUSTERS"] +OPTS = _p["OPTS"] +SECTOR_OPTS = _p["SECTOR_OPTS"] +PLANNING_HORIZON = _p["PLANNING_HORIZON"] +WC = _p["WC"] +RES = _p["RES"] +RESULTS = _p["RESULTS"] + +_shared_policy = ( + _args.shared_resources + or _p["cfg"].get("run", {}).get("shared_resources", {}).get("policy") + or RDIR +) +SHARED_RES_POLICY = str(_shared_policy) +SHARED_RES = BASE_DIR / "resources" / SHARED_RES_POLICY + +# ── Helpers ────────────────────────────────────────────────────────────────── +OK = " [OK]" +FAIL = " [MISSING]" +WARN = " [WARN]" + +def check_file(path: Path, label: str) -> bool: + exists = path.exists() + status = OK if exists else FAIL + size = f" ({path.stat().st_size / 1e6:.1f} MB)" if exists else "" + print(f"{status} {label}{size}") + print(f" {path}") + return exists + + +def section(title: str): + print(f"\n{'='*60}") + print(f" {title}") + print('='*60) + + +# ── 1. Retrieve rule outputs ───────────────────────────────────────────────── +section("1. RETRIEVE: eurostat crops + perennial yields") + +check_file( + SHARED_RES / "eurostat_apro_cpshr_nuts2_raw.csv", + "Eurostat NUTS2 crops" +) +check_file( + SHARED_RES / "eurostat_apro_cpshr_nuts0_raw.csv", + "Eurostat NUTS0 crops" +) +check_file( + SHARED_RES / "perennials_yields_1G_biofuels.csv", + "Perennials yields (all NUTS)" +) + +# ── 2. Build rule outputs ──────────────────────────────────────────────────── +section("2. BUILD: perennial potentials (clustered)") + +yields_clustered = SHARED_RES / f"perennials_yields_1G_biofuels_s_{CLUSTERS}.csv" +if check_file(yields_clustered, f"Perennials yields clustered s_{CLUSTERS}"): + df = pd.read_csv(yields_clustered) + print(f" Rows: {len(df)} | Columns: {list(df.columns)}") + if 'perennials' in df.columns: + print(f" perennials [tDM/ha] — min: {df['perennials'].min():.2f}, mean: {df['perennials'].mean():.2f}, max: {df['perennials'].max():.2f}") + else: + print(f" {WARN} no 'perennials' column found!") + biofuels_1G_cols = [c for c in df.columns if 'biofuels_1G' in c] + if biofuels_1G_cols: + print(f" biofuels_1G columns: {biofuels_1G_cols}") + for col in biofuels_1G_cols: + print(f" {col}: mean={df[col].mean():.3f} MWh/ha") + else: + print(f" {WARN} No 'biofuels_1G_*' columns found — biomass.classes in config.default.yaml") + print(f" {WARN} must use biofuels_1G_* names, NOT 'not included', for perennials to work!") + +# ── 3. Pre-network (prepare_sector_network output) ─────────────────────────── +section("3. PRE-NETWORK: sector-coupled (prepare_sector_network)") + +prenet_path = SHARED_RES / "networks" / f"{WC}.nc" +prenet_ok = check_file(prenet_path, f"Pre-network {WC}.nc") + +if prenet_ok: + try: + import pypsa + n = pypsa.Network(str(prenet_path)) + + # Check carriers + perenn_carriers = [c for c in n.carriers.index if "perennial" in c.lower()] + print(f"\n Carriers with 'perennial': {perenn_carriers}") + + # Check links + perenn_links = n.links[n.links.carrier.str.contains("perennial", case=False, na=False)] + print(f" Links (carrier=perennial): {len(perenn_links)}") + if not perenn_links.empty: + print(perenn_links[["bus0", "bus1", "carrier", "p_nom"]].to_string(index=True)) + + # Check stores + perenn_stores = n.stores[n.stores.carrier.str.contains("perennial", case=False, na=False)] + print(f" Stores (carrier=perennial store): {len(perenn_stores)}") + if not perenn_stores.empty: + print(perenn_stores[["bus", "carrier", "e_nom_max"]].head(10).to_string(index=True)) + if "e_nom_max" in perenn_stores.columns: + finite_max = perenn_stores["e_nom_max"][perenn_stores["e_nom_max"] < 1e18] + total_cap = finite_max.sum() + print(f"\n Total store capacity (e_nom_max): {total_cap:,.0f} tCO2 ({total_cap/1e6:.3f} MtCO2)") + if total_cap == 0: + print(f" {WARN} ALL stores have e_nom_max=0!") + print(f" {WARN} This usually means biomass.classes in config.default.yaml") + print(f" {WARN} is missing biofuels_1G_* entries — check and re-run prepare_sector_network.") + + if not perenn_carriers: + print(f"\n{WARN} No perennial carriers found — add_perennials may NOT have run!") + else: + print(f"\n{OK} Perennial components found in pre-network.") + + except Exception as e: + print(f"\n{WARN} Could not load network: {e}") +else: + print(f"\n{FAIL} Pre-network missing — prepare_sector_network has not run yet.") + +# ── 4. Optimal solution ────────────────────────────────────────────────────── +section("4. OPTIMAL SOLUTION: solved network (solve_sector_network)") + +opt_path = RESULTS / "networks" / f"{WC}.nc" +opt_ok = check_file(opt_path, f"Optimal network {WC}.nc") + +if opt_ok: + try: + import pypsa + n_opt = pypsa.Network(str(opt_path)) + + perenn_links = n_opt.links[n_opt.links.carrier.str.contains("perennial", case=False, na=False)] + perenn_stores = n_opt.stores[n_opt.stores.carrier.str.contains("perennial", case=False, na=False)] + + print(f"\n Links (carrier=perennial): {len(perenn_links)}") + if not perenn_links.empty: + cols = ["carrier", "p_nom_opt"] if "p_nom_opt" in perenn_links.columns else ["carrier", "p_nom"] + print(perenn_links[cols].to_string(index=True)) + if "p_nom_opt" in perenn_links.columns: + active = perenn_links[perenn_links["p_nom_opt"] > 0] + print(f"\n Links with p_nom_opt > 0: {len(active)}") + if active.empty: + print(f"{WARN} Perennial links exist but have zero optimal capacity.") + else: + print(f"{OK} Perennial links are deployed in the optimal solution.") + + print(f"\n Stores (carrier=perennial store): {len(perenn_stores)}") + if not perenn_stores.empty: + cols = ["carrier", "e_nom_opt"] if "e_nom_opt" in perenn_stores.columns else ["carrier", "e_nom"] + print(perenn_stores[cols].to_string(index=True)) + if "e_nom_opt" in perenn_stores.columns: + active = perenn_stores[perenn_stores["e_nom_opt"] > 0] + total = perenn_stores["e_nom_opt"].sum() + print(f"\n Stores with e_nom_opt > 0: {len(active)}") + print(f" Total e_nom_opt: {total:,.0f} tCO2 ({total / 1e6:.3f} MtCO2)") + if active.empty: + print(f"{WARN} All perennial stores have e_nom_opt = 0 (not deployed).") + else: + print(f"\n Per-node e_nom_opt [tCO2] stats:") + print(f" min: {perenn_stores['e_nom_opt'].min():,.0f}") + print(f" mean: {perenn_stores['e_nom_opt'].mean():,.0f}") + print(f" max: {perenn_stores['e_nom_opt'].max():,.0f}") + + if perenn_links.empty and perenn_stores.empty: + print(f"\n{WARN} No perennial components in optimal network!") + + except Exception as e: + print(f"\n{WARN} Could not load optimal network: {e}") +else: + print(f"\n{FAIL} Optimal network missing — solve_sector_network has not run yet.") + +# ── Summary ────────────────────────────────────────────────────────────────── +section("SUMMARY") +print(f" Run: {RDIR}") +print(f" Wildcard: {WC}") +print(f" Shared resources:{SHARED_RES.resolve()}") +print(f" Results dir: {RESULTS.resolve()}") diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 1f84f038b4..f4b7b262c1 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -1247,6 +1247,111 @@ def add_dac(n, costs): ) +def add_perennials(n, costs): + """ + Add perennialisation (CDR via 1st-generation biofuel cropland conversion) + to the network as Bus, Store, and Link components. + + Perennial grasses sequester more soil carbon than the annual 1st-generation + (1G) biofuel crops (cereals, sugar beet, rapeseed) they replace. The land + area available for conversion at each node is backed out from the biomass + potential already allocated to 1G biofuels (``biomass_potentials``, + MWh/y) divided by the 1G crop yield (MWh/ha/y) at that node, giving a + displaced area in ha; multiplying by a fixed CO2 sequestration rate per + hectare (``perennials.potential_co2``) gives the store's CO2 potential. + A single "co2 perennials" Link models the harvesting process: CO2 drawn + from the atmosphere (bus0) is converted into biogas (bus3) and stored CO2 + (bus1), with capacity restricted to the April-October harvesting season + via ``p_max_pu``. + + Parameters + ---------- + n : pypsa.Network + The PyPSA network container object + costs : pd.DataFrame + Costs and parameters for different technologies. Must contain a + 'perennials gbr' entry with 'electricity-input', 'biogas-output', + 'capital_cost', 'VOM', and 'lifetime' parameters + + Returns + ------- + None + Modifies the network object in-place by adding the perennials Bus, + Store, and Link + + Notes + ----- + Reads ``snakemake.input.biomass_potentials`` and + ``snakemake.input.perennials_yields_1G_biofuels`` (NUTS2-derived crop + yields aggregated to clustered network regions, see + ``build_perennials_crop_yields_nuts2.py`` and + ``build_perennials_potentials.py``), and + ``snakemake.config["perennials"]["potential_co2"]``. + """ + + logger.info("Adding perennials.") + + # load resources + biomass_potentials = pd.read_csv(snakemake.input.biomass_potentials, index_col=0) + perennials_yields_1G_biofuels = pd.read_csv(snakemake.input.perennials_yields_1G_biofuels).set_index("name") + + # calculate CO2 sequestration per tDM perennials + perennial_CO2_seq = perennials_yields_1G_biofuels["perennials"] / snakemake.config["perennials"]["potential_co2"] # (tDM/tCO2 seq) + + # calculate perennials potential based on conversion on first generation biofuels + perennials_area_spatial = (biomass_potentials.filter(regex='biofuels_1G') / perennials_yields_1G_biofuels.filter(regex='biofuels_1G')).sum(axis=1) + # (MWh/y) / (MWh / ha / y) = (ha) returns the area used by sum of the 3 biofuels_1G classes which can be assigned for perennials + perennials_potentials_spatial = perennials_area_spatial * snakemake.config["perennials"]["potential_co2"] # (tCO2seq) = (ha) * (tCO2 seq/ha) + + nodes = pop_layout.index + n.add("Carrier", "co2 perennials") + + n.add( + "Bus", + nodes + " perennials co2 store", + location=nodes, + carrier="co2 perennials", + unit="t_co2", + ) + + # calculate biogas production based on harvesting time (in month) + df_gbr = pd.DataFrame(index=n.snapshots, columns=["harvest"]) + df_gbr["harvest"] = df_gbr.index.month.isin([4, 5, 6, 7, 8, 9, 10]).astype(int) + p_max_pu = pd.DataFrame(index=n.snapshots, columns=nodes) + for node in nodes: + p_max_pu[node] = df_gbr["harvest"] + + n.add( + "Link", + nodes, + suffix=" perennials GBR", + bus0="co2 atmosphere", + bus1=nodes + " perennials co2 store", + bus2=nodes.values, + bus3=spatial.gas.biogas, + efficiency=1, + efficiency2=-costs.at["perennials gbr", "electricity-input"] * perennial_CO2_seq, + efficiency3=costs.at["perennials gbr", "biogas-output"] * perennial_CO2_seq, + carrier="co2 perennials", + p_nom_extendable=True, + p_max_pu=p_max_pu, + capital_cost=costs.at["perennials gbr", "capital_cost"] * perennial_CO2_seq, + marginal_cost=costs.at["perennials gbr", "VOM"] * perennial_CO2_seq, + lifetime=costs.at["perennials gbr", "lifetime"], + ) + + n.add( + "Store", + nodes, + suffix=" CO2s_perennials", + bus=nodes + " perennials co2 store", + e_nom_extendable=True, + e_nom_max=perennials_potentials_spatial.values, + carrier="co2 perennials", + e_cyclic=False, + ) + + def add_co2limit(n, options, co2_totals_file, countries, nyears, limit): """ Add a global CO2 emissions constraint to the network. @@ -6506,6 +6611,9 @@ def add_import_options( if options["dac"]: add_dac(n, costs) + if options.get("perennials"): + add_perennials(n, costs) + if not options["electricity_transmission_grid"]: decentral(n) From 0e1e7beb43de4084e62d885c56d7541f69b84349 Mon Sep 17 00:00:00 2001 From: BertoGBG Date: Tue, 23 Jun 2026 10:09:45 +0200 Subject: [PATCH 02/10] Remove check_perennials_pipeline.py diagnostic script from PR --- scripts/_check_utils.py | 218 ---------------------------------------- 1 file changed, 218 deletions(-) delete mode 100644 scripts/_check_utils.py diff --git a/scripts/_check_utils.py b/scripts/_check_utils.py deleted file mode 100644 index 2dc9a6f1a0..0000000000 --- a/scripts/_check_utils.py +++ /dev/null @@ -1,218 +0,0 @@ -""" -Shared utilities for carbon dioxide removal pipeline check scripts. - -Loads run parameters from the saved pypsa-eur config so check scripts need -only a single input: the run name. - -Usage in any check script: - - from _check_utils import parse_check_args, load_check_params - - args = parse_check_args() # positional run_name, optional overrides - p = load_check_params(args) # returns dict with RDIR, WC, paths, cfg, … - -Run from the pypsa-eur root: - - python scripts/check_EW_pipeline.py EW_2050 - python scripts/check_EW_pipeline.py --config results/EW_2050/EW_2050/configs/config.base_s_90__168h_2050.yaml -""" - -import argparse -import sys -from pathlib import Path - -import yaml - -BASE_DIR = Path(".") # check scripts are run from the pypsa-eur root - - -# ── Saved-config discovery ───────────────────────────────────────────────────── - -def find_saved_config(run_name: str, base_dir: Path) -> Path: - """Find the config saved by pypsa-eur under results//**/configs/.""" - pattern = f"results/{run_name}/**/configs/config.*.yaml" - hits = sorted(base_dir.glob(pattern)) - if not hits: - sys.exit( - f"ERROR: no saved config found for run '{run_name}'.\n" - f" Looked for: {base_dir / pattern}\n" - f" Use --config to point directly at the config file." - ) - if len(hits) > 1: - print( - "WARNING: multiple saved configs found; using first:\n " - + "\n ".join(str(h) for h in hits) - ) - return hits[0] - - -def _load_yaml(path: Path) -> dict: - with open(path) as f: - return yaml.safe_load(f) or {} - - -def load_config(args: argparse.Namespace, base_dir: Path = BASE_DIR) -> tuple: - """ - Load the saved run config. Returns (cfg_dict, config_path). - - Priority: - 1. --config → load that file directly - 2. run_name (positional) → auto-discover under results//**/configs/ - """ - if getattr(args, "config", None): - config_path = Path(args.config) - if not config_path.exists(): - sys.exit(f"ERROR: config file not found: {config_path}") - elif getattr(args, "run_name", None): - config_path = find_saved_config(args.run_name, base_dir) - else: - sys.exit("ERROR: provide a run_name or --config .") - - cfg = _load_yaml(config_path) - print(f"Config: {config_path}") - return cfg, config_path - - -# ── CLI ──────────────────────────────────────────────────────────────────────── - -def parse_check_args(extra_args=None) -> argparse.Namespace: - """ - Common CLI parser for all carbon dioxide removal check scripts. - - Positional: - run_name Run directory name (e.g. rock_weathering_2050). Script finds the saved - config under results//**/configs/ automatically. - - Optional overrides (take precedence over config values): - --config Direct path to a saved config YAML (skips auto-discovery). - --base-dir pypsa-eur root directory (default: current directory). - --run-name Override run name from config. - --clusters Override cluster count. - --opts Override opts wildcard. - --sector-opts Override sector opts wildcard. - --horizon Override planning horizon year. - - extra_args: list of ([flags], kwargs) for script-specific arguments. - """ - p = argparse.ArgumentParser( - description="Carbon dioxide removal pipeline check — reads wildcards from the saved run config.", - formatter_class=argparse.ArgumentDefaultsHelpFormatter, - ) - p.add_argument( - "run_name", nargs="?", - help="Run name (e.g. EW_2050). Auto-discovers saved config under results//.", - ) - p.add_argument( - "--config", "-c", default=None, metavar="PATH", - help="Direct path to a saved config YAML (overrides run_name auto-discovery).", - ) - p.add_argument( - "--base-dir", default=".", metavar="DIR", - help="pypsa-eur root directory (default: current directory).", - ) - p.add_argument("--run-name", default=None, metavar="NAME", - help="Override run name (= results/resources sub-directory).") - p.add_argument("--clusters", default=None, metavar="N", - help="Override cluster count (e.g. 90).") - p.add_argument("--opts", default=None, metavar="OPTS", - help="Override opts wildcard.") - p.add_argument("--sector-opts", default=None, metavar="OPTS", - help="Override sector opts wildcard (e.g. 168h).") - p.add_argument("--horizon", default=None, metavar="YEAR", - help="Override planning horizon year (e.g. 2050).") - - if extra_args: - for flags, kwargs in extra_args: - p.add_argument(*flags, **kwargs) - - args = p.parse_args() - if not args.run_name and not args.config: - p.print_help() - sys.exit(1) - return args - - -# ── Parameter extraction ─────────────────────────────────────────────────────── - -def load_check_params(args: argparse.Namespace) -> dict: - """ - Load the saved config and return a dict with all derived parameters: - - RDIR – run directory name - CLUSTERS – cluster count string (e.g. "90") - OPTS – opts wildcard (e.g. "") - SECTOR_OPTS – sector opts wildcard (e.g. "168h") - PLANNING_HORIZON – planning horizon string (e.g. "2050") - WC – full wildcard stem (e.g. "base_s_90__168h_2050") - SHARED_RES – shared-resources Path (= RES_RUN if no shared policy) - BASE_DIR – Path to pypsa-eur root - RES – Path to resources/ - RES_RUN – Path to resources// - RESULTS – Path to results// - cfg – full config dict (for script-specific lookups) - """ - base_dir = Path(getattr(args, "base_dir", ".")) - cfg, _ = load_config(args, base_dir) - - run_cfg = cfg.get("run", {}) - scenario_cfg = cfg.get("scenario", {}) - - # CLI --run-name > config run.name > positional run_name - RDIR = ( - getattr(args, "run_name_override", None) # --run-name flag (argparse stores as run_name_override below) - or run_cfg.get("name") - or getattr(args, "run_name", None) - or "run" - ) - # Note: argparse stores --run-name as args.run_name which clashes with the - # positional. We use dest="run_name_cli" to separate them. - # In practice, parse_check_args stores --run-name in args.run_name (the flag) - # and the positional in args.run_name too — last one wins in argparse. - # Simpler: just use config value as primary, CLI flags as overrides. - RDIR = run_cfg.get("name") or getattr(args, "run_name", None) or "run" - if getattr(args, "run_name", None) and not run_cfg.get("name"): - RDIR = args.run_name - - CLUSTERS = str( - args.clusters - or (scenario_cfg.get("clusters") or [90])[0] - ) - OPTS = ( - args.opts if args.opts is not None - else str((scenario_cfg.get("opts") or [""])[0]) - ) - SECTOR_OPTS = str( - args.sector_opts - or (scenario_cfg.get("sector_opts") or ["168h"])[0] - ) - PLANNING_HORIZON = str( - args.horizon - or (scenario_cfg.get("planning_horizons") or [2050])[-1] - ) - - WC = f"base_s_{CLUSTERS}_{OPTS}_{SECTOR_OPTS}_{PLANNING_HORIZON}" - - RES = base_dir / "resources" - RES_RUN = base_dir / "resources" / RDIR - RESULTS = base_dir / "results" / RDIR - - # shared resources: if run.shared_resources.policy is a string, resources live there - shared_policy = run_cfg.get("shared_resources", {}).get("policy", False) - SHARED_RES = (base_dir / "resources" / shared_policy) if shared_policy else RES_RUN - - print(f"Run: {RDIR} | WC: {WC}") - - return dict( - RDIR=RDIR, - CLUSTERS=CLUSTERS, - OPTS=OPTS, - SECTOR_OPTS=SECTOR_OPTS, - PLANNING_HORIZON=PLANNING_HORIZON, - WC=WC, - BASE_DIR=base_dir, - RES=RES, - RES_RUN=RES_RUN, - RESULTS=RESULTS, - SHARED_RES=SHARED_RES, - cfg=cfg, - ) From 569568431b83a61db5989fd56b45b525484632b9 Mon Sep 17 00:00:00 2001 From: BertoGBG Date: Tue, 23 Jun 2026 10:11:02 +0200 Subject: [PATCH 03/10] delete check_pipeline_files --- scripts/check_perennials_pipeline.py | 207 --------------------------- 1 file changed, 207 deletions(-) delete mode 100644 scripts/check_perennials_pipeline.py diff --git a/scripts/check_perennials_pipeline.py b/scripts/check_perennials_pipeline.py deleted file mode 100644 index 05336eca46..0000000000 --- a/scripts/check_perennials_pipeline.py +++ /dev/null @@ -1,207 +0,0 @@ -""" -Check script for the perennials pipeline in pypsa-eur. -Run from the pypsa-eur root directory: - - python scripts/check_perennials_pipeline.py peren_2050 - -Wildcards are read from the saved run config automatically. -Override any value on the CLI: - --run-name NAME --clusters N --horizon YEAR --sector-opts OPTS - --shared-resources NAME - --config PATH (direct path to a saved config YAML) -""" - -import sys -from pathlib import Path - -import pandas as pd - -from _check_utils import parse_check_args, load_check_params - -# ── Configuration (from config file + CLI overrides) ────────────────────────── -_args = parse_check_args(extra_args=[ - (["--shared-resources"], - {"default": None, - "metavar": "NAME", - "help": "Override shared resources directory name. " - "Reads run.shared_resources.policy from config if not set, " - "falls back to run name."}), -]) -_p = load_check_params(_args) - -BASE_DIR = _p["BASE_DIR"] -RDIR = _p["RDIR"] -CLUSTERS = _p["CLUSTERS"] -OPTS = _p["OPTS"] -SECTOR_OPTS = _p["SECTOR_OPTS"] -PLANNING_HORIZON = _p["PLANNING_HORIZON"] -WC = _p["WC"] -RES = _p["RES"] -RESULTS = _p["RESULTS"] - -_shared_policy = ( - _args.shared_resources - or _p["cfg"].get("run", {}).get("shared_resources", {}).get("policy") - or RDIR -) -SHARED_RES_POLICY = str(_shared_policy) -SHARED_RES = BASE_DIR / "resources" / SHARED_RES_POLICY - -# ── Helpers ────────────────────────────────────────────────────────────────── -OK = " [OK]" -FAIL = " [MISSING]" -WARN = " [WARN]" - -def check_file(path: Path, label: str) -> bool: - exists = path.exists() - status = OK if exists else FAIL - size = f" ({path.stat().st_size / 1e6:.1f} MB)" if exists else "" - print(f"{status} {label}{size}") - print(f" {path}") - return exists - - -def section(title: str): - print(f"\n{'='*60}") - print(f" {title}") - print('='*60) - - -# ── 1. Retrieve rule outputs ───────────────────────────────────────────────── -section("1. RETRIEVE: eurostat crops + perennial yields") - -check_file( - SHARED_RES / "eurostat_apro_cpshr_nuts2_raw.csv", - "Eurostat NUTS2 crops" -) -check_file( - SHARED_RES / "eurostat_apro_cpshr_nuts0_raw.csv", - "Eurostat NUTS0 crops" -) -check_file( - SHARED_RES / "perennials_yields_1G_biofuels.csv", - "Perennials yields (all NUTS)" -) - -# ── 2. Build rule outputs ──────────────────────────────────────────────────── -section("2. BUILD: perennial potentials (clustered)") - -yields_clustered = SHARED_RES / f"perennials_yields_1G_biofuels_s_{CLUSTERS}.csv" -if check_file(yields_clustered, f"Perennials yields clustered s_{CLUSTERS}"): - df = pd.read_csv(yields_clustered) - print(f" Rows: {len(df)} | Columns: {list(df.columns)}") - if 'perennials' in df.columns: - print(f" perennials [tDM/ha] — min: {df['perennials'].min():.2f}, mean: {df['perennials'].mean():.2f}, max: {df['perennials'].max():.2f}") - else: - print(f" {WARN} no 'perennials' column found!") - biofuels_1G_cols = [c for c in df.columns if 'biofuels_1G' in c] - if biofuels_1G_cols: - print(f" biofuels_1G columns: {biofuels_1G_cols}") - for col in biofuels_1G_cols: - print(f" {col}: mean={df[col].mean():.3f} MWh/ha") - else: - print(f" {WARN} No 'biofuels_1G_*' columns found — biomass.classes in config.default.yaml") - print(f" {WARN} must use biofuels_1G_* names, NOT 'not included', for perennials to work!") - -# ── 3. Pre-network (prepare_sector_network output) ─────────────────────────── -section("3. PRE-NETWORK: sector-coupled (prepare_sector_network)") - -prenet_path = SHARED_RES / "networks" / f"{WC}.nc" -prenet_ok = check_file(prenet_path, f"Pre-network {WC}.nc") - -if prenet_ok: - try: - import pypsa - n = pypsa.Network(str(prenet_path)) - - # Check carriers - perenn_carriers = [c for c in n.carriers.index if "perennial" in c.lower()] - print(f"\n Carriers with 'perennial': {perenn_carriers}") - - # Check links - perenn_links = n.links[n.links.carrier.str.contains("perennial", case=False, na=False)] - print(f" Links (carrier=perennial): {len(perenn_links)}") - if not perenn_links.empty: - print(perenn_links[["bus0", "bus1", "carrier", "p_nom"]].to_string(index=True)) - - # Check stores - perenn_stores = n.stores[n.stores.carrier.str.contains("perennial", case=False, na=False)] - print(f" Stores (carrier=perennial store): {len(perenn_stores)}") - if not perenn_stores.empty: - print(perenn_stores[["bus", "carrier", "e_nom_max"]].head(10).to_string(index=True)) - if "e_nom_max" in perenn_stores.columns: - finite_max = perenn_stores["e_nom_max"][perenn_stores["e_nom_max"] < 1e18] - total_cap = finite_max.sum() - print(f"\n Total store capacity (e_nom_max): {total_cap:,.0f} tCO2 ({total_cap/1e6:.3f} MtCO2)") - if total_cap == 0: - print(f" {WARN} ALL stores have e_nom_max=0!") - print(f" {WARN} This usually means biomass.classes in config.default.yaml") - print(f" {WARN} is missing biofuels_1G_* entries — check and re-run prepare_sector_network.") - - if not perenn_carriers: - print(f"\n{WARN} No perennial carriers found — add_perennials may NOT have run!") - else: - print(f"\n{OK} Perennial components found in pre-network.") - - except Exception as e: - print(f"\n{WARN} Could not load network: {e}") -else: - print(f"\n{FAIL} Pre-network missing — prepare_sector_network has not run yet.") - -# ── 4. Optimal solution ────────────────────────────────────────────────────── -section("4. OPTIMAL SOLUTION: solved network (solve_sector_network)") - -opt_path = RESULTS / "networks" / f"{WC}.nc" -opt_ok = check_file(opt_path, f"Optimal network {WC}.nc") - -if opt_ok: - try: - import pypsa - n_opt = pypsa.Network(str(opt_path)) - - perenn_links = n_opt.links[n_opt.links.carrier.str.contains("perennial", case=False, na=False)] - perenn_stores = n_opt.stores[n_opt.stores.carrier.str.contains("perennial", case=False, na=False)] - - print(f"\n Links (carrier=perennial): {len(perenn_links)}") - if not perenn_links.empty: - cols = ["carrier", "p_nom_opt"] if "p_nom_opt" in perenn_links.columns else ["carrier", "p_nom"] - print(perenn_links[cols].to_string(index=True)) - if "p_nom_opt" in perenn_links.columns: - active = perenn_links[perenn_links["p_nom_opt"] > 0] - print(f"\n Links with p_nom_opt > 0: {len(active)}") - if active.empty: - print(f"{WARN} Perennial links exist but have zero optimal capacity.") - else: - print(f"{OK} Perennial links are deployed in the optimal solution.") - - print(f"\n Stores (carrier=perennial store): {len(perenn_stores)}") - if not perenn_stores.empty: - cols = ["carrier", "e_nom_opt"] if "e_nom_opt" in perenn_stores.columns else ["carrier", "e_nom"] - print(perenn_stores[cols].to_string(index=True)) - if "e_nom_opt" in perenn_stores.columns: - active = perenn_stores[perenn_stores["e_nom_opt"] > 0] - total = perenn_stores["e_nom_opt"].sum() - print(f"\n Stores with e_nom_opt > 0: {len(active)}") - print(f" Total e_nom_opt: {total:,.0f} tCO2 ({total / 1e6:.3f} MtCO2)") - if active.empty: - print(f"{WARN} All perennial stores have e_nom_opt = 0 (not deployed).") - else: - print(f"\n Per-node e_nom_opt [tCO2] stats:") - print(f" min: {perenn_stores['e_nom_opt'].min():,.0f}") - print(f" mean: {perenn_stores['e_nom_opt'].mean():,.0f}") - print(f" max: {perenn_stores['e_nom_opt'].max():,.0f}") - - if perenn_links.empty and perenn_stores.empty: - print(f"\n{WARN} No perennial components in optimal network!") - - except Exception as e: - print(f"\n{WARN} Could not load optimal network: {e}") -else: - print(f"\n{FAIL} Optimal network missing — solve_sector_network has not run yet.") - -# ── Summary ────────────────────────────────────────────────────────────────── -section("SUMMARY") -print(f" Run: {RDIR}") -print(f" Wildcard: {WC}") -print(f" Shared resources:{SHARED_RES.resolve()}") -print(f" Results dir: {RESULTS.resolve()}") From 9de8c9eb6a021dfdeabff1bbe1b1276f9d0c83b2 Mon Sep 17 00:00:00 2001 From: BertoGBG Date: Thu, 9 Jul 2026 12:02:21 +0200 Subject: [PATCH 04/10] 1G biofuels are substituted is are located in the not included biomass class --- config/config.default.yaml | 4 --- rules/build_sector.smk | 3 +- scripts/build_perennials_crop_yields_nuts2.py | 28 +++++++++++++------ scripts/prepare_sector_network.py | 3 +- 4 files changed, 21 insertions(+), 17 deletions(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index 05bcb35da7..6a28a57ec1 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -422,10 +422,6 @@ renewable: # alternative carbon dioxide removal technologies perennials: potential_co2: 2 # tCO2e/(ha y) sequestered - biofuel_conversion: # Biofuel conversion factors from JRC Technical Report doi:10.2760/69179 - MINBIOCRP11: 0.295 # t_ethanol / t_wheat grain (13.5% moisture), Table 93 - MINBIOCRP21: 0.07777 # t_ethanol / t_sugar beet (16% sugar content), Table 133 - MINBIORPS1: 0.4176 # t_crude_oil / t_rapeseed (9% moisture) x crude-to-FAME (1/1.0063), Tables 155+159 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#conventional conventional: diff --git a/rules/build_sector.smk b/rules/build_sector.smk index c876b3f3f6..352fd8257a 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -921,12 +921,11 @@ rule build_biomass_potentials: rule build_perennials_yields_nuts_file: - params: - biofuel_conversion=config_provider("perennials", "biofuel_conversion"), input: nuts2021=rules.retrieve_eu_nuts_2021.output.shapes_level_2, crops_nuts2=rules.retrieve_co2_removal_data.output.eurostat_crops_nuts2, crops_nuts0=rules.retrieve_co2_removal_data.output.eurostat_crops_nuts0, + costs=resources(f"costs_{config['costs']['year']}_processed.csv"), output: yields_all=resources("perennials_yields_1G_biofuels.csv"), log: diff --git a/scripts/build_perennials_crop_yields_nuts2.py b/scripts/build_perennials_crop_yields_nuts2.py index 43d47f82ee..e505a48b3a 100644 --- a/scripts/build_perennials_crop_yields_nuts2.py +++ b/scripts/build_perennials_crop_yields_nuts2.py @@ -10,21 +10,23 @@ Outputs a single CSV with columns for each crop class (cereals, sugar beet, rapeseed, perennials) indexed by NUTS2 region. -Biofuel conversion efficiencies (t_biofuel / t_feedstock) are read from -``config["perennials"]["biofuel_conversion"]`` and sourced from: +Biofuel conversion efficiencies (t_biofuel / t_feedstock) are read from the +``efficiency`` parameter of the ``ethanol from wheat``, ``ethanol from sugar +beet``, and ``biodiesel from rapeseed`` technologies in the technology-data +cost assumptions, sourced from: Banja et al. (2013), "Biofuels in the European Union - A general overview", JRC Technical Report, doi:10.2760/69179, Tables 93, 133, 155, 159. """ import logging -import os from pathlib import Path import geopandas as gpd import numpy as np import pandas as pd -import requests + +from scripts._helpers import load_costs logger = logging.getLogger(__name__) @@ -216,7 +218,7 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b * unsustainable_biofuels_yields.index.get_level_values("mapping").map(biofuel_yields) ) - # yields of perennials per hectar in ton/ha + # yields of perennials per hectare in ton/ha # standard humidity for perennials = 0.65 (tH2O/t_fresh) -> note production is for fresh until 2025 std_moist_perennials = 0.65 @@ -255,6 +257,11 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b CROPS_CSV_NUTS2.parent.mkdir(parents=True, exist_ok=True) OUT_CSV_YIELDS_ALL.parent.mkdir(parents=True, exist_ok=True) + # G0000: total green plants (kept – often the only G-code with NUTS2 coverage, e.g. DK) + # G1000: temporary grasses and grazings + # G2000: aggregate legumes (G2100 + G2900) – included for NUTS0 fallback coverage; + # does NOT inflate MAX because it is always ≤ max(G2100, G2900) + # G2100: lucerne/alfalfa; G2900: clover and other leguminous plants perennial_codes = ["G0000", "G1000", "G2000", "G2100", "G2900"] crops_mapping = dict( @@ -263,15 +270,16 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b MINBIORPS1=["I1110", "I1120", "I1130", "I1110-1130", "I0000"], PERENNIALS=perennial_codes, ) + costs = load_costs(snakemake.input.costs) - conv = snakemake.params.biofuel_conversion + #conv = snakemake.params.biofuel_conversion # LHV values are fixed physical constants, not parameters: JRC Technical Report doi:10.2760/69179 LHV_fuels = {"ethanol": 7.447, "biodiesel": 10.194} # MWh/t (26.81 MJ/kg, 36.7 MJ/kg) biofuel_yields = { - "MINBIOCRP11": conv["MINBIOCRP11"] * LHV_fuels["ethanol"], - "MINBIOCRP21": conv["MINBIOCRP21"] * LHV_fuels["ethanol"], - "MINBIORPS1": conv["MINBIORPS1"] * LHV_fuels["biodiesel"], + "MINBIOCRP11": costs.at['ethanol from wheat', 'efficiency'] * LHV_fuels["ethanol"], + "MINBIOCRP21": costs.at['ethanol from sugar beet', 'efficiency'] * LHV_fuels["ethanol"], + "MINBIORPS1": costs.at['biodiesel from rapeseed', 'efficiency'] * LHV_fuels["biodiesel"], } other_crops_codes = [ @@ -339,6 +347,8 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b "MINBIOCRP11": yield_MINBIOCRP11_full["energy_yields_(MWh/ha)"], "MINBIOCRP21": yield_MINBIOCRP21_full["energy_yields_(MWh/ha)"], "MINBIORPS1": yield_MINBIORPS1_full["energy_yields_(MWh/ha)"], + # Max yield across G-codes: models the best-available perennial crop + # choice in each region when substituting 1G biofuel crops. "PERENNIALS_MAX": yields_perennials_max_full["YL_(t/ha)"], }, axis=1, diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index f4b7b262c1..8b94eb81f6 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -1345,8 +1345,7 @@ def add_perennials(n, costs): nodes, suffix=" CO2s_perennials", bus=nodes + " perennials co2 store", - e_nom_extendable=True, - e_nom_max=perennials_potentials_spatial.values, + e_nom=perennials_potentials_spatial.values, carrier="co2 perennials", e_cyclic=False, ) From e17bcd578a6ca1d6400260cd6162d8437c0abce1 Mon Sep 17 00:00:00 2001 From: BertoGBG Date: Wed, 19 Aug 2026 18:37:46 +0200 Subject: [PATCH 05/10] pre PR changes --- config/config.default.yaml | 2 +- rules/build_sector.smk | 12 +- scripts/build_biomass_potentials.py | 10 +- ...tentials.py => build_perennials_yields.py} | 129 ++++++++++++------ ...ild_perennials_yields_eurostat_average.py} | 18 ++- scripts/prepare_sector_network.py | 51 +++---- 6 files changed, 137 insertions(+), 85 deletions(-) rename scripts/{build_perennials_potentials.py => build_perennials_yields.py} (56%) rename scripts/{build_perennials_crop_yields_nuts2.py => build_perennials_yields_eurostat_average.py} (94%) diff --git a/config/config.default.yaml b/config/config.default.yaml index 6a28a57ec1..b43c33cf0b 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -421,7 +421,7 @@ renewable: # alternative carbon dioxide removal technologies perennials: - potential_co2: 2 # tCO2e/(ha y) sequestered + sequestration_co2: 2 # tCO2e/(ha y) sequestered # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#conventional conventional: diff --git a/rules/build_sector.smk b/rules/build_sector.smk index 352fd8257a..b775620743 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -920,7 +920,7 @@ rule build_biomass_potentials: scripts("build_biomass_potentials.py") -rule build_perennials_yields_nuts_file: +rule build_perennials_yields_eurostat_average: input: nuts2021=rules.retrieve_eu_nuts_2021.output.shapes_level_2, crops_nuts2=rules.retrieve_co2_removal_data.output.eurostat_crops_nuts2, @@ -929,12 +929,12 @@ rule build_perennials_yields_nuts_file: output: yields_all=resources("perennials_yields_1G_biofuels.csv"), log: - logs("build_perennials_yields_nuts_file.log"), + logs("build_perennials_yields_eurostat_average.log"), script: - scripts("build_perennials_crop_yields_nuts2.py") + scripts("build_perennials_yields_eurostat_average.py") -rule build_perennial_potentials: +rule build_perennials_yields: params: biomass=config_provider("biomass"), input: @@ -945,11 +945,11 @@ rule build_perennial_potentials: output: csv_file=resources("perennials_yields_1G_biofuels_s_{clusters}.csv"), log: - logs("build_perennial_potentials_s_{clusters}.log"), + logs("build_perennials_yields_s_{clusters}.log"), resources: mem_mb=8000, script: - scripts("build_perennials_potentials.py") + scripts("build_perennials_yields.py") rule build_biomass_transport_costs: diff --git a/scripts/build_biomass_potentials.py b/scripts/build_biomass_potentials.py index 70ab27d80c..78aa1e91e6 100755 --- a/scripts/build_biomass_potentials.py +++ b/scripts/build_biomass_potentials.py @@ -184,17 +184,15 @@ def disaggregate_nuts0(bio): return bio -def build_nuts2_shapes(): +def build_nuts2_shapes(nuts2_fn, country_shapes_fn): """ - 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("NUTS_ID").geometry - ) + nuts2 = gpd.GeoDataFrame(gpd.read_file(nuts2_fn).set_index("NUTS_ID").geometry) - countries = gpd.read_file(snakemake.input.country_shapes).set_index("name") + countries = gpd.read_file(country_shapes_fn).set_index("name") missing_iso2 = countries.index.intersection(["AL", "RS", "XK", "BA"]) missing = countries.loc[missing_iso2] @@ -378,7 +376,7 @@ def add_unsustainable_potentials(df, input_eurostat): enspreso = disaggregate_nuts0(enspreso) - nuts2 = build_nuts2_shapes() + nuts2 = build_nuts2_shapes(snakemake.input.nuts2, snakemake.input.country_shapes) df_nuts2 = gpd.GeoDataFrame(nuts2.geometry).join(enspreso) diff --git a/scripts/build_perennials_potentials.py b/scripts/build_perennials_yields.py similarity index 56% rename from scripts/build_perennials_potentials.py rename to scripts/build_perennials_yields.py index 747e50ac57..fec978c9cd 100644 --- a/scripts/build_perennials_potentials.py +++ b/scripts/build_perennials_yields.py @@ -2,23 +2,37 @@ # # SPDX-License-Identifier: MIT """ -Reproject NUTS2-level 1st-generation (1G) biofuel and perennial crop yields -(from ``build_perennials_crop_yields_nuts2.py``) onto the clustered network +Reproject NUTS2-level 1st-generation (1G) biofuel and perennial crop YIELDS +(from ``build_perennials_yields_eurostat_average.py``) onto the clustered network regions, for use by ``add_perennials()`` in ``prepare_sector_network.py``. -NUTS2 yields (MWh/ha/y for 1G biofuels, t/ha/y for perennials) are mapped to -clustered regions via an area-weighted overlay (NUTS2 geometries intersected -with cluster region geometries, weighted by intersection area); NUTS2 -regions not covered by Eurostat data (non-EU countries, small islands, -city-states) are first filled from the nearest valid NUTS2 centroid. -Resulting columns are then grouped into biomass classes -(``resolve_biomass_classes()``) to match the class structure used by -``build_biomass_potentials.py``. - -``add_perennials()`` later combines this script's clustered-region 1G yields -with the ENSPRESO 1G biomass potential (MWh/y) to back out the land area -available for conversion to perennial grasses, and converts that area into a -CO2 sequestration potential via ``perennials.potential_co2`` (tCO2/ha/y). +Output units, NOT an area or a potential +----------------------------------------- +The output CSV holds per-hectare yield RATES, not land area or a CO2 +potential: MWh/ha/y for the 1G biofuel crop columns (cereals, sugar beet, +rapeseed) and t/ha/y for the ``perennials`` column. The land area available +for conversion to perennial grasses, and the resulting CO2 sequestration +potential, are only derived later - see "Downstream usage" below. + +Method +------ +NUTS2 yields are mapped to clustered regions via an area-weighted average +over the NUTS2/region overlay (NUTS2 geometries intersected with cluster +region geometries, weighted by intersection area) - see +``convert_nuts2_to_regions_yields()``. NUTS2 regions not covered by +Eurostat data (non-EU countries, small islands, city-states) are first +filled from the nearest valid NUTS2 centroid - see +``impute_missing_values()``. Resulting columns are then grouped into +biomass classes (``resolve_biomass_classes()``) to match the class +structure used by ``build_biomass_potentials.py``. + +Downstream usage +----------------- +``add_perennials()`` divides the ENSPRESO 1G biomass potential (MWh/y, an +extensive quantity from ``build_biomass_potentials.py``) by this script's +1G yields (MWh/ha/y) to back out the land area (ha) available for +conversion to perennial grasses, then converts that area into a CO2 +sequestration potential via ``perennials.sequestration_co2`` (tCO2e/ha/y). Outputs a single CSV with one column per crop class (cereals, sugar beet, rapeseed, perennials) indexed by clustered region name. @@ -28,36 +42,42 @@ import geopandas as gpd import pandas as pd -from _helpers import configure_logging, resolve_biomass_classes, set_scenario_config + +from scripts._helpers import ( + configure_logging, + resolve_biomass_classes, + set_scenario_config, +) +from scripts.build_biomass_potentials import build_nuts2_shapes logger = logging.getLogger(__name__) -def build_nuts2_shapes(): +def impute_missing_values(df_nuts2, missing_shapes, yield_cols): """ - - load NUTS2 geometries - - add RS, AL, BA country shapes (not covered in NUTS 2013) - - consistently name ME, MK + Fill missing yield values for NUTS2 regions not covered by Eurostat data + (non-EU countries, small islands, city-states) by copying the values from + the nearest NUTS2 region that does have valid (non-NaN) yields. + + Nearest is determined by centroid distance in an equal-area CRS (EPSG:3035). + + Parameters + ---------- + df_nuts2 : gpd.GeoDataFrame + NUTS2 geometries joined with yield columns; rows for regions without + Eurostat coverage are entirely NaN in ``yield_cols``. + missing_shapes : gpd.GeoDataFrame + Geometries of the NUTS2 (or country-level substitute) regions to + impute, indexed the same way as ``df_nuts2``. + yield_cols : list of str + Columns in ``df_nuts2`` to impute. + + Returns + ------- + gpd.GeoDataFrame + One row per entry in ``missing_shapes``, with ``yield_cols`` filled + from the nearest valid NUTS2 neighbour. """ - nuts2 = gpd.GeoDataFrame( - gpd.read_file(snakemake.input.nuts2).set_index("NUTS_ID").geometry - ) - - countries = gpd.read_file(snakemake.input.country_shapes).set_index("name") - missing_iso2 = countries.index.intersection(["AL", "RS", "XK", "BA"]) - missing = countries.loc[missing_iso2] - - nuts2.rename(index={"ME00": "ME", "MK00": "MK"}, inplace=True) - - return pd.concat([nuts2, missing]) - - -def area(gdf): - return gdf.to_crs(epsg=3035).area.div(1e6) - - -def impute_missing_values(df_nuts2, missing_shapes, yield_cols): - # Keep only rows that have valid yields (drop NaN rows!) df_valid = df_nuts2.dropna(subset=yield_cols).copy() @@ -91,12 +111,35 @@ def impute_missing_values(df_nuts2, missing_shapes, yield_cols): def convert_nuts2_to_regions_yields(df_nuts2, regions, yield_cols=None): """ - Convert NUTS2-level yields (intensive) to PyPSA regions using: + Convert NUTS2-level yields (intensive, e.g. MWh/ha/y or t/ha/y) to + PyPSA-Eur clustered regions via an area-weighted average over the + NUTS2/region overlay: y_n = Σ_i (y_i * A_i∩n) / Σ_i A_i∩n + Unlike an extensive quantity (a total, e.g. MWh), a yield is a rate and + must be area-weighted-averaged rather than redistributed by area share + - see ``convert_nuts2_to_regions`` in ``build_biomass_potentials.py`` for + the extensive-quantity equivalent. + Only NUTS2 rows with non-NaN yields are used. Regions with no overlapping valid NUTS2 get NaN. + + Parameters + ---------- + df_nuts2 : gpd.GeoDataFrame + NUTS2 geometries joined with yield columns. + regions : gpd.GeoDataFrame + PyPSA-Eur clustered onshore regions, with a ``name`` column. + yield_cols : list of str, optional + Columns in ``df_nuts2`` to convert. Defaults to all columns except + ``geometry`` and ``NUTS_ID``. + + Returns + ------- + pd.DataFrame + Area-weighted-average yields indexed by region name, one row per + entry in ``regions``. """ nuts = df_nuts2.copy() @@ -147,10 +190,10 @@ def convert_nuts2_to_regions_yields(df_nuts2, regions, yield_cols=None): if __name__ == "__main__": if "snakemake" not in globals(): - from _helpers import mock_snakemake + from scripts._helpers import mock_snakemake snakemake = mock_snakemake( - "build_perennials_potentials", + "build_perennials_yields", clusters="39", planning_horizons=2050, ) @@ -159,7 +202,7 @@ def convert_nuts2_to_regions_yields(df_nuts2, regions, yield_cols=None): set_scenario_config(snakemake) regions = gpd.read_file(snakemake.input.regions_onshore) - nuts2 = build_nuts2_shapes() + nuts2 = build_nuts2_shapes(snakemake.input.nuts2, snakemake.input.country_shapes) yields = pd.read_csv(snakemake.input.perennials_yields_1G_biofuels, index_col=0) diff --git a/scripts/build_perennials_crop_yields_nuts2.py b/scripts/build_perennials_yields_eurostat_average.py similarity index 94% rename from scripts/build_perennials_crop_yields_nuts2.py rename to scripts/build_perennials_yields_eurostat_average.py index e505a48b3a..31e06996d9 100644 --- a/scripts/build_perennials_crop_yields_nuts2.py +++ b/scripts/build_perennials_yields_eurostat_average.py @@ -3,9 +3,19 @@ # SPDX-License-Identifier: MIT """ Retrieve crop harvest data from the Eurostat API (dataset ``apro_cpshr``) at -NUTS2 and NUTS0 resolution, compute area-weighted yields for 1st-generation -biofuel feedstocks and perennial grasses, and harmonize the results to the -NUTS2021 region definitions used by PyPSA-Eur. +NUTS2 and NUTS0 resolution, and compute production-weighted average yields +(t/ha, then converted to MWh/ha for 1G biofuel crops) per NUTS2 region for +1st-generation biofuel feedstocks and perennial grasses. + +"Eurostat average" in the name refers to this production-weighted averaging +across crop codes and years - not the geometric area-weighted overlay used +downstream in ``build_perennials_yields.py`` to reproject these NUTS2-level +yields onto clustered network regions; do not confuse the two. + +Missing NUTS2 coverage is filled via a three-tier fallback (NUTS2 data -> +NUTS0 country-level data -> spatial neighbor mean or nearest valid NUTS2 for +islands - see ``harmonize_to_nuts2021()``), and results are harmonized to +the NUTS2021 region definitions used by PyPSA-Eur. Outputs a single CSV with columns for each crop class (cereals, sugar beet, rapeseed, perennials) indexed by NUTS2 region. @@ -242,7 +252,7 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b if "snakemake" not in globals(): from scripts._helpers import mock_snakemake - snakemake = mock_snakemake("build_perennials_yields_nuts_file") + snakemake = mock_snakemake("build_perennials_yields_eurostat_average") from scripts._helpers import configure_logging, set_scenario_config diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 8b94eb81f6..241d5c5650 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -1258,7 +1258,7 @@ def add_perennials(n, costs): potential already allocated to 1G biofuels (``biomass_potentials``, MWh/y) divided by the 1G crop yield (MWh/ha/y) at that node, giving a displaced area in ha; multiplying by a fixed CO2 sequestration rate per - hectare (``perennials.potential_co2``) gives the store's CO2 potential. + hectare (``perennials.sequestration_co2``) gives the store's CO2 potential. A single "co2 perennials" Link models the harvesting process: CO2 drawn from the atmosphere (bus0) is converted into biogas (bus3) and stored CO2 (bus1), with capacity restricted to the April-October harvesting season @@ -1270,7 +1270,7 @@ def add_perennials(n, costs): The PyPSA network container object costs : pd.DataFrame Costs and parameters for different technologies. Must contain a - 'perennials gbr' entry with 'electricity-input', 'biogas-output', + 'perennials refining' entry with 'electricity-input', 'biogas-output', 'capital_cost', 'VOM', and 'lifetime' parameters Returns @@ -1284,9 +1284,9 @@ def add_perennials(n, costs): Reads ``snakemake.input.biomass_potentials`` and ``snakemake.input.perennials_yields_1G_biofuels`` (NUTS2-derived crop yields aggregated to clustered network regions, see - ``build_perennials_crop_yields_nuts2.py`` and - ``build_perennials_potentials.py``), and - ``snakemake.config["perennials"]["potential_co2"]``. + ``build_perennials_yields_eurostat_average.py`` and + ``build_perennials_yields.py``), and + ``snakemake.config["perennials"]["sequestration_co2"]``. """ logger.info("Adding perennials.") @@ -1295,57 +1295,58 @@ def add_perennials(n, costs): biomass_potentials = pd.read_csv(snakemake.input.biomass_potentials, index_col=0) perennials_yields_1G_biofuels = pd.read_csv(snakemake.input.perennials_yields_1G_biofuels).set_index("name") - # calculate CO2 sequestration per tDM perennials - perennial_CO2_seq = perennials_yields_1G_biofuels["perennials"] / snakemake.config["perennials"]["potential_co2"] # (tDM/tCO2 seq) - - # calculate perennials potential based on conversion on first generation biofuels - perennials_area_spatial = (biomass_potentials.filter(regex='biofuels_1G') / perennials_yields_1G_biofuels.filter(regex='biofuels_1G')).sum(axis=1) + # calculate perennials potential based on the conversion on first generation biofuels for equal area + perennials_area = (biomass_potentials.filter(regex='biofuels_1G') / perennials_yields_1G_biofuels.filter(regex='biofuels_1G')).sum(axis=1) # (MWh/y) / (MWh / ha / y) = (ha) returns the area used by sum of the 3 biofuels_1G classes which can be assigned for perennials - perennials_potentials_spatial = perennials_area_spatial * snakemake.config["perennials"]["potential_co2"] # (tCO2seq) = (ha) * (tCO2 seq/ha) + perennials_potentials = perennials_area * snakemake.config["perennials"]["sequestration_co2"] # (tCO2seq) = (ha) * (tCO2 seq/ha) nodes = pop_layout.index n.add("Carrier", "co2 perennials") n.add( "Bus", - nodes + " perennials co2 store", + nodes, + suffix=" co2 perennials", location=nodes, carrier="co2 perennials", unit="t_co2", ) + # calculate CO2 sequestration per tDM perennials + perennial_CO2_seq = perennials_yields_1G_biofuels["perennials"] / snakemake.config["perennials"]["sequestration_co2"] # (tDM/tCO2 seq) + # calculate biogas production based on harvesting time (in month) - df_gbr = pd.DataFrame(index=n.snapshots, columns=["harvest"]) - df_gbr["harvest"] = df_gbr.index.month.isin([4, 5, 6, 7, 8, 9, 10]).astype(int) + df_harvest = pd.DataFrame(index=n.snapshots, columns=["harvest"]) + df_harvest["harvest"] = df_harvest.index.month.isin([4, 5, 6, 7, 8, 9, 10]).astype(int) p_max_pu = pd.DataFrame(index=n.snapshots, columns=nodes) for node in nodes: - p_max_pu[node] = df_gbr["harvest"] + p_max_pu[node] = df_harvest["harvest"] n.add( "Link", nodes, - suffix=" perennials GBR", + suffix=" perennials refining", bus0="co2 atmosphere", - bus1=nodes + " perennials co2 store", + bus1=nodes + " co2 perennials", bus2=nodes.values, bus3=spatial.gas.biogas, efficiency=1, - efficiency2=-costs.at["perennials gbr", "electricity-input"] * perennial_CO2_seq, - efficiency3=costs.at["perennials gbr", "biogas-output"] * perennial_CO2_seq, + efficiency2=-costs.at["perennials refining", "electricity-input"] * perennial_CO2_seq, + efficiency3=costs.at["perennials refining", "biogas-output"] * perennial_CO2_seq, carrier="co2 perennials", p_nom_extendable=True, p_max_pu=p_max_pu, - capital_cost=costs.at["perennials gbr", "capital_cost"] * perennial_CO2_seq, - marginal_cost=costs.at["perennials gbr", "VOM"] * perennial_CO2_seq, - lifetime=costs.at["perennials gbr", "lifetime"], + capital_cost=costs.at["perennials refining", "capital_cost"] * perennial_CO2_seq, + marginal_cost=costs.at["perennials refining", "VOM"] * perennial_CO2_seq, + lifetime=costs.at["perennials refining", "lifetime"], ) n.add( "Store", nodes, - suffix=" CO2s_perennials", - bus=nodes + " perennials co2 store", - e_nom=perennials_potentials_spatial.values, + suffix=" CO2s perennials", + bus=nodes + " co2 perennials", + e_nom=perennials_potentials.values, carrier="co2 perennials", e_cyclic=False, ) From 7e24967f1521921daf82707c8e89aa41963ce1e4 Mon Sep 17 00:00:00 2001 From: BertoGBG Date: Thu, 20 Aug 2026 14:33:17 +0200 Subject: [PATCH 06/10] Add config schema validation for perennials; add release note Mirrors the pattern set by rock_weathering_PR (scripts/lib/validation/config/rock_weathering.py): - New scripts/lib/validation/config/perennials.py with PerennialsConfig (sequestration_co2), registered on the top-level ConfigSchema in _schema.py. - sector.py: added sector.perennials: bool (the CDR technology toggle). - data.py: added data.co2_removal_data as a _DataSourceConfig (feeds the retrieve_co2_removal_data rule's Eurostat crop yield inputs). - Regenerated config/config.default.yaml and config/schema.default.json from the new models; both config-schema sync tests now pass (they were failing before, same root cause as biochar_PR's still-open gap: perennials config wasn't covered by any pydantic model). Also added a release_notes.md entry matching rock_weathering_PR's "Upcoming Release" bullet style, referencing tracking issue #2143. This commit was preceded by a merge of upstream/master (57 commits behind at the time) to pull in the validation framework this depends on. That merge had exactly one conflict, in scripts/_helpers.py: both branches had independently appended unrelated functions to the end of the file (resolve_biomass_classes() here vs. _simplify_polys()/ load_data_versions() upstream) - resolved by keeping both. Full test suite (29 tests) and ruff both pass. Verified via Snakemake dry-run with sector.perennials temporarily set to true that the full 77-job DAG resolves end-to-end. Co-Authored-By: Claude Sonnet 5 Claude-Session: https://claude.ai/code/session_01Mh2EaBSg63Gsi4epA1Rmky --- config/config.default.yaml | 4 +- config/schema.default.json | 70 +++++++++++++++++++++ doc/release_notes.md | 1 + scripts/lib/validation/config/_schema.py | 5 ++ scripts/lib/validation/config/data.py | 4 ++ scripts/lib/validation/config/perennials.py | 22 +++++++ scripts/lib/validation/config/sector.py | 4 ++ 7 files changed, 108 insertions(+), 2 deletions(-) create mode 100644 scripts/lib/validation/config/perennials.py diff --git a/config/config.default.yaml b/config/config.default.yaml index 58eab0487d..0b66f0d7b5 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -420,9 +420,9 @@ renewable: eia_correct_by_capacity: false eia_approximate_missing: false -# alternative carbon dioxide removal technologies +# docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#perennials perennials: - sequestration_co2: 2 # tCO2e/(ha y) sequestered + sequestration_co2: 2 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#conventional conventional: diff --git a/config/schema.default.json b/config/schema.default.json index a442c12f17..72ac2f3e37 100644 --- a/config/schema.default.json +++ b/config/schema.default.json @@ -1044,6 +1044,26 @@ } } }, + "co2_removal_data": { + "description": "Configuration for a single data source.", + "properties": { + "source": { + "default": "archive", + "description": "Source of the data. 'archive' retrieves pre-built data, 'primary' retrieves from primary source.", + "enum": [ + "archive", + "primary", + "build" + ], + "type": "string" + }, + "version": { + "default": "latest", + "description": "Version of the data to use. Uses the specific 'version' for the selected 'source' or the dataset tagged 'latest' for this source.", + "type": "string" + } + } + }, "co2stop": { "description": "Configuration for a single data source.", "properties": { @@ -2925,6 +2945,16 @@ } } }, + "PerennialsConfig": { + "description": "Configuration for `perennials` settings.", + "properties": { + "sequestration_co2": { + "default": 2, + "description": "Tonnes of CO2 equivalent sequestered per hectare per year when 1st-generation biofuel cropland is converted to perennial grasses.", + "type": "number" + } + } + }, "PypsaEurConfig": { "description": "Configuration for `pypsa_eur` settings.", "properties": { @@ -4647,6 +4677,11 @@ "description": "Add option for Direct Air Capture (DAC).", "type": "boolean" }, + "perennials": { + "default": false, + "description": "Add option for perennialisation (converting 1st-generation biofuel cropland to perennial grasses) as a carbon dioxide removal (CDR) technology.", + "type": "boolean" + }, "co2_vent": { "default": false, "description": "Add option for vent out CO2 from storages to the atmosphere.", @@ -10091,6 +10126,16 @@ } } }, + "perennials": { + "description": "Configuration for `perennials` settings.", + "properties": { + "sequestration_co2": { + "default": 2, + "description": "Tonnes of CO2 equivalent sequestered per hectare per year when 1st-generation biofuel cropland is converted to perennial grasses.", + "type": "number" + } + } + }, "conventional": { "additionalProperties": true, "description": "Configuration for `conventional` settings.", @@ -11139,6 +11184,11 @@ "description": "Add option for Direct Air Capture (DAC).", "type": "boolean" }, + "perennials": { + "default": false, + "description": "Add option for perennialisation (converting 1st-generation biofuel cropland to perennial grasses) as a carbon dioxide removal (CDR) technology.", + "type": "boolean" + }, "co2_vent": { "default": false, "description": "Add option for vent out CO2 from storages to the atmosphere.", @@ -12901,6 +12951,26 @@ } } }, + "co2_removal_data": { + "description": "Configuration for a single data source.", + "properties": { + "source": { + "default": "archive", + "description": "Source of the data. 'archive' retrieves pre-built data, 'primary' retrieves from primary source.", + "enum": [ + "archive", + "primary", + "build" + ], + "type": "string" + }, + "version": { + "default": "latest", + "description": "Version of the data to use. Uses the specific 'version' for the selected 'source' or the dataset tagged 'latest' for this source.", + "type": "string" + } + } + }, "co2stop": { "description": "Configuration for a single data source.", "properties": { diff --git a/doc/release_notes.md b/doc/release_notes.md index 8fcdee22ae..c699ab3e2b 100644 --- a/doc/release_notes.md +++ b/doc/release_notes.md @@ -5,6 +5,7 @@ +* feat: Add perennialisation as a carbon dioxide removal (CDR) technology, converting 1st-generation biofuel cropland to perennial grasses, with node-level potential derived from NUTS2-resolved Eurostat crop yields ([#2143](https://github.com/PyPSA/pypsa-eur/issues/2143)). ## PyPSA-Eur v2026.08.0 (19th August 2026) diff --git a/scripts/lib/validation/config/_schema.py b/scripts/lib/validation/config/_schema.py index ad9e941c00..9f982ba4a7 100644 --- a/scripts/lib/validation/config/_schema.py +++ b/scripts/lib/validation/config/_schema.py @@ -27,6 +27,7 @@ from scripts.lib.validation.config.links import LinksConfig from scripts.lib.validation.config.load import LoadConfig from scripts.lib.validation.config.overpass_api import OverpassApiConfig +from scripts.lib.validation.config.perennials import PerennialsConfig from scripts.lib.validation.config.pypsa_eur import PypsaEurConfig from scripts.lib.validation.config.renewable import RenewableConfig from scripts.lib.validation.config.run import RunConfig @@ -140,6 +141,10 @@ class ConfigSchema(BaseModel): default_factory=RenewableConfig, description="Renewable energy technologies configuration.", ) + perennials: PerennialsConfig = Field( + default_factory=PerennialsConfig, + description="Perennialisation (carbon dioxide removal) configuration.", + ) conventional: ConventionalConfig = Field( default_factory=ConventionalConfig, description="Conventional power plants configuration.", diff --git a/scripts/lib/validation/config/data.py b/scripts/lib/validation/config/data.py index 18477518fd..bdc0dfddf4 100644 --- a/scripts/lib/validation/config/data.py +++ b/scripts/lib/validation/config/data.py @@ -201,6 +201,10 @@ def check_version_files_are_correct_suffix( default_factory=lambda: _DataSourceConfig(source="primary"), description="Instrat CO2 prices data source configuration.", ) + co2_removal_data: _DataSourceConfig = Field( + default_factory=lambda: _DataSourceConfig(source="primary"), + description="Carbon dioxide removal (afforestation, perennialisation) input data source configuration.", + ) co2stop: _DataSourceConfig = Field( default_factory=_DataSourceConfig, description="CO2Stop data source configuration.", diff --git a/scripts/lib/validation/config/perennials.py b/scripts/lib/validation/config/perennials.py new file mode 100644 index 0000000000..a7abc3f631 --- /dev/null +++ b/scripts/lib/validation/config/perennials.py @@ -0,0 +1,22 @@ +# SPDX-FileCopyrightText: Contributors to PyPSA-Eur +# +# SPDX-License-Identifier: MIT + +""" +Perennialisation configuration. + +See docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#perennials +""" + +from pydantic import Field + +from scripts.lib.validation.config._base import ConfigModel + + +class PerennialsConfig(ConfigModel): + """Configuration for `perennials` settings.""" + + sequestration_co2: float = Field( + 2, + description="Tonnes of CO2 equivalent sequestered per hectare per year when 1st-generation biofuel cropland is converted to perennial grasses.", + ) diff --git a/scripts/lib/validation/config/sector.py b/scripts/lib/validation/config/sector.py index fa561f2adf..df7370a8e5 100644 --- a/scripts/lib/validation/config/sector.py +++ b/scripts/lib/validation/config/sector.py @@ -669,6 +669,10 @@ class SectorConfig(BaseModel): False, description="Add option for coal CHPs with carbon capture." ) dac: bool = Field(True, description="Add option for Direct Air Capture (DAC).") + perennials: bool = Field( + False, + description="Add option for perennialisation (converting 1st-generation biofuel cropland to perennial grasses) as a carbon dioxide removal (CDR) technology.", + ) co2_vent: bool = Field( False, description="Add option for vent out CO2 from storages to the atmosphere.", From 9d746fcddd7c059b138d6a6f03fa906a877280fc Mon Sep 17 00:00:00 2001 From: BertoGBG Date: Thu, 20 Aug 2026 15:06:56 +0200 Subject: [PATCH 07/10] new rules documented in doc.md --- doc/sector.md | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/doc/sector.md b/doc/sector.md index 549a5b98b4..edc48bb26b 100644 --- a/doc/sector.md +++ b/doc/sector.md @@ -35,6 +35,14 @@ Having downloaded the necessary data, ::: build_biomass_potentials +## Rule `build_perennials_yields_eurostat_average` + +::: build_perennials_yields_eurostat_average + +## Rule `build_perennials_yields` + +::: build_perennials_yields + ## Rule `build_egs_potentials` ::: build_egs_potentials From 4f0dc374fc3afa00147b7e2ff76ab8242bd95a95 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 20 Aug 2026 13:07:31 +0000 Subject: [PATCH 08/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- rules/build_sector.smk | 4 +- rules/retrieve.smk | 46 +++++++--- scripts/build_biomass_potentials.py | 6 +- scripts/build_perennials_yields.py | 6 +- ...uild_perennials_yields_eurostat_average.py | 72 ++++++++++------ scripts/prepare_sector_network.py | 83 +++++++++++-------- 6 files changed, 139 insertions(+), 78 deletions(-) diff --git a/rules/build_sector.smk b/rules/build_sector.smk index 240ff5e2bc..6622b6fcb7 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -968,8 +968,6 @@ rule build_perennials_yields_eurostat_average: rule build_perennials_yields: - params: - biomass=config_provider("biomass"), input: nuts2=rules.retrieve_eu_nuts_2021.output.shapes_level_2, country_shapes=resources("country_shapes.geojson"), @@ -981,6 +979,8 @@ rule build_perennials_yields: logs("build_perennials_yields_s_{clusters}.log"), resources: mem_mb=8000, + params: + biomass=config_provider("biomass"), script: scripts("build_perennials_yields.py") diff --git a/rules/retrieve.smk b/rules/retrieve.smk index b47f985ae7..968a1d4049 100755 --- a/rules/retrieve.smk +++ b/rules/retrieve.smk @@ -1637,17 +1637,23 @@ if (CO2_REMOVAL_DATASET := dataset_version("co2_removal_data"))["source"] in [ ]: rule retrieve_co2_removal_data: - message: - "Downloading carbon dioxide removal data (afforestation, perennialisation inputs)" input: zip=storage(CO2_REMOVAL_DATASET["url"]), output: - afforestation_nuts_biomass_densities=resources("afforestation_nuts_biomass_densities.xlsx"), - afforestation_nuts2_afforestation_rates=resources("afforestation_rates_nuts2_full.csv"), - afforestation_nuts2_monthly_weights=resources("afforestation_nuts2_monthly_weights.csv"), + afforestation_nuts_biomass_densities=resources( + "afforestation_nuts_biomass_densities.xlsx" + ), + afforestation_nuts2_afforestation_rates=resources( + "afforestation_rates_nuts2_full.csv" + ), + afforestation_nuts2_monthly_weights=resources( + "afforestation_nuts2_monthly_weights.csv" + ), eurostat_crops_nuts2=resources("eurostat_apro_cpshr_nuts2_raw.csv"), eurostat_crops_nuts0=resources("eurostat_apro_cpshr_nuts0_raw.csv"), retries: 2 + message: + "Downloading carbon dioxide removal data (afforestation, perennialisation inputs)" run: with ZipFile(input.zip) as z: # GitHub's release archive nests everything under a single @@ -1655,11 +1661,29 @@ if (CO2_REMOVAL_DATASET := dataset_version("co2_removal_data"))["source"] in [ # changes with every release, so resolve it at runtime. top_dir = z.namelist()[0].split("/")[0] for src_path, dest in [ - ("outputs/afforestation/afforestation_nuts_biomass_densities.xlsx", output.afforestation_nuts_biomass_densities), - ("outputs/afforestation/afforestation_rates_nuts2_full.csv", output.afforestation_nuts2_afforestation_rates), - ("outputs/afforestation/afforestation_nuts2_monthly_weights.csv", output.afforestation_nuts2_monthly_weights), - ("outputs/perennialisation/eurostat_apro_cpshr_nuts2_raw.csv", output.eurostat_crops_nuts2), - ("outputs/perennialisation/eurostat_apro_cpshr_nuts0_raw.csv", output.eurostat_crops_nuts0), + ( + "outputs/afforestation/afforestation_nuts_biomass_densities.xlsx", + output.afforestation_nuts_biomass_densities, + ), + ( + "outputs/afforestation/afforestation_rates_nuts2_full.csv", + output.afforestation_nuts2_afforestation_rates, + ), + ( + "outputs/afforestation/afforestation_nuts2_monthly_weights.csv", + output.afforestation_nuts2_monthly_weights, + ), + ( + "outputs/perennialisation/eurostat_apro_cpshr_nuts2_raw.csv", + output.eurostat_crops_nuts2, + ), + ( + "outputs/perennialisation/eurostat_apro_cpshr_nuts0_raw.csv", + output.eurostat_crops_nuts0, + ), ]: - with z.open(f"{top_dir}/{src_path}") as src, open(dest, "wb") as dst: + with ( + z.open(f"{top_dir}/{src_path}") as src, + open(dest, "wb") as dst, + ): dst.write(src.read()) diff --git a/scripts/build_biomass_potentials.py b/scripts/build_biomass_potentials.py index bfe643d42e..aa163338cd 100755 --- a/scripts/build_biomass_potentials.py +++ b/scripts/build_biomass_potentials.py @@ -12,7 +12,11 @@ import numpy as np import pandas as pd -from scripts._helpers import configure_logging, resolve_biomass_classes, set_scenario_config +from scripts._helpers import ( + configure_logging, + resolve_biomass_classes, + set_scenario_config, +) logger = logging.getLogger(__name__) AVAILABLE_BIOMASS_YEARS = [2010, 2020, 2030, 2040, 2050] diff --git a/scripts/build_perennials_yields.py b/scripts/build_perennials_yields.py index fec978c9cd..2188321b22 100644 --- a/scripts/build_perennials_yields.py +++ b/scripts/build_perennials_yields.py @@ -91,7 +91,6 @@ def impute_missing_values(df_nuts2, missing_shapes, yield_cols): imputed_rows = [] for missing_id, c_geom in missing_centroids.items(): - # Distance to only VALID NUTS2 rows dists = valid_centroids.distance(c_geom) nearest_nuts2 = dists.idxmin() @@ -220,10 +219,7 @@ def convert_nuts2_to_regions_yields(df_nuts2, regions, yield_cols=None): "perennials", ] imputed_missing = impute_missing_values(df_nuts2, missing_shapes, yield_cols) - df_nuts2 = pd.concat([ - df_nuts2.drop(index=missing_countries), - imputed_missing - ]) + df_nuts2 = pd.concat([df_nuts2.drop(index=missing_countries), imputed_missing]) # convert nuts2 yields to regions df = convert_nuts2_to_regions_yields(df_nuts2, regions) diff --git a/scripts/build_perennials_yields_eurostat_average.py b/scripts/build_perennials_yields_eurostat_average.py index 31e06996d9..3d80818659 100644 --- a/scripts/build_perennials_yields_eurostat_average.py +++ b/scripts/build_perennials_yields_eurostat_average.py @@ -40,9 +40,10 @@ logger = logging.getLogger(__name__) + def harmonize_to_nuts2021(df, keep_col, nuts2021_n2): """ - df : DataFrame indexed by ['geo', 'TIME_PERIOD', 'mapping'] + Df : DataFrame indexed by ['geo', 'TIME_PERIOD', 'mapping'] contains both NUTS2 and NUTS0 rows keep_col : column to harmonize (e.g. 'weighted_YL_(t/ha)') nuts2021_n2 : GeoDataFrame with index=NUTS2_ID and geometry @@ -108,7 +109,10 @@ def harmonize_to_nuts2021(df, keep_col, nuts2021_n2): result.sort_index() return result -def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, biofuel_yields): + +def calculate_yields( + filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, biofuel_yields +): # filter columns - keep only relevant df_crops_raw_nuts2 = pd.read_csv(filepath_nuts2) df_crops_raw_nuts2["TIME_PERIOD"] = df_crops_raw_nuts2["TIME_PERIOD"].astype(int) @@ -116,7 +120,9 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b df_crops_raw_nuts0 = pd.read_csv(filepath_nuts0) df_crops_raw_nuts0["TIME_PERIOD"] = df_crops_raw_nuts0["TIME_PERIOD"].astype(int) - df_crops_raw = pd.concat([df_crops_raw_nuts0, df_crops_raw_nuts2], ignore_index=True) + df_crops_raw = pd.concat( + [df_crops_raw_nuts0, df_crops_raw_nuts2], ignore_index=True + ) # drop empty and irrelevant columns columns_to_drop = [ @@ -146,7 +152,8 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b # Step 1: Filter to relevant rows for 2023 and strucpro of interest df_sub = df_crops[ - (df_crops["strucpro"].isin(["AR", "PR_HU_EU"])) & (df_crops["crops"].isin(crops_sel)) + (df_crops["strucpro"].isin(["AR", "PR_HU_EU"])) + & (df_crops["crops"].isin(crops_sel)) ][["crops", "geo", "TIME_PERIOD", "strucpro", "OBS_VALUE"]] # Pivot so AR and PR_HU_EU are columns for each (crop, geo, year) @@ -169,10 +176,9 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b ) # Average data across years - df_avg_yield = ( - df_pivot.groupby(["crops", "geo"], as_index=False)[["AR", "PR_HU_EU", "YL_(t/ha)"]] - .mean() - ) + df_avg_yield = df_pivot.groupby(["crops", "geo"], as_index=False)[ + ["AR", "PR_HU_EU", "YL_(t/ha)"] + ].mean() min_year = df_pivot["TIME_PERIOD"].min() max_year = df_pivot["TIME_PERIOD"].max() @@ -193,7 +199,9 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b df_avg_yield["PR_share"] = df_avg_yield["PR_share"].fillna(0) # calculated average weighted yield - df_avg_yield["weighted_YL_(t/ha)"] = df_avg_yield["YL_(t/ha)"] * df_avg_yield["PR_share"] + df_avg_yield["weighted_YL_(t/ha)"] = ( + df_avg_yield["YL_(t/ha)"] * df_avg_yield["PR_share"] + ) # sanity check for very low yields due to small productions thresholds = { @@ -205,9 +213,11 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b # Apply crop-specific minimum threshold df_avg_yield["weighted_YL_(t/ha)"] = df_avg_yield.apply( - lambda row: row["weighted_YL_(t/ha)"] - if row["weighted_YL_(t/ha)"] >= thresholds.get(row["mapping"], 0) - else 0, + lambda row: ( + row["weighted_YL_(t/ha)"] + if row["weighted_YL_(t/ha)"] >= thresholds.get(row["mapping"], 0) + else 0 + ), axis=1, ) @@ -225,7 +235,9 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b unsustainable_biofuels_yields["energy_yields_(MWh/ha)"] = ( unsustainable_biofuels_yields["weighted_YL_(t/ha)"] - * unsustainable_biofuels_yields.index.get_level_values("mapping").map(biofuel_yields) + * unsustainable_biofuels_yields.index.get_level_values("mapping").map( + biofuel_yields + ) ) # yields of perennials per hectare in ton/ha @@ -238,7 +250,9 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b ] * (1 - std_moist_perennials) # max yields from current production : applies to perennials for green biorefining - max_yields = df_avg_yield.groupby(["geo", "TIME_PERIOD", "mapping"])["YL_(t/ha)"].max() + max_yields = df_avg_yield.groupby(["geo", "TIME_PERIOD", "mapping"])[ + "YL_(t/ha)" + ].max() perennial_yields_max = pd.DataFrame(max_yields) perennial_yields_max = perennial_yields_max[ @@ -282,14 +296,20 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b ) costs = load_costs(snakemake.input.costs) - #conv = snakemake.params.biofuel_conversion + # conv = snakemake.params.biofuel_conversion # LHV values are fixed physical constants, not parameters: JRC Technical Report doi:10.2760/69179 - LHV_fuels = {"ethanol": 7.447, "biodiesel": 10.194} # MWh/t (26.81 MJ/kg, 36.7 MJ/kg) + LHV_fuels = { + "ethanol": 7.447, + "biodiesel": 10.194, + } # MWh/t (26.81 MJ/kg, 36.7 MJ/kg) biofuel_yields = { - "MINBIOCRP11": costs.at['ethanol from wheat', 'efficiency'] * LHV_fuels["ethanol"], - "MINBIOCRP21": costs.at['ethanol from sugar beet', 'efficiency'] * LHV_fuels["ethanol"], - "MINBIORPS1": costs.at['biodiesel from rapeseed', 'efficiency'] * LHV_fuels["biodiesel"], + "MINBIOCRP11": costs.at["ethanol from wheat", "efficiency"] + * LHV_fuels["ethanol"], + "MINBIOCRP21": costs.at["ethanol from sugar beet", "efficiency"] + * LHV_fuels["ethanol"], + "MINBIORPS1": costs.at["biodiesel from rapeseed", "efficiency"] + * LHV_fuels["biodiesel"], } other_crops_codes = [ @@ -300,12 +320,14 @@ def calculate_yields(filepath_nuts2, filepath_nuts0, crops_sel, crops_mapping, b crops_sel = perennial_codes + other_crops_codes logger.info("Computing crop yields...") - unsustainable_biofuels_yields, perennial_yields, perennial_yields_max = calculate_yields( - filepath_nuts0=CROPS_CSV_NUTS0, - filepath_nuts2=CROPS_CSV_NUTS2, - crops_sel=crops_sel, - crops_mapping=crops_mapping, - biofuel_yields=biofuel_yields, + unsustainable_biofuels_yields, perennial_yields, perennial_yields_max = ( + calculate_yields( + filepath_nuts0=CROPS_CSV_NUTS0, + filepath_nuts2=CROPS_CSV_NUTS2, + crops_sel=crops_sel, + crops_mapping=crops_mapping, + biofuel_yields=biofuel_yields, + ) ) yield_MINBIOCRP11 = unsustainable_biofuels_yields[ diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 6b638e792b..baf47299e0 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -1368,62 +1368,77 @@ def add_perennials(n, costs): # load resources biomass_potentials = pd.read_csv(snakemake.input.biomass_potentials, index_col=0) - perennials_yields_1G_biofuels = pd.read_csv(snakemake.input.perennials_yields_1G_biofuels).set_index("name") + perennials_yields_1G_biofuels = pd.read_csv( + snakemake.input.perennials_yields_1G_biofuels + ).set_index("name") # calculate perennials potential based on the conversion on first generation biofuels for equal area - perennials_area = (biomass_potentials.filter(regex='biofuels_1G') / perennials_yields_1G_biofuels.filter(regex='biofuels_1G')).sum(axis=1) + perennials_area = ( + biomass_potentials.filter(regex="biofuels_1G") + / perennials_yields_1G_biofuels.filter(regex="biofuels_1G") + ).sum(axis=1) # (MWh/y) / (MWh / ha / y) = (ha) returns the area used by sum of the 3 biofuels_1G classes which can be assigned for perennials - perennials_potentials = perennials_area * snakemake.config["perennials"]["sequestration_co2"] # (tCO2seq) = (ha) * (tCO2 seq/ha) + perennials_potentials = ( + perennials_area * snakemake.config["perennials"]["sequestration_co2"] + ) # (tCO2seq) = (ha) * (tCO2 seq/ha) nodes = pop_layout.index n.add("Carrier", "co2 perennials") n.add( - "Bus", - nodes, - suffix=" co2 perennials", - location=nodes, - carrier="co2 perennials", - unit="t_co2", + "Bus", + nodes, + suffix=" co2 perennials", + location=nodes, + carrier="co2 perennials", + unit="t_co2", ) # calculate CO2 sequestration per tDM perennials - perennial_CO2_seq = perennials_yields_1G_biofuels["perennials"] / snakemake.config["perennials"]["sequestration_co2"] # (tDM/tCO2 seq) + perennial_CO2_seq = ( + perennials_yields_1G_biofuels["perennials"] + / snakemake.config["perennials"]["sequestration_co2"] + ) # (tDM/tCO2 seq) # calculate biogas production based on harvesting time (in month) df_harvest = pd.DataFrame(index=n.snapshots, columns=["harvest"]) - df_harvest["harvest"] = df_harvest.index.month.isin([4, 5, 6, 7, 8, 9, 10]).astype(int) + df_harvest["harvest"] = df_harvest.index.month.isin([4, 5, 6, 7, 8, 9, 10]).astype( + int + ) p_max_pu = pd.DataFrame(index=n.snapshots, columns=nodes) for node in nodes: p_max_pu[node] = df_harvest["harvest"] n.add( - "Link", - nodes, - suffix=" perennials refining", - bus0="co2 atmosphere", - bus1=nodes + " co2 perennials", - bus2=nodes.values, - bus3=spatial.gas.biogas, - efficiency=1, - efficiency2=-costs.at["perennials refining", "electricity-input"] * perennial_CO2_seq, - efficiency3=costs.at["perennials refining", "biogas-output"] * perennial_CO2_seq, - carrier="co2 perennials", - p_nom_extendable=True, - p_max_pu=p_max_pu, - capital_cost=costs.at["perennials refining", "capital_cost"] * perennial_CO2_seq, - marginal_cost=costs.at["perennials refining", "VOM"] * perennial_CO2_seq, - lifetime=costs.at["perennials refining", "lifetime"], + "Link", + nodes, + suffix=" perennials refining", + bus0="co2 atmosphere", + bus1=nodes + " co2 perennials", + bus2=nodes.values, + bus3=spatial.gas.biogas, + efficiency=1, + efficiency2=-costs.at["perennials refining", "electricity-input"] + * perennial_CO2_seq, + efficiency3=costs.at["perennials refining", "biogas-output"] + * perennial_CO2_seq, + carrier="co2 perennials", + p_nom_extendable=True, + p_max_pu=p_max_pu, + capital_cost=costs.at["perennials refining", "capital_cost"] + * perennial_CO2_seq, + marginal_cost=costs.at["perennials refining", "VOM"] * perennial_CO2_seq, + lifetime=costs.at["perennials refining", "lifetime"], ) n.add( - "Store", - nodes, - suffix=" CO2s perennials", - bus=nodes + " co2 perennials", - e_nom=perennials_potentials.values, - carrier="co2 perennials", - e_cyclic=False, + "Store", + nodes, + suffix=" CO2s perennials", + bus=nodes + " co2 perennials", + e_nom=perennials_potentials.values, + carrier="co2 perennials", + e_cyclic=False, ) From f6490e36ace7a6542f45a39e7a850b00cb10813b Mon Sep 17 00:00:00 2001 From: BertoGBG Date: Tue, 25 Aug 2026 10:04:20 +0200 Subject: [PATCH 09/10] Fix add_perennials cost lookup: "perennials refining" -> "perennials gbr" The technology-data source (BertoGBG/technology-data pypsa-eur_AA branch) names this row "perennials gbr" (Grass BioRefinery), not "perennials refining". add_perennials() was looking up the wrong name, which would crash with KeyError: 'perennials refining' as soon as a run actually exercised this path. Already fixed on the `perennialisation` branch; ported here and to a_CDRs/heat_industry. --- scripts/prepare_sector_network.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index baf47299e0..cc59543a5d 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -1345,7 +1345,7 @@ def add_perennials(n, costs): The PyPSA network container object costs : pd.DataFrame Costs and parameters for different technologies. Must contain a - 'perennials refining' entry with 'electricity-input', 'biogas-output', + 'perennials gbr' entry with 'electricity-input', 'biogas-output', 'capital_cost', 'VOM', and 'lifetime' parameters Returns @@ -1418,17 +1418,17 @@ def add_perennials(n, costs): bus2=nodes.values, bus3=spatial.gas.biogas, efficiency=1, - efficiency2=-costs.at["perennials refining", "electricity-input"] + efficiency2=-costs.at["perennials gbr", "electricity-input"] * perennial_CO2_seq, - efficiency3=costs.at["perennials refining", "biogas-output"] + efficiency3=costs.at["perennials gbr", "biogas-output"] * perennial_CO2_seq, carrier="co2 perennials", p_nom_extendable=True, p_max_pu=p_max_pu, - capital_cost=costs.at["perennials refining", "capital_cost"] + capital_cost=costs.at["perennials gbr", "capital_cost"] * perennial_CO2_seq, - marginal_cost=costs.at["perennials refining", "VOM"] * perennial_CO2_seq, - lifetime=costs.at["perennials refining", "lifetime"], + marginal_cost=costs.at["perennials gbr", "VOM"] * perennial_CO2_seq, + lifetime=costs.at["perennials gbr", "lifetime"], ) n.add( From 4c85cf584c434efad50a4679ece93ddcb5ea05cb Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 25 Aug 2026 08:04:38 +0000 Subject: [PATCH 10/10] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- scripts/prepare_sector_network.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index cc59543a5d..ae1fc07d7b 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -1420,13 +1420,11 @@ def add_perennials(n, costs): efficiency=1, efficiency2=-costs.at["perennials gbr", "electricity-input"] * perennial_CO2_seq, - efficiency3=costs.at["perennials gbr", "biogas-output"] - * perennial_CO2_seq, + efficiency3=costs.at["perennials gbr", "biogas-output"] * perennial_CO2_seq, carrier="co2 perennials", p_nom_extendable=True, p_max_pu=p_max_pu, - capital_cost=costs.at["perennials gbr", "capital_cost"] - * perennial_CO2_seq, + capital_cost=costs.at["perennials gbr", "capital_cost"] * perennial_CO2_seq, marginal_cost=costs.at["perennials gbr", "VOM"] * perennial_CO2_seq, lifetime=costs.at["perennials gbr", "lifetime"], )