diff --git a/config/config.default.yaml b/config/config.default.yaml index ea95ad6d58..259b799d39 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -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 diff --git a/scripts/build_powerplants.py b/scripts/build_powerplants.py index 364b286356..347fe7638b 100755 --- a/scripts/build_powerplants.py +++ b/scripts/build_powerplants.py @@ -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: @@ -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) ) diff --git a/scripts/build_renewable_profiles.py b/scripts/build_renewable_profiles.py index 724d002e16..3b007c74a5 100644 --- a/scripts/build_renewable_profiles.py +++ b/scripts/build_renewable_profiles.py @@ -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) regions = regions.loc[offshore_regions] regions = regions.map(lambda g: _simplify_polys(g, minarea=1)).set_crs( regions.crs @@ -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) + continue row = layoutmatrix.sel(bus=bus, bin=bin).data nz_b = row != 0 row = row[nz_b] diff --git a/scripts/cluster_network.py b/scripts/cluster_network.py index 9eb9bd2fdf..c2cca8de98 100644 --- a/scripts/cluster_network.py +++ b/scripts/cluster_network.py @@ -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 @@ -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]) diff --git a/scripts/determine_availability_matrix_MD_UA.py b/scripts/determine_availability_matrix_MD_UA.py index 27a091db60..bcab169ad2 100644 --- a/scripts/determine_availability_matrix_MD_UA.py +++ b/scripts/determine_availability_matrix_MD_UA.py @@ -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 @@ -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": @@ -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: @@ -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: diff --git a/scripts/simplify_network.py b/scripts/simplify_network.py index 6ba67f0700..ac5bba10c8 100644 --- a/scripts/simplify_network.py +++ b/scripts/simplify_network.py @@ -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)