diff --git a/scripts/build_bus_regions.py b/scripts/build_bus_regions.py index d91d0575..8003d370 100644 --- a/scripts/build_bus_regions.py +++ b/scripts/build_bus_regions.py @@ -47,9 +47,10 @@ from _helpers import configure_logging import pypsa import os import pandas as pd +import numpy as np import geopandas as gpd - -from vresutils.graph import voronoi_partition_pts +from shapely.geometry import Polygon +from scipy.spatial import Voronoi logger = logging.getLogger(__name__) @@ -61,6 +62,53 @@ def save_to_geojson(s, fn): s.to_file(fn, driver='GeoJSON', schema=schema) +def voronoi_partition_pts(points, outline): + """ + Compute the polygons of a voronoi partition of `points` within the + polygon `outline`. Taken from + https://github.com/FRESNA/vresutils/blob/master/vresutils/graph.py + Attributes + ---------- + points : Nx2 - ndarray[dtype=float] + outline : Polygon + Returns + ------- + polygons : N - ndarray[dtype=Polygon|MultiPolygon] + """ + + points = np.asarray(points) + + if len(points) == 1: + polygons = [outline] + else: + xmin, ymin = np.amin(points, axis=0) + xmax, ymax = np.amax(points, axis=0) + xspan = xmax - xmin + yspan = ymax - ymin + + # to avoid any network positions outside all Voronoi cells, append + # the corners of a rectangle framing these points + vor = Voronoi(np.vstack((points, + [[xmin-3.*xspan, ymin-3.*yspan], + [xmin-3.*xspan, ymax+3.*yspan], + [xmax+3.*xspan, ymin-3.*yspan], + [xmax+3.*xspan, ymax+3.*yspan]]))) + + polygons = [] + for i in range(len(points)): + poly = Polygon(vor.vertices[vor.regions[vor.point_region[i]]]) + + if not poly.is_valid: + poly = poly.buffer(0) + + poly = poly.intersection(outline) + + polygons.append(poly) + + + return np.array(polygons, dtype=object) + + if __name__ == "__main__": if 'snakemake' not in globals(): from _helpers import mock_snakemake diff --git a/scripts/build_natura_raster.py b/scripts/build_natura_raster.py index 71d2c45e..7fa9d544 100644 --- a/scripts/build_natura_raster.py +++ b/scripts/build_natura_raster.py @@ -40,7 +40,7 @@ Description """ import logging -from _helpers import configure_logging, retrieve_snakemake_keys +from _helpers import configure_logging import atlite import geopandas as gpd @@ -73,20 +73,17 @@ if __name__ == "__main__": snakemake = mock_snakemake('build_natura_raster') configure_logging(snakemake) - paths, config, wildcards, logs, out = retrieve_snakemake_keys(snakemake) - - cutouts = paths.cutouts + cutouts = snakemake.input.cutouts xs, Xs, ys, Ys = zip(*(determine_cutout_xXyY(cutout) for cutout in cutouts)) bounds = transform_bounds(4326, 3035, min(xs), min(ys), max(Xs), max(Ys)) transform, out_shape = get_transform_and_shape(bounds, res=100) # adjusted boundaries - shapes = gpd.read_file(paths.natura).to_crs(3035) + shapes = gpd.read_file(snakemake.input.natura).to_crs(3035) raster = ~geometry_mask(shapes.geometry, out_shape[::-1], transform) raster = raster.astype(rio.uint8) - with rio.open(out[0], 'w', driver='GTiff', dtype=rio.uint8, + with rio.open(snakemake.output[0], 'w', driver='GTiff', dtype=rio.uint8, count=1, transform=transform, crs=3035, compress='lzw', width=raster.shape[1], height=raster.shape[0]) as dst: dst.write(raster, indexes=1) - diff --git a/scripts/plot_p_nom_max.py b/scripts/plot_p_nom_max.py index ea66d612..e79ad274 100644 --- a/scripts/plot_p_nom_max.py +++ b/scripts/plot_p_nom_max.py @@ -19,7 +19,7 @@ Description """ import logging -from _helpers import configure_logging, retrieve_snakemake_keys +from _helpers import configure_logging import pypsa import pandas as pd @@ -53,13 +53,11 @@ if __name__ == "__main__": clusts= '5,full', country= 'all') configure_logging(snakemake) - paths, config, wildcards, logs, out = retrieve_snakemake_keys(snakemake) - plot_kwds = dict(drawstyle="steps-post") - clusters = wildcards.clusts.split(',') - techs = wildcards.techs.split(',') - country = wildcards.country + clusters = snakemake.wildcards.clusts.split(',') + techs = snakemake.wildcards.techs.split(',') + country = snakemake.wildcards.country if country == 'all': country = None else: @@ -68,7 +66,7 @@ if __name__ == "__main__": fig, axes = plt.subplots(1, len(techs)) for j, cluster in enumerate(clusters): - net = pypsa.Network(paths[j]) + net = pypsa.Network(snakemake.input[j]) for i, tech in enumerate(techs): cum_p_nom_max(net, tech, country).plot(x="p_max_pu", y="cum_p_nom_max", @@ -81,4 +79,4 @@ if __name__ == "__main__": plt.legend(title="Cluster level") - fig.savefig(out[0], transparent=True, bbox_inches='tight') + fig.savefig(snakemake.output[0], transparent=True, bbox_inches='tight') diff --git a/scripts/plot_summary.py b/scripts/plot_summary.py index 48f064b0..bc2bd30c 100644 --- a/scripts/plot_summary.py +++ b/scripts/plot_summary.py @@ -21,7 +21,7 @@ Description import os import logging -from _helpers import configure_logging, retrieve_snakemake_keys +from _helpers import configure_logging import pandas as pd import matplotlib.pyplot as plt @@ -170,12 +170,12 @@ if __name__ == "__main__": attr='', ext='png', country='all') configure_logging(snakemake) - paths, config, wildcards, logs, out = retrieve_snakemake_keys(snakemake) + config = snakemake.config - summary = wildcards.summary + summary = snakemake.wildcards.summary try: func = globals()[f"plot_{summary}"] except KeyError: raise RuntimeError(f"plotting function for {summary} has not been defined") - func(os.path.join(paths[0], f"{summary}.csv"), config, out[0]) + func(os.path.join(snakemake.input[0], f"{summary}.csv"), config, snakemake.output[0]) diff --git a/scripts/prepare_links_p_nom.py b/scripts/prepare_links_p_nom.py index 6bd4bca4..b83089d6 100644 --- a/scripts/prepare_links_p_nom.py +++ b/scripts/prepare_links_p_nom.py @@ -37,7 +37,7 @@ Description """ import logging -from _helpers import configure_logging, retrieve_snakemake_keys +from _helpers import configure_logging import pandas as pd @@ -63,8 +63,6 @@ if __name__ == "__main__": snakemake = mock_snakemake('prepare_links_p_nom', simpl='', network='elec') configure_logging(snakemake) - paths, config, wildcards, logs, out = retrieve_snakemake_keys(snakemake) - links_p_nom = pd.read_html('https://en.wikipedia.org/wiki/List_of_HVDC_projects', header=0, match="SwePol")[0] mw = "Power (MW)" @@ -76,4 +74,4 @@ if __name__ == "__main__": links_p_nom['x1'], links_p_nom['y1'] = extract_coordinates(links_p_nom['Converterstation 1']) links_p_nom['x2'], links_p_nom['y2'] = extract_coordinates(links_p_nom['Converterstation 2']) - links_p_nom.dropna(subset=['x1', 'y1', 'x2', 'y2']).to_csv(out[0], index=False) + links_p_nom.dropna(subset=['x1', 'y1', 'x2', 'y2']).to_csv(snakemake.output[0], index=False)