249 lines
9.6 KiB
Python
249 lines
9.6 KiB
Python
#!/usr/bin/env python
|
|
"""
|
|
Create ``.csv`` files and plots for comparing per country full load hours of renewable time series.
|
|
|
|
Relevant Settings
|
|
-----------------
|
|
|
|
.. code:: yaml
|
|
|
|
snapshots:
|
|
|
|
renewable:
|
|
{technology}:
|
|
cutout:
|
|
resource:
|
|
correction_factor:
|
|
|
|
Inputs
|
|
------
|
|
|
|
- ``data/bundle/corine/g250_clc06_V18_5.tif``:
|
|
|
|
.. image:: img/corine.png
|
|
:scale: 33 %
|
|
|
|
- ``data/bundle/GEBCO_2014_2D.nc``:
|
|
- ``data/pietzker2014.xlsx``: NOT IN DATA BUNDLE!
|
|
- ``resources/natura.tiff``: confer :ref:`natura`
|
|
- ``resources/country_shapes.geojson``: confer :ref:`shapes`
|
|
- ``resources/offshore_shapes.geojson``: confer :ref:`shapes`
|
|
- ``resources/regions_onshore.geojson``: (if not offshore wind), confer :ref:`busregions`
|
|
- ``resources/regions_offshore.geojson``: (if offshore wind), :ref:`busregions`
|
|
- ``"cutouts/" + config["renewable"][{technology}]['cutout']``: :ref:`cutout`
|
|
- ``networks/base.nc``: :ref:`base`
|
|
|
|
Outputs
|
|
-------
|
|
|
|
- ``resources/country_flh_area_{technology}.csv``:
|
|
- ``resources/country_flh_aggregated_{technology}.csv``:
|
|
- ``resources/country_flh_uncorrected_{technology}.csv``:
|
|
- ``resources/country_flh_{technology}.pdf``:
|
|
- ``resources/country_exclusion_{technology}``:
|
|
|
|
Description
|
|
-----------
|
|
|
|
"""
|
|
|
|
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 None:
|
|
agg_a['Overall'] = float((np.asarray((areamatrix.sum(axis=0) / areamatrix.sum())
|
|
.dot(flh.values)).squeeze()))
|
|
|
|
agg = pd.DataFrame({'PyPSA-Eur': agg_a})
|
|
else:
|
|
# 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)
|