diff --git a/Snakefile b/Snakefile index 6c6a7c88..3e26051a 100644 --- a/Snakefile +++ b/Snakefile @@ -300,6 +300,30 @@ rule plot_p_nom_max: output: "results/plots/{network}_s{simpl}_cum_p_nom_max_{clusters}_{technology}_{country}.{ext}" script: "scripts/plot_p_nom_max.py" +rule build_country_flh: + input: + base_network="networks/base.nc", + corine="data/bundle/corine/g250_clc06_V18_5.tif", + natura="resources/natura.tiff", + gebco="data/bundle/GEBCO_2014_2D.nc", + country_shapes='resources/country_shapes.geojson', + offshore_shapes='resources/offshore_shapes.geojson', + pietzker="data/pietzker2014.xlsx", + regions=lambda w: ("resources/country_shapes.geojson" + if w.technology in ('onwind', 'solar') + else "resources/offshore_shapes.geojson"), + cutout=lambda w: "cutouts/" + config["renewable"][w.technology]['cutout'] + output: + area="resources/country_flh_area_{technology}.csv", + aggregated="resources/country_flh_aggregated_{technology}.csv", + uncorrected="resources/country_flh_uncorrected_{technology}.csv", + plot="resources/country_flh_{technology}.pdf", + exclusion=directory("resources/country_exclusion_{technology}") + resources: mem=10000 + benchmark: "benchmarks/build_country_flh_{technology}" + # group: 'feedin_preparation' + script: "scripts/build_country_flh.py" + # Local Variables: # mode: python # End: diff --git a/scripts/build_country_flh.py b/scripts/build_country_flh.py new file mode 100644 index 00000000..273e6551 --- /dev/null +++ b/scripts/build_country_flh.py @@ -0,0 +1,197 @@ +#!/usr/bin/env python + +import os +import atlite +import numpy as np +import xarray as xr +import pandas as pd + +import geokit as gk +from osgeo import gdal +from scipy.sparse import vstack +import pycountry as pyc +import matplotlib.pyplot as plt + +from vresutils import landuse as vlanduse +from vresutils.array import spdiag + +import progressbar as pgb +import logging +logger = logging.getLogger(__name__) + +from build_renewable_profiles import init_globals, calculate_potential + +def build_area(flh, countries, areamatrix, breaks, fn): + area_unbinned = xr.DataArray(areamatrix.todense(), [countries, capacity_factor.coords['spatial']]) + bins = xr.DataArray(pd.cut(flh.to_series(), bins=breaks), flh.coords, name="bins") + area = area_unbinned.groupby(bins).sum(dim="spatial").to_pandas() + area.loc[:,slice(*area.sum()[lambda s: s > 0].index[[0,-1]])].to_csv(fn) + area.columns = area.columns.map(lambda s: s.left) + return area + +def plot_area_not_solar(area, countries): + # onshore wind/offshore wind + a = area.T + + fig, axes = plt.subplots(nrows=len(countries), sharex=True) + for c, ax in zip(countries, axes): + d = a[[c]] / 1e3 + d.plot.bar(ax=ax, legend=False, align='edge', width=1.) + ax.set_ylabel(f"Potential {c} / GW") + ax.set_title(c) + ax.legend() + ax.set_xlabel("Full-load hours") + fig.savefig(snakemake.output.plot, transparent=True, bbox_inches='tight') + +def plot_area_solar(area, p_area, countries): + # onshore wind/offshore wind + p = p_area.T + a = area.T + + fig, axes = plt.subplots(nrows=len(countries), sharex=True, squeeze=False) + for c, ax in zip(countries, axes.flat): + d = pd.concat([a[c], p[c]], keys=['PyPSA-Eur', 'Pietzker'], axis=1) / 1e3 + d.plot.bar(ax=ax, legend=False, align='edge', width=1.) + # ax.set_ylabel(f"Potential {c} / GW") + ax.set_title(c) + ax.legend() + ax.set_xlabel("Full-load hours") + + fig.savefig(snakemake.output.plot, transparent=True, bbox_inches='tight') + + +def build_aggregate(flh, countries, areamatrix, breaks, p_area, fn): + agg_a = pd.Series(np.ravel((areamatrix / areamatrix.sum(axis=1)).dot(flh.values)), + countries, name="PyPSA-Eur") + + if p_area is not None: + # Determine indices of countries which are also in Pietzcker + inds = pd.Index(countries).get_indexer(p_area.index) + areamatrix = areamatrix[inds] + + agg_a['Overall'] = float((np.asarray((areamatrix.sum(axis=0) / areamatrix.sum()) + .dot(flh.values)).squeeze())) + + + midpoints = (breaks[1:] + breaks[:-1])/2. + p = p_area.T + + # Per-country FLH comparison + agg_p = pd.Series((p / p.sum()).multiply(midpoints, axis=0).sum(), name="Pietzker") + agg_p['Overall'] = float((p.sum(axis=1) / p.sum().sum()).multiply(midpoints, axis=0).sum()) + + agg = pd.DataFrame({'PyPSA-Eur': agg_a, 'Pietzcker': agg_p, 'Ratio': agg_p / agg_a}) + + agg.to_csv(fn) + +if __name__ == '__main__': + # Detect running outside of snakemake and mock snakemake for testing + if 'snakemake' not in globals(): + from vresutils.snakemake import MockSnakemake, Dict + snakemake = MockSnakemake( + wildcards=Dict(technology='solar'), + input=Dict( + base_network="networks/base.nc", + corine="data/bundle/corine/g250_clc06_V18_5.tif", + natura="resources/natura.tiff", + gebco="data/bundle/GEBCO_2014_2D.nc", + country_shapes='resources/country_shapes.geojson', + offshore_shapes='resources/offshore_shapes.geojson', + pietzker="data/pietzker2014.xlsx" + ), + output=Dict( + area="resources/country_flh_area_{technology}.csv", + aggregated="resources/country_flh_aggregated_{technology}.csv", + uncorrected="resources/country_flh_uncorrected_{technology}.csv", + plot="resources/country_flh_{technology}.pdf", + exclusion="resources/country_exclusion_{technology}" + ) + ) + snakemake.input['regions'] = os.path.join(snakemake.path, "resources", + "country_shapes.geojson" + if snakemake.wildcards.technology in ('onwind', 'solar') + else "offshore_shapes.geojson") + snakemake.input['cutout'] = os.path.join(snakemake.path, "cutouts", + snakemake.config["renewable"][snakemake.wildcards.technology]['cutout']) + + pgb.streams.wrap_stderr() + logging.basicConfig(level=snakemake.config['logging_level']) + + config = snakemake.config['renewable'][snakemake.wildcards.technology] + + time = pd.date_range(freq='m', **snakemake.config['snapshots']) + params = dict(years=slice(*time.year[[0, -1]]), months=slice(*time.month[[0, -1]])) + + cutout = atlite.Cutout(config['cutout'], + cutout_dir=os.path.dirname(snakemake.input.cutout), + **params) + + minx, maxx, miny, maxy = cutout.extent + dx = (maxx - minx) / (cutout.shape[1] - 1) + dy = (maxy - miny) / (cutout.shape[0] - 1) + bounds = gk.Extent.from_xXyY((minx - dx/2., maxx + dx/2., + miny - dy/2., maxy + dy/2.)) + + # Use GLAES to compute available potentials and the transition matrix + paths = dict(snakemake.input) + + init_globals(bounds, dx, dy, config, paths) + regions = gk.vector.extractFeatures(paths["regions"], onlyAttr=True) + countries = pd.Index(regions["name"], name="country") + + widgets = [ + pgb.widgets.Percentage(), + ' ', pgb.widgets.SimpleProgress(format='(%s)' % pgb.widgets.SimpleProgress.DEFAULT_FORMAT), + ' ', pgb.widgets.Bar(), + ' ', pgb.widgets.Timer(), + ' ', pgb.widgets.ETA() + ] + progressbar = pgb.ProgressBar(prefix='Compute GIS potentials: ', widgets=widgets, max_value=len(countries)) + + if not os.path.isdir(snakemake.output.exclusion): + os.makedirs(snakemake.output.exclusion) + + matrix = vstack([calculate_potential(i, save_map=os.path.join(snakemake.output.exclusion, countries[i])) + for i in progressbar(regions.index)]) + + areamatrix = matrix * spdiag(vlanduse._cutout_cell_areas(cutout).ravel()) + areamatrix.data[areamatrix.data < 1.] = 0 # ignore weather cells where only less than 1 km^2 can be installed + areamatrix.eliminate_zeros() + + resource = config['resource'] + func = getattr(cutout, resource.pop('method')) + correction_factor = config.get('correction_factor', 1.) + + capacity_factor = func(capacity_factor=True, show_progress='Compute capacity factors: ', **resource).stack(spatial=('y', 'x')) + flh_uncorr = capacity_factor * 8760 + flh_corr = correction_factor * flh_uncorr + + if snakemake.wildcards.technology == 'solar': + pietzcker = pd.read_excel(snakemake.input.pietzker, sheet_name="PV on all area", skiprows=2, header=[0,1]).iloc[1:177] + p_area1_50 = pietzcker['Usable Area at given FLh in 1-50km distance to settlement '].dropna(axis=1) + p_area1_50.columns = p_area1_50.columns.str.split(' ').str[0] + + p_area50_100 = pietzcker['Usable Area at given FLh in 50-100km distance to settlement '] + + p_area = p_area1_50 + p_area50_100 + cols = p_area.columns + breaks = cols.str.split('-').str[0].append(pd.Index([cols[-1].split('-')[1]])).astype(int) + p_area.columns = breaks[:-1] + + p_area = p_area.reindex(countries.map(lambda c: pyc.countries.get(alpha_2=c).name)) + p_area.index = countries + p_area = p_area.dropna() # Pietzcker does not have data for CZ and MK + else: + breaks = np.r_[0:8000:50] + p_area = None + + + area = build_area(flh_corr, countries, areamatrix, breaks, snakemake.output.area) + + if snakemake.wildcards.technology == 'solar': + plot_area_solar(area, p_area, p_area.index) + else: + plot_area_not_solar(area, countries) + + build_aggregate(flh_uncorr, countries, areamatrix, breaks, p_area, snakemake.output.uncorrected) + build_aggregate(flh_corr, countries, areamatrix, breaks, p_area, snakemake.output.aggregated)