Skip to content

Fix corine bool and offshore keyerror #2113

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
MarcoAnarmo wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Fix corine bool and offshore keyerror #2113

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
ff70a33
fix: handle corine=false in determine_availability_matrix_MD_UA
MarcoAnarmo Mar 4, 2026
999711d
fix: skip offshore buses missing from clustered regions
MarcoAnarmo Mar 4, 2026
3b86799
fix: add to_380 option to skip 380kV voltage-level collapsing
MarcoAnarmo Mar 19, 2026
ac6c09b
fix: preserve voltage levels and country info during clustering
MarcoAnarmo Mar 19, 2026
2fb49f5
fix: support non-aggregated bus names in powerplant assignment
MarcoAnarmo Mar 19, 2026
9704e1c
fix: handle empty MD/UA regions in availability matrix computation
MarcoAnarmo Mar 19, 2026
fdb9918
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Mar 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions config/config.default.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1100,6 +1100,7 @@ clustering:
remove_islands: false
aggregate_to_tyndp: false
simplify_network:
to_380: true
to_substations: false
remove_stubs: true
remove_stubs_across_borders: false
Expand Down
24 changes: 17 additions & 7 deletions scripts/build_powerplants.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -175,12 +175,16 @@ def map_to_country_bus(
unmatched = []

for country, plants in ppl.groupby("Country"):
country_regions = regions[regions.index.str[:2] == country]
joined = (
plants.sjoin(country_regions)
.rename(columns={"name": "bus"})
.reindex(plants.index)
if "country" in regions.columns:
country_regions = regions[regions["country"] == country]
else:
country_regions = regions[regions.index.str[:2] == country]
joined = plants.sjoin(country_regions[["geometry"]]).rename(
columns={"name": "bus"}
)
# Drop duplicate matches (plant in overlapping onshore/offshore regions)
joined = joined[~joined.index.duplicated(keep="first")]
joined = joined.reindex(plants.index)
assigned.append(joined.dropna(subset=["bus"]))
missing = joined[joined["bus"].isna()]
if not missing.empty:
Expand All @@ -189,10 +193,16 @@ def map_to_country_bus(
if unmatched:
unmatched = pd.concat(unmatched)
for country, plants in unmatched.groupby("Country"):
country_regions = regions[regions.index.str[:2] == country]
if "country" in regions.columns:
country_regions = regions[regions["country"] == country]
else:
country_regions = regions[regions.index.str[:2] == country]
nearest = (
plants.to_crs(3035)
.sjoin_nearest(country_regions.to_crs(3035), max_distance=max_distance)
.sjoin_nearest(
country_regions[["geometry"]].to_crs(3035),
max_distance=max_distance,
)
.rename(columns={"name": "bus"})
.to_crs(4326)
)
Expand Down
4 changes: 4 additions & 0 deletions scripts/build_renewable_profiles.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,6 +157,7 @@
if snakemake.wildcards.technology.startswith("offwind"):
# for offshore regions, the shortest distance to the shoreline is used
offshore_regions = availability.coords["bus"].values
offshore_regions = np.intersect1d(offshore_regions, regions.index)
Copy link

Copilot AI Mar 4, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Filtering offshore_regions only when selecting regions leaves availability (and later resource_regions/class_masks) still containing buses that were dropped during clustering. This can produce profiles/potentials for buses that have no corresponding offshore distance geometry and makes outputs inconsistent (e.g. class_regions written only for the filtered bus set). Consider subsetting availability to the same offshore_regions intersection here (and similarly filter resource_regions after loading) so all downstream computations operate on a consistent bus index.

Suggested change
offshore_regions = np.intersect1d(offshore_regions, regions.index)
offshore_regions = np.intersect1d(offshore_regions, regions.index)
# align availability with the filtered offshore regions to keep a consistent bus index
availability = availability.sel(bus=offshore_regions)

Copilot uses AI. Check for mistakes.
regions = regions.loc[offshore_regions]
regions = regions.map(lambda g: _simplify_polys(g, minarea=1)).set_crs(
regions.crs
Expand Down Expand Up @@ -294,6 +295,9 @@
average_distance = []
bus_bins = layoutmatrix.indexes["bus_bin"]
for bus, bin in bus_bins:
if bus not in regions.index:
average_distance.append(0.0)
Comment on lines 297 to +299
Copy link

Copilot AI Mar 4, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Appending 0.0 for buses missing from regions.index silently produces a plausible-looking distance and can mask data issues (and leave inconsistent outputs if those buses still have nonzero potential). Prefer to ensure those buses are removed upstream (e.g. by subsetting availability/class_masks to the regions bus set), or use NaN plus a clear warning/error so downstream consumers don’t interpret missing geometry as zero distance.

Suggested change
for bus, bin in bus_bins:
if bus not in regions.index:
average_distance.append(0.0)
# Identify buses present in the layout but missing geometry in regions
missing_buses = {bus for bus, _ in bus_bins} - set(regions.index)
if missing_buses:
logger.warning(
"Average distance could not be computed for %d bus(es) missing from "
"`regions.index`: %s. Distances for these buses will be set to NaN.",
len(missing_buses),
", ".join(sorted(map(str, missing_buses))),
)
for bus, bin in bus_bins:
if bus not in regions.index:
# Mark missing geometry explicitly instead of assuming zero distance
average_distance.append(np.nan)

Copilot uses AI. Check for mistakes.
continue
row = layoutmatrix.sel(bus=bus, bin=bin).data
nz_b = row != 0
row = row[nz_b]
Expand Down
6 changes: 4 additions & 2 deletions scripts/cluster_network.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -425,7 +425,7 @@ def clustering_for_n_clusters(
busmap,
bus_strategies=bus_strategies,
line_strategies=line_strategies,
custom_line_groupers=["build_year"],
custom_line_groupers=["build_year", "v_nom"],
)

return clustering
Expand Down Expand Up @@ -703,7 +703,9 @@ def update_bus_coordinates(
# nc.shapes = n.shapes.copy()
for which in ["regions_onshore", "regions_offshore"]:
regions = gpd.read_file(snakemake.input[which])
clustered_regions = cluster_regions((clustering.busmap,), regions)
clustered_regions = cluster_regions(
(clustering.busmap,), regions, with_country=True
)
clustered_regions.to_file(snakemake.output[which])
# append_bus_shapes(nc, clustered_regions, type=which.split("_")[1])

Expand Down
58 changes: 42 additions & 16 deletions scripts/determine_availability_matrix_MD_UA.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@
import fiona
import geopandas as gpd
import numpy as np
import xarray as xr

from scripts._helpers import configure_logging, load_cutout, set_scenario_config

Expand Down Expand Up @@ -57,6 +58,8 @@ def get_wdpa_layer_name(wdpa_fn, layer_substring):
excluder = atlite.ExclusionContainer(crs=3035, res=100)

corine = config.get("corine", {})
if not isinstance(corine, dict):
corine = {}
if "grid_codes" in corine:
# Land cover codes to emulate CORINE results
if snakemake.wildcards.technology == "solar":
Expand Down Expand Up @@ -85,18 +88,35 @@ def get_wdpa_layer_name(wdpa_fn, layer_substring):
snakemake.input.copernicus, codes=codes, buffer=buffer, crs="EPSG:4326"
)

has_valid_geometry = not regions.empty and not regions.geometry.isna().all()

if not has_valid_geometry:
logger.info("No valid MD/UA regions; writing empty availability matrix.")
availability = xr.DataArray(
np.zeros((len(buses), len(cutout.data.y), len(cutout.data.x))),
dims=["bus", "y", "x"],
coords={"bus": buses, "y": cutout.data.y, "x": cutout.data.x},
)
availability.to_netcdf(snakemake.output.availability_matrix)
raise SystemExit(0)

if config["natura"]:
wdpa_fn = (
snakemake.input.wdpa_marine
if "offwind" in snakemake.wildcards.technology
else snakemake.input.wdpa
)
layer = get_wdpa_layer_name(wdpa_fn, "polygons")
wdpa = gpd.read_file(
wdpa_fn,
bbox=regions.geometry,
layer=layer,
).to_crs(3035)

if has_valid_geometry:
layer = get_wdpa_layer_name(wdpa_fn, "polygons")
wdpa = gpd.read_file(
wdpa_fn,
bbox=regions.geometry,
layer=layer,
).to_crs(3035)
else:
logger.info("No valid region geometries for MD/UA; skipping WDPA polygons.")
wdpa = gpd.GeoDataFrame()

# temporary file needed for parallelization
with NamedTemporaryFile(suffix=".geojson", delete=False) as f:
Expand All @@ -107,17 +127,23 @@ def get_wdpa_layer_name(wdpa_fn, layer_substring):
time.sleep(1)
excluder.add_geometry(plg_tmp_fn)

layer = get_wdpa_layer_name(wdpa_fn, "points")
wdpa_pts = gpd.read_file(
wdpa_fn,
bbox=regions.geometry,
layer=layer,
).to_crs(3035)
if has_valid_geometry:
layer = get_wdpa_layer_name(wdpa_fn, "points")
wdpa_pts = gpd.read_file(
wdpa_fn,
bbox=regions.geometry,
layer=layer,
).to_crs(3035)
else:
logger.info("No valid region geometries for MD/UA; skipping WDPA points.")
wdpa_pts = gpd.GeoDataFrame(columns=["REP_AREA", "geometry"])

wdpa_pts = wdpa_pts[wdpa_pts["REP_AREA"] > 1]
wdpa_pts["buffer_radius"] = np.sqrt(wdpa_pts["REP_AREA"] / np.pi) * 1000
wdpa_pts = wdpa_pts.set_geometry(
wdpa_pts["geometry"].buffer(wdpa_pts["buffer_radius"])
)
if not wdpa_pts.empty:
wdpa_pts["buffer_radius"] = np.sqrt(wdpa_pts["REP_AREA"] / np.pi) * 1000
wdpa_pts = wdpa_pts.set_geometry(
wdpa_pts["geometry"].buffer(wdpa_pts["buffer_radius"])
)

# temporary file needed for parallelization
with NamedTemporaryFile(suffix=".geojson", delete=False) as f:
Expand Down
11 changes: 9 additions & 2 deletions scripts/simplify_network.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -432,8 +432,15 @@ def remove_converters(n: pypsa.Network) -> pypsa.Network:
Nyears = n.snapshot_weightings.objective.sum() / 8760
buses_prev, lines_prev, links_prev = len(n.buses), len(n.lines), len(n.links)

linetype_380 = snakemake.config["lines"]["types"][380]
n, trafo_map = simplify_network_to_380(n, linetype_380)
if params.simplify_network.get("to_380", True):
linetype_380 = snakemake.config["lines"]["types"][380]
n, trafo_map = simplify_network_to_380(n, linetype_380)
else:
logger.info(
"Skipping 380kV simplification — preserving original voltage levels"
" and transformers"
)
trafo_map = n.buses.index.to_series()
busmaps = [trafo_map]

n, converter_map = remove_converters(n)
Expand Down
Loading