Merge branch 'master' into technology-data
This commit is contained in:
lisazeyen
commit
42a78262d8
@ -29,6 +29,9 @@ before_install:
|
||||
# list packages for easier debugging
|
||||
- conda list
|
||||
|
||||
before_script:
|
||||
- 'echo -ne "url: ${CDSAPI_URL}\nkey: ${CDSAPI_TOKEN}\n" > ~/.cdsapirc'
|
||||
|
||||
script:
|
||||
- cp ./test/config.test1.yaml ./config.yaml
|
||||
- snakemake -j all solve_all_networks
|
||||
|
||||
55
Snakefile
55
Snakefile
@ -136,10 +136,12 @@ rule build_bus_regions:
|
||||
resources: mem=1000
|
||||
script: "scripts/build_bus_regions.py"
|
||||
|
||||
|
||||
if config['enable'].get('build_cutout', False):
|
||||
rule build_cutout:
|
||||
output: directory("cutouts/{cutout}")
|
||||
input:
|
||||
regions_onshore="resources/regions_onshore.geojson",
|
||||
regions_offshore="resources/regions_offshore.geojson"
|
||||
output: "cutouts/{cutout}.nc"
|
||||
log: "logs/build_cutout/{cutout}.log"
|
||||
benchmark: "benchmarks/build_cutout_{cutout}"
|
||||
threads: ATLITE_NPROCESSES
|
||||
@ -149,16 +151,16 @@ if config['enable'].get('build_cutout', False):
|
||||
|
||||
if config['enable'].get('retrieve_cutout', True):
|
||||
rule retrieve_cutout:
|
||||
output: directory(expand("cutouts/{cutouts}", **config['atlite'])),
|
||||
log: "logs/retrieve_cutout.log"
|
||||
script: 'scripts/retrieve_cutout.py'
|
||||
input: HTTP.remote("zenodo.org/record/4709858/files/{cutout}.nc", keep_local=True)
|
||||
output: "cutouts/{cutout}.nc"
|
||||
shell: "mv {input} {output}"
|
||||
|
||||
|
||||
if config['enable'].get('build_natura_raster', False):
|
||||
rule build_natura_raster:
|
||||
input:
|
||||
natura="data/bundle/natura/Natura2000_end2015.shp",
|
||||
cutouts=expand("cutouts/{cutouts}", **config['atlite'])
|
||||
cutouts=expand("cutouts/{cutouts}.nc", **config['atlite'])
|
||||
output: "resources/natura.tiff"
|
||||
log: "logs/build_natura_raster.log"
|
||||
script: "scripts/build_natura_raster.py"
|
||||
@ -166,9 +168,9 @@ if config['enable'].get('build_natura_raster', False):
|
||||
|
||||
if config['enable'].get('retrieve_natura_raster', True):
|
||||
rule retrieve_natura_raster:
|
||||
input: HTTP.remote("zenodo.org/record/4706686/files/natura.tiff", keep_local=True)
|
||||
output: "resources/natura.tiff"
|
||||
log: "logs/retrieve_natura_raster.log"
|
||||
script: 'scripts/retrieve_natura_raster.py'
|
||||
shell: "mv {input} {output}"
|
||||
|
||||
if config['enable'].get('retrieve_cost_data', True):
|
||||
rule retrieve_cost_data:
|
||||
@ -190,11 +192,10 @@ rule build_renewable_profiles:
|
||||
country_shapes='resources/country_shapes.geojson',
|
||||
offshore_shapes='resources/offshore_shapes.geojson',
|
||||
regions=lambda w: ("resources/regions_onshore.geojson"
|
||||
if w.technology in ('onwind', 'solar')
|
||||
else "resources/regions_offshore.geojson"),
|
||||
cutout=lambda w: "cutouts/" + config["renewable"][w.technology]['cutout']
|
||||
output:
|
||||
profile="resources/profile_{technology}.nc",
|
||||
if w.technology in ('onwind', 'solar')
|
||||
else "resources/regions_offshore.geojson"),
|
||||
cutout=lambda w: "cutouts/" + config["renewable"][w.technology]['cutout'] + ".nc"
|
||||
output: profile="resources/profile_{technology}.nc",
|
||||
log: "logs/build_renewable_profile_{technology}.log"
|
||||
benchmark: "benchmarks/build_renewable_profiles_{technology}"
|
||||
threads: ATLITE_NPROCESSES
|
||||
@ -207,7 +208,7 @@ if 'hydro' in config['renewable'].keys():
|
||||
input:
|
||||
country_shapes='resources/country_shapes.geojson',
|
||||
eia_hydro_generation='data/bundle/EIA_hydro_generation_2000_2014.csv',
|
||||
cutout="cutouts/" + config["renewable"]['hydro']['cutout']
|
||||
cutout="cutouts/" + config["renewable"]['hydro']['cutout'] + ".nc"
|
||||
output: 'resources/profile_hydro.nc'
|
||||
log: "logs/build_hydro_profile.log"
|
||||
resources: mem=5000
|
||||
@ -397,29 +398,3 @@ rule plot_p_nom_max:
|
||||
log: "logs/plot_p_nom_max/elec_s{simpl}_{clusts}_{techs}_{country}_{ext}.log"
|
||||
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=lambda w: ("data/bundle/GEBCO_2014_2D.nc"
|
||||
if "max_depth" in config["renewable"][w.technology].keys()
|
||||
else []),
|
||||
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}")
|
||||
log: "logs/build_country_flh_{technology}.log"
|
||||
resources: mem=10000
|
||||
benchmark: "benchmarks/build_country_flh_{technology}"
|
||||
script: "scripts/build_country_flh.py"
|
||||
|
||||
@ -63,18 +63,28 @@ electricity:
|
||||
atlite:
|
||||
nprocesses: 4
|
||||
cutouts:
|
||||
# use 'base' to determine geographical bounds and time span from config
|
||||
# base:
|
||||
# module: era5
|
||||
europe-2013-era5:
|
||||
module: era5
|
||||
xs: [-12., 35.]
|
||||
ys: [72., 33.]
|
||||
years: [2013, 2013]
|
||||
module: era5 # in priority order
|
||||
x: [-12., 35.]
|
||||
y: [33., 72]
|
||||
dx: 0.3
|
||||
dy: 0.3
|
||||
time: ['2013', '2013']
|
||||
europe-2013-sarah:
|
||||
module: sarah
|
||||
resolution: 0.2
|
||||
xs: [-12., 42.]
|
||||
ys: [65., 33.]
|
||||
years: [2013, 2013]
|
||||
module: [sarah, era5] # in priority order
|
||||
x: [-12., 45.]
|
||||
y: [33., 65]
|
||||
dx: 0.2
|
||||
dy: 0.2
|
||||
time: ['2013', '2013']
|
||||
sarah_interpolate: false
|
||||
sarah_dir:
|
||||
features: [influx, temperature]
|
||||
|
||||
|
||||
renewable:
|
||||
onwind:
|
||||
cutout: europe-2013-era5
|
||||
|
||||
@ -55,16 +55,15 @@ electricity:
|
||||
atlite:
|
||||
nprocesses: 4
|
||||
cutouts:
|
||||
europe-2013-era5:
|
||||
europe-2013-era5-tutorial:
|
||||
module: era5
|
||||
xs: [4., 15.]
|
||||
ys: [56., 46.]
|
||||
months: [3, 3]
|
||||
years: [2013, 2013]
|
||||
x: [4., 15.]
|
||||
y: [46., 56.]
|
||||
time: ["2013-03", "2013-03"]
|
||||
|
||||
renewable:
|
||||
onwind:
|
||||
cutout: europe-2013-era5
|
||||
cutout: europe-2013-era5-tutorial
|
||||
resource:
|
||||
method: wind
|
||||
turbine: Vestas_V112_3MW
|
||||
@ -81,7 +80,7 @@ renewable:
|
||||
potential: simple # or conservative
|
||||
clip_p_max_pu: 1.e-2
|
||||
offwind-ac:
|
||||
cutout: europe-2013-era5
|
||||
cutout: europe-2013-era5-tutorial
|
||||
resource:
|
||||
method: wind
|
||||
turbine: NREL_ReferenceTurbine_5MW_offshore
|
||||
@ -93,7 +92,7 @@ renewable:
|
||||
potential: simple # or conservative
|
||||
clip_p_max_pu: 1.e-2
|
||||
offwind-dc:
|
||||
cutout: europe-2013-era5
|
||||
cutout: europe-2013-era5-tutorial
|
||||
resource:
|
||||
method: wind
|
||||
turbine: NREL_ReferenceTurbine_5MW_offshore
|
||||
@ -106,7 +105,7 @@ renewable:
|
||||
potential: simple # or conservative
|
||||
clip_p_max_pu: 1.e-2
|
||||
solar:
|
||||
cutout: europe-2013-era5
|
||||
cutout: europe-2013-era5-tutorial
|
||||
resource:
|
||||
method: pv
|
||||
panel: CSi
|
||||
|
||||
@ -1,8 +1,9 @@
|
||||
,Unit,Values,Description
|
||||
nprocesses,--,int,"Number of parallel processes in cutout preparation"
|
||||
cutouts,,,
|
||||
-- {name},--,"Convention is to name cutouts like ``<region>-<year>-<source>`` (e.g. ``europe-2013-era5``).","Directory to write cutout data to. The user may specify multiple cutouts under configuration ``atlite: cutouts:``. Reference is used in configuration ``renewable: {technology}: cutout:``"
|
||||
-- -- module,--,"One of {'era5','sarah'}","Source of the reanalysis weather dataset (e.g. `ERA5 <https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5>`_ or `SARAH-2 <https://wui.cmsaf.eu/safira/action/viewDoiDetails?acronym=SARAH_V002>`_)"
|
||||
-- -- xs,°,"Float interval within [-180, 180]","Range of longitudes to download weather data for."
|
||||
-- -- ys,°,"Float interval within [-90, 90]","Range of latitudes to download weather data for."
|
||||
-- -- years,--,"Integer interval within [1979,2018]","Range of years to download weather data for."
|
||||
-- {name},--,"Convention is to name cutouts like ``<region>-<year>-<source>`` (e.g. ``europe-2013-era5``).","Name of the cutout netcdf file. The user may specify multiple cutouts under configuration ``atlite: cutouts:``. Reference is used in configuration ``renewable: {technology}: cutout:``. The cutout ``base`` may be used to automatically calculate temporal and spatial bounds of the network."
|
||||
-- -- module,--,"Subset of {'era5','sarah'}","Source of the reanalysis weather dataset (e.g. `ERA5 <https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5>`_ or `SARAH-2 <https://wui.cmsaf.eu/safira/action/viewDoiDetails?acronym=SARAH_V002>`_)"
|
||||
-- -- x,°,"Float interval within [-180, 180]","Range of longitudes to download weather data for. If not defined, it defaults to the spatial bounds of all bus shapes."
|
||||
-- -- y,°,"Float interval within [-90, 90]","Range of latitudes to download weather data for. If not defined, it defaults to the spatial bounds of all bus shapes."
|
||||
-- -- time,,"Time interval within ['1979', '2018'] (with valid pandas date time strings)","Time span to download weather data for. If not defined, it defaults to the time interval spanned by the snapshots."
|
||||
-- -- features,,"String or list of strings with valid cutout features ('inlfux', 'wind').","When freshly building a cutout, retrieve data only for those features. If not defined, it defaults to all available features."
|
||||
|
||||
|
@ -95,9 +95,12 @@ Specifies the temporal range to build an energy system model for as arguments to
|
||||
``atlite``
|
||||
=============
|
||||
|
||||
Define and specify the ``atlite.Cutout`` used for calculating renewable potentials and time-series. All options except for ``features`` are directly used as `cutout parameters <https://atlite.readthedocs.io/en/latest/ref_api.html#cutout>`_.
|
||||
|
||||
.. literalinclude:: ../config.default.yaml
|
||||
:language: yaml
|
||||
:lines: 62-75
|
||||
:start-at: atlite:
|
||||
:end-before: renewable:
|
||||
|
||||
.. csv-table::
|
||||
:header-rows: 1
|
||||
|
||||
@ -9,50 +9,6 @@ Plotting and Summary
|
||||
|
||||
.. warning:: The corresponding code is currently under revision and has only minimal documentation.
|
||||
|
||||
.. _flh:
|
||||
|
||||
Rule ``build_country_flh``
|
||||
=============================
|
||||
|
||||
.. graphviz::
|
||||
:align: center
|
||||
|
||||
digraph snakemake_dag {
|
||||
graph [bgcolor=white,
|
||||
margin=0,
|
||||
size="8,5"
|
||||
];
|
||||
node [fontname=sans,
|
||||
fontsize=10,
|
||||
penwidth=2,
|
||||
shape=box,
|
||||
style=rounded
|
||||
];
|
||||
edge [color=grey,
|
||||
penwidth=2
|
||||
];
|
||||
0 [color="0.31 0.6 0.85",
|
||||
fillcolor=gray,
|
||||
label=build_country_flh,
|
||||
style=filled];
|
||||
1 [color="0.06 0.6 0.85",
|
||||
label=base_network];
|
||||
1 -> 0;
|
||||
2 [color="0.42 0.6 0.85",
|
||||
label=build_natura_raster];
|
||||
2 -> 0;
|
||||
3 [color="0.58 0.6 0.85",
|
||||
label=build_shapes];
|
||||
3 -> 0;
|
||||
4 [color="0.14 0.6 0.85",
|
||||
label=build_cutout];
|
||||
4 -> 0;
|
||||
}
|
||||
|
||||
|
|
||||
|
||||
.. automodule:: build_country_flh
|
||||
|
||||
.. _plot_potentials:
|
||||
|
||||
Rule ``plot_p_nom_max``
|
||||
|
||||
@ -21,9 +21,59 @@ Rule ``retrieve_databundle``
|
||||
Rule ``retrieve_cutout``
|
||||
------------------------
|
||||
|
||||
.. automodule:: retrieve_cutout
|
||||
.. image:: https://zenodo.org/badge/DOI/10.5281/zenodo.3517949.svg
|
||||
:target: https://doi.org/10.5281/zenodo.3517949
|
||||
|
||||
Cutouts are spatio-temporal subsets of the European weather data from the `ECMWF ERA5 <https://software.ecmwf.int/wiki/display/CKB/ERA5+data+documentation>`_ reanalysis dataset and the `CMSAF SARAH-2 <https://wui.cmsaf.eu/safira/action/viewDoiDetails?acronym=SARAH_V002>`_ solar surface radiation dataset for the year 2013.
|
||||
They have been prepared by and are for use with the `atlite <https://github.com/PyPSA/atlite>`_ tool. You can either generate them yourself using the ``build_cutouts`` rule or retrieve them directly from `zenodo <https://doi.org/10.5281/zenodo.3517949>`_ through the rule ``retrieve_cutout``.
|
||||
The :ref:`tutorial` uses a smaller cutout than required for the full model (30 MB), which is also automatically downloaded.
|
||||
|
||||
.. note::
|
||||
To download cutouts yourself from the `ECMWF ERA5 <https://software.ecmwf.int/wiki/display/CKB/ERA5+data+documentation>`_ you need to `set up the CDS API <https://cds.climate.copernicus.eu/api-how-to>`_.
|
||||
|
||||
|
||||
**Relevant Settings**
|
||||
|
||||
.. code:: yaml
|
||||
|
||||
tutorial:
|
||||
enable:
|
||||
build_cutout:
|
||||
|
||||
.. seealso::
|
||||
Documentation of the configuration file ``config.yaml`` at
|
||||
:ref:`toplevel_cf`
|
||||
|
||||
**Outputs**
|
||||
|
||||
- ``cutouts/{cutout}``: weather data from either the `ERA5 <https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5>`_ reanalysis weather dataset or `SARAH-2 <https://wui.cmsaf.eu/safira/action/viewProduktSearch>`_ satellite-based historic weather data.
|
||||
|
||||
.. seealso::
|
||||
For details see :mod:`build_cutout` and read the `atlite documentation <https://atlite.readthedocs.io>`_.
|
||||
|
||||
|
||||
Rule ``retrieve_natura_raster``
|
||||
-------------------------------
|
||||
|
||||
.. automodule:: retrieve_natura_raster
|
||||
.. image:: https://zenodo.org/badge/DOI/10.5281/zenodo.4706686.svg
|
||||
:target: https://doi.org/10.5281/zenodo.4706686
|
||||
|
||||
This rule, as a substitute for :mod:`build_natura_raster`, downloads an already rasterized version (`natura.tiff <https://zenodo.org/record/4706686/files/natura.tiff>`_) of `Natura 2000 <https://en.wikipedia.org/wiki/Natura_2000>`_ natural protection areas to reduce computation times. The file is placed into the ``resources`` sub-directory.
|
||||
|
||||
**Relevant Settings**
|
||||
|
||||
.. code:: yaml
|
||||
|
||||
enable:
|
||||
build_natura_raster:
|
||||
|
||||
.. seealso::
|
||||
Documentation of the configuration file ``config.yaml`` at
|
||||
:ref:`toplevel_cf`
|
||||
|
||||
**Outputs**
|
||||
|
||||
- ``resources/natura.tiff``: Rasterized version of `Natura 2000 <https://en.wikipedia.org/wiki/Natura_2000>`_ natural protection areas to reduce computation times.
|
||||
|
||||
.. seealso::
|
||||
For details see :mod:`build_natura_raster`.
|
||||
|
||||
@ -11,8 +11,12 @@ Release Notes
|
||||
Upcoming Release
|
||||
================
|
||||
|
||||
* Switch to new major release, ``>=v0.2.1`` of ``atlite``. The version upgrade comes along with significant speed up for the rule ``build_renewable_profiles.py`` (~factor 2). A lot of the code which calculated the landuse availability is now outsourced and does not rely on ``glaes``, ``geokit`` anymore. This facilitates the environment building and version compatibility of ``gdal``, ``libgdal`` with other packages.
|
||||
* The minimum python version was set to ``3.8``.
|
||||
* The rule and script ``build_country_flh`` are removed as they're no longer used or maintained.
|
||||
* The flag ``keep_all_available_areas`` in the configuration for renewable potentials (config.yaml -> renewable -> {technology}) was deprecated and now defaults to ``True``.
|
||||
* The tutorial cutout was renamed from ``cutouts/europe-2013-era5.nc`` to ``cutouts/europe-2013-era5-tutorial.nc`` to accomodate tutorial and productive cutouts side-by-side.
|
||||
* Fix: Value for ``co2base`` in ``config.yaml`` adjusted to 1.487e9 t CO2-eq (from 3.1e9 t CO2-eq). The new value represents emissions related to the electricity sector for EU+UK. The old value was ~2x too high and used when the emissions wildcard in ``{opts}`` was used.
|
||||
|
||||
* Add option to include marginal costs of links representing fuel cells, electrolysis, and battery inverters
|
||||
[`#232 <https://github.com/PyPSA/pypsa-eur/pull/232>`_].
|
||||
|
||||
|
||||
@ -130,8 +130,7 @@ It can take the values ``onwind``, ``offwind-ac``, ``offwind-dc``, and ``solar``
|
||||
|
||||
The wildcard can moreover be used to create technology specific figures and summaries.
|
||||
For instance ``{technology}`` can be used to plot regionally disaggregated potentials
|
||||
with the rule :mod:`plot_p_nom_max` or to summarize a particular technology's
|
||||
full load hours in various countries with the rule :mod:`build_country_flh`.
|
||||
with the rule :mod:`plot_p_nom_max`.
|
||||
|
||||
.. _attr:
|
||||
|
||||
|
||||
@ -9,7 +9,8 @@ dependencies:
|
||||
- python<=3.7
|
||||
- pip
|
||||
- pypsa>=0.17.1
|
||||
- atlite=0.0.3
|
||||
- atlite>=0.2.2
|
||||
- dask<=2021.3.1 # until https://github.com/dask/dask/issues/7583 is solved
|
||||
- pre-commit
|
||||
|
||||
# Dependencies of the workflow itself
|
||||
@ -19,27 +20,13 @@ dependencies:
|
||||
- memory_profiler
|
||||
- yaml
|
||||
- pytables
|
||||
- powerplantmatching>=0.4.3
|
||||
|
||||
# Second order dependencies which should really be deps of atlite
|
||||
- xarray
|
||||
- progressbar2
|
||||
- pyyaml>=5.1.0
|
||||
- powerplantmatching>=0.4.8
|
||||
|
||||
# GIS dependencies have to come all from conda-forge
|
||||
- cartopy
|
||||
- fiona
|
||||
- proj
|
||||
- pyshp
|
||||
- geopandas
|
||||
- rasterio
|
||||
- shapely
|
||||
- libgdal
|
||||
- descartes
|
||||
|
||||
- pip:
|
||||
- vresutils==0.3.1
|
||||
- git+https://github.com/PyPSA/glaes.git#egg=glaes
|
||||
- git+https://github.com/PyPSA/geokit.git#egg=geokit
|
||||
- cdsapi
|
||||
- sphinx
|
||||
- sphinx_rtd_theme
|
||||
|
||||
@ -8,12 +8,13 @@ channels:
|
||||
- bioconda
|
||||
- http://conda.anaconda.org/gurobi
|
||||
dependencies:
|
||||
- python
|
||||
- python>=3.8
|
||||
- pip
|
||||
- mamba # esp for windows build
|
||||
|
||||
- pypsa>=0.17.1
|
||||
- atlite=0.0.3
|
||||
- pypsa>=0.17.1
|
||||
- atlite>=0.2.2
|
||||
- dask<=2021.3.1 # until https://github.com/dask/dask/issues/7583 is solved
|
||||
|
||||
# Dependencies of the workflow itself
|
||||
- xlrd
|
||||
@ -29,32 +30,14 @@ dependencies:
|
||||
- powerplantmatching>=0.4.8
|
||||
- numpy<=1.19.0 # otherwise macos fails
|
||||
|
||||
# Second order dependencies which should really be deps of atlite
|
||||
- xarray
|
||||
- netcdf4
|
||||
- bottleneck
|
||||
- toolz
|
||||
- dask
|
||||
- progressbar2
|
||||
- pyyaml>=5.1.0
|
||||
|
||||
# Keep in conda environment when calling ipython
|
||||
- ipython
|
||||
|
||||
# GIS dependencies:
|
||||
- cartopy
|
||||
- fiona
|
||||
- proj
|
||||
- pyshp
|
||||
- geopandas
|
||||
- rasterio
|
||||
- shapely
|
||||
- libgdal<=3.0.4
|
||||
- descartes
|
||||
|
||||
- pip:
|
||||
- vresutils==0.3.1
|
||||
- tsam>=1.1.0
|
||||
- git+https://github.com/PyPSA/glaes.git#egg=glaes
|
||||
- git+https://github.com/PyPSA/geokit.git#egg=geokit
|
||||
- cdsapi
|
||||
|
||||
@ -152,7 +152,6 @@ def attach_stores(n, costs):
|
||||
bus1=buses_i,
|
||||
carrier='battery discharger',
|
||||
efficiency=costs.at['battery inverter','efficiency'],
|
||||
capital_cost=costs.at['battery inverter', 'capital_cost'],
|
||||
p_nom_extendable=True,
|
||||
marginal_cost=costs.at["battery inverter", "marginal_cost"])
|
||||
|
||||
|
||||
@ -1,243 +0,0 @@
|
||||
#!/usr/bin/env python
|
||||
|
||||
# SPDX-FileCopyrightText: : 2017-2020 The PyPSA-Eur Authors
|
||||
#
|
||||
# SPDX-License-Identifier: GPL-3.0-or-later
|
||||
|
||||
"""
|
||||
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:
|
||||
|
||||
.. seealso::
|
||||
Documentation of the configuration file ``config.yaml`` at
|
||||
:ref:`snapshots_cf`, :ref:`renewable_cf`
|
||||
|
||||
Inputs
|
||||
------
|
||||
|
||||
- ``data/bundle/corine/g250_clc06_V18_5.tif``: `CORINE Land Cover (CLC) <https://land.copernicus.eu/pan-european/corine-land-cover>`_ inventory on `44 classes <https://wiki.openstreetmap.org/wiki/Corine_Land_Cover#Tagging>`_ of land use (e.g. forests, arable land, industrial, urban areas).
|
||||
|
||||
.. image:: img/corine.png
|
||||
:scale: 33 %
|
||||
|
||||
- ``data/bundle/GEBCO_2014_2D.nc``: A `bathymetric <https://en.wikipedia.org/wiki/Bathymetry>`_ data set with a global terrain model for ocean and land at 15 arc-second intervals by the `General Bathymetric Chart of the Oceans (GEBCO) <https://www.gebco.net/data_and_products/gridded_bathymetry_data/>`_.
|
||||
|
||||
.. image:: img/gebco_2019_grid_image.jpg
|
||||
:scale: 50 %
|
||||
|
||||
**Source:** `GEBCO <https://www.gebco.net/data_and_products/images/gebco_2019_grid_image.jpg>`_
|
||||
|
||||
- ``data/pietzker2014.xlsx``: `Supplementary material 2 <https://ars.els-cdn.com/content/image/1-s2.0-S0306261914008149-mmc2.xlsx>`_ from `Pietzcker et al. <https://doi.org/10.1016/j.apenergy.2014.08.011>`_; not part of the data bundle; download and place here yourself.
|
||||
- ``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 logging
|
||||
from _helpers import configure_logging
|
||||
|
||||
import os
|
||||
import atlite
|
||||
import numpy as np
|
||||
import xarray as xr
|
||||
import pandas as pd
|
||||
|
||||
import geokit as gk
|
||||
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
|
||||
|
||||
from build_renewable_profiles import init_globals, calculate_potential
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
|
||||
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__':
|
||||
if 'snakemake' not in globals():
|
||||
from _helpers import mock_snakemake
|
||||
snakemake = mock_snakemake('build_country_flh', technology='solar')
|
||||
configure_logging(snakemake)
|
||||
|
||||
pgb.streams.wrap_stderr()
|
||||
|
||||
|
||||
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.xXyY, 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)
|
||||
@ -1,4 +1,4 @@
|
||||
# SPDX-FileCopyrightText: : 2017-2020 The PyPSA-Eur Authors
|
||||
# SPDX-FileCopyrightText: : 2017-2021 The PyPSA-Eur Authors
|
||||
#
|
||||
# SPDX-License-Identifier: GPL-3.0-or-later
|
||||
|
||||
@ -92,10 +92,11 @@ Description
|
||||
"""
|
||||
|
||||
import logging
|
||||
import atlite
|
||||
import geopandas as gpd
|
||||
import pandas as pd
|
||||
from _helpers import configure_logging
|
||||
|
||||
import os
|
||||
import atlite
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
@ -106,14 +107,24 @@ if __name__ == "__main__":
|
||||
configure_logging(snakemake)
|
||||
|
||||
cutout_params = snakemake.config['atlite']['cutouts'][snakemake.wildcards.cutout]
|
||||
for p in ('xs', 'ys', 'years', 'months'):
|
||||
if p in cutout_params:
|
||||
cutout_params[p] = slice(*cutout_params[p])
|
||||
|
||||
cutout = atlite.Cutout(snakemake.wildcards.cutout,
|
||||
cutout_dir=os.path.dirname(snakemake.output[0]),
|
||||
**cutout_params)
|
||||
snapshots = pd.date_range(freq='h', **snakemake.config['snapshots'])
|
||||
time = [snapshots[0], snapshots[-1]]
|
||||
cutout_params['time'] = slice(*cutout_params.get('time', time))
|
||||
|
||||
nprocesses = snakemake.config['atlite'].get('nprocesses', 4)
|
||||
if {'x', 'y', 'bounds'}.isdisjoint(cutout_params):
|
||||
# Determine the bounds from bus regions with a buffer of two grid cells
|
||||
onshore = gpd.read_file(snakemake.input.regions_onshore)
|
||||
offshore = gpd.read_file(snakemake.input.regions_offshore)
|
||||
regions = onshore.append(offshore)
|
||||
d = max(cutout_params.get('dx', 0.25), cutout_params.get('dy', 0.25))*2
|
||||
cutout_params['bounds'] = regions.total_bounds + [-d, -d, d, d]
|
||||
elif {'x', 'y'}.issubset(cutout_params):
|
||||
cutout_params['x'] = slice(*cutout_params['x'])
|
||||
cutout_params['y'] = slice(*cutout_params['y'])
|
||||
|
||||
cutout.prepare(nprocesses=nprocesses)
|
||||
|
||||
logging.info(f"Preparing cutout with parameters {cutout_params}.")
|
||||
features = cutout_params.pop('features', None)
|
||||
cutout = atlite.Cutout(snakemake.output[0], **cutout_params)
|
||||
cutout.prepare(features=features)
|
||||
|
||||
@ -62,7 +62,6 @@ Description
|
||||
import logging
|
||||
from _helpers import configure_logging
|
||||
|
||||
import os
|
||||
import atlite
|
||||
import geopandas as gpd
|
||||
from vresutils import hydro as vhydro
|
||||
@ -76,20 +75,21 @@ if __name__ == "__main__":
|
||||
configure_logging(snakemake)
|
||||
|
||||
config = snakemake.config['renewable']['hydro']
|
||||
cutout_dir = os.path.dirname(snakemake.input.cutout)
|
||||
cutout = atlite.Cutout(config['cutout'], cutout_dir=cutout_dir)
|
||||
cutout = atlite.Cutout(snakemake.input.cutout)
|
||||
|
||||
countries = snakemake.config['countries']
|
||||
country_shapes = gpd.read_file(snakemake.input.country_shapes).set_index('name')['geometry'].reindex(countries)
|
||||
country_shapes = (gpd.read_file(snakemake.input.country_shapes)
|
||||
.set_index('name')['geometry'].reindex(countries))
|
||||
country_shapes.index.name = 'countries'
|
||||
|
||||
eia_stats = vhydro.get_eia_annual_hydro_generation(snakemake.input.eia_hydro_generation).reindex(columns=countries)
|
||||
eia_stats = vhydro.get_eia_annual_hydro_generation(
|
||||
snakemake.input.eia_hydro_generation).reindex(columns=countries)
|
||||
inflow = cutout.runoff(shapes=country_shapes,
|
||||
smooth=True,
|
||||
lower_threshold_quantile=True,
|
||||
normalize_using_yearly=eia_stats)
|
||||
|
||||
if 'clip_min_inflow' in config:
|
||||
inflow.values[inflow.values < config['clip_min_inflow']] = 0.
|
||||
inflow = inflow.where(inflow > config['clip_min_inflow'], 0)
|
||||
|
||||
inflow.to_netcdf(snakemake.output[0])
|
||||
|
||||
@ -43,30 +43,49 @@ import logging
|
||||
from _helpers import configure_logging
|
||||
|
||||
import atlite
|
||||
import geokit as gk
|
||||
from pathlib import Path
|
||||
import geopandas as gpd
|
||||
import rasterio as rio
|
||||
from rasterio.features import geometry_mask
|
||||
from rasterio.warp import transform_bounds
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
|
||||
def determine_cutout_xXyY(cutout_name):
|
||||
cutout = atlite.Cutout(cutout_name, cutout_dir=cutout_dir)
|
||||
cutout = atlite.Cutout(cutout_name)
|
||||
assert cutout.crs.to_epsg() == 4326
|
||||
x, X, y, Y = cutout.extent
|
||||
dx = (X - x) / (cutout.shape[1] - 1)
|
||||
dy = (Y - y) / (cutout.shape[0] - 1)
|
||||
dx, dy = cutout.dx, cutout.dy
|
||||
return [x - dx/2., X + dx/2., y - dy/2., Y + dy/2.]
|
||||
|
||||
|
||||
def get_transform_and_shape(bounds, res):
|
||||
left, bottom = [(b // res)* res for b in bounds[:2]]
|
||||
right, top = [(b // res + 1) * res for b in bounds[2:]]
|
||||
shape = int((top - bottom) // res), int((right - left) / res)
|
||||
transform = rio.Affine(res, 0, left, 0, -res, top)
|
||||
return transform, shape
|
||||
|
||||
|
||||
if __name__ == "__main__":
|
||||
if 'snakemake' not in globals():
|
||||
from _helpers import mock_snakemake
|
||||
snakemake = mock_snakemake('build_natura_raster')
|
||||
configure_logging(snakemake)
|
||||
|
||||
cutout_dir = Path(snakemake.input.cutouts[0]).parent.resolve()
|
||||
cutout_names = {res['cutout'] for res in snakemake.config['renewable'].values()}
|
||||
xs, Xs, ys, Ys = zip(*(determine_cutout_xXyY(cutout) for cutout in cutout_names))
|
||||
xXyY = min(xs), max(Xs), min(ys), max(Ys)
|
||||
|
||||
natura = gk.vector.loadVector(snakemake.input.natura)
|
||||
extent = gk.Extent.from_xXyY(xXyY).castTo(3035).fit(100)
|
||||
extent.rasterize(natura, pixelWidth=100, pixelHeight=100, output=snakemake.output[0])
|
||||
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(snakemake.input.natura).to_crs(3035)
|
||||
raster = ~geometry_mask(shapes.geometry, out_shape[::-1], transform)
|
||||
raster = raster.astype(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)
|
||||
|
||||
|
||||
@ -60,7 +60,6 @@ Inputs
|
||||
**Source:** `GEBCO <https://www.gebco.net/data_and_products/images/gebco_2019_grid_image.jpg>`_
|
||||
|
||||
- ``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`
|
||||
@ -180,220 +179,154 @@ node (`p_nom_max`): ``simple`` and ``conservative``:
|
||||
reached.
|
||||
|
||||
"""
|
||||
import progressbar as pgb
|
||||
import geopandas as gpd
|
||||
import xarray as xr
|
||||
import numpy as np
|
||||
import atlite
|
||||
import logging
|
||||
from pypsa.geo import haversine
|
||||
from shapely.geometry import LineString
|
||||
import time
|
||||
|
||||
from _helpers import configure_logging
|
||||
|
||||
import os
|
||||
import atlite
|
||||
|
||||
import numpy as np
|
||||
import xarray as xr
|
||||
import pandas as pd
|
||||
import multiprocessing as mp
|
||||
import matplotlib.pyplot as plt
|
||||
import progressbar as pgb
|
||||
|
||||
from scipy.sparse import csr_matrix, vstack
|
||||
from pypsa.geo import haversine
|
||||
from vresutils import landuse as vlanduse
|
||||
from vresutils.array import spdiag
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
bounds = dx = dy = config = paths = gebco = clc = natura = None
|
||||
|
||||
|
||||
def init_globals(bounds_xXyY, n_dx, n_dy, n_config, n_paths):
|
||||
# Late import so that the GDAL Context is only created in the new processes
|
||||
global gl, gk, gdal
|
||||
import glaes as gl
|
||||
import geokit as gk
|
||||
from osgeo import gdal as gdal
|
||||
|
||||
# global in each process of the multiprocessing.Pool
|
||||
global bounds, dx, dy, config, paths, gebco, clc, natura
|
||||
|
||||
bounds = gk.Extent.from_xXyY(bounds_xXyY)
|
||||
dx = n_dx
|
||||
dy = n_dy
|
||||
config = n_config
|
||||
paths = n_paths
|
||||
|
||||
if "max_depth" in config:
|
||||
gebco = gk.raster.loadRaster(paths["gebco"])
|
||||
gebco.SetProjection(gk.srs.loadSRS(4326).ExportToWkt())
|
||||
|
||||
clc = gk.raster.loadRaster(paths["corine"])
|
||||
clc.SetProjection(gk.srs.loadSRS(3035).ExportToWkt())
|
||||
|
||||
natura = gk.raster.loadRaster(paths["natura"])
|
||||
|
||||
|
||||
def downsample_to_coarse_grid(bounds, dx, dy, mask, data):
|
||||
# The GDAL warp function with the 'average' resample algorithm needs a band of zero values of at least
|
||||
# the size of one coarse cell around the original raster or it produces erroneous results
|
||||
orig = mask.createRaster(data=data)
|
||||
padded_extent = mask.extent.castTo(bounds.srs).pad(max(dx, dy)).castTo(mask.srs)
|
||||
padded = padded_extent.fit((mask.pixelWidth, mask.pixelHeight)).warp(orig, mask.pixelWidth, mask.pixelHeight)
|
||||
orig = None # free original raster
|
||||
average = bounds.createRaster(dx, dy, dtype=gdal.GDT_Float32)
|
||||
assert gdal.Warp(average, padded, resampleAlg='average') == 1, "gdal warp failed: %s" % gdal.GetLastErrorMsg()
|
||||
return average
|
||||
|
||||
|
||||
def calculate_potential(gid, save_map=None):
|
||||
feature = gk.vector.extractFeature(paths["regions"], where=gid)
|
||||
ec = gl.ExclusionCalculator(feature.geom)
|
||||
|
||||
corine = config.get("corine", {})
|
||||
if isinstance(corine, list):
|
||||
corine = {'grid_codes': corine}
|
||||
if "grid_codes" in corine:
|
||||
ec.excludeRasterType(clc, value=corine["grid_codes"], invert=True)
|
||||
if corine.get("distance", 0.) > 0.:
|
||||
ec.excludeRasterType(clc, value=corine["distance_grid_codes"], buffer=corine["distance"])
|
||||
|
||||
if config.get("natura", False):
|
||||
ec.excludeRasterType(natura, value=1)
|
||||
if "max_depth" in config:
|
||||
ec.excludeRasterType(gebco, (None, -config["max_depth"]))
|
||||
|
||||
# TODO compute a distance field as a raster beforehand
|
||||
if 'max_shore_distance' in config:
|
||||
ec.excludeVectorType(paths["country_shapes"], buffer=config['max_shore_distance'], invert=True)
|
||||
if 'min_shore_distance' in config:
|
||||
ec.excludeVectorType(paths["country_shapes"], buffer=config['min_shore_distance'])
|
||||
|
||||
if save_map is not None:
|
||||
ec.draw()
|
||||
plt.savefig(save_map, transparent=True)
|
||||
plt.close()
|
||||
|
||||
availability = downsample_to_coarse_grid(bounds, dx, dy, ec.region, np.where(ec.region.mask, ec._availability, 0))
|
||||
|
||||
return csr_matrix(gk.raster.extractMatrix(availability).flatten() / 100.)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
if 'snakemake' not in globals():
|
||||
from _helpers import mock_snakemake
|
||||
snakemake = mock_snakemake('build_renewable_profiles', technology='solar')
|
||||
configure_logging(snakemake)
|
||||
|
||||
pgb.streams.wrap_stderr()
|
||||
|
||||
paths = snakemake.input
|
||||
nprocesses = snakemake.config['atlite'].get('nprocesses')
|
||||
noprogress = not snakemake.config['atlite'].get('show_progress', True)
|
||||
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_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)
|
||||
|
||||
# Use the following for testing the default windows method on linux
|
||||
# mp.set_start_method('spawn')
|
||||
with mp.Pool(initializer=init_globals, initargs=(bounds_xXyY, dx, dy, config, paths),
|
||||
maxtasksperchild=20, processes=snakemake.config['atlite'].get('nprocesses', 2)) as pool:
|
||||
|
||||
# The GDAL library creates a GDAL context on module import, which may not be shared over multiple
|
||||
# processes or the PROJ4 library has a hickup, so we import only after forking.
|
||||
import geokit as gk
|
||||
|
||||
regions = gk.vector.extractFeatures(paths["regions"], onlyAttr=True)
|
||||
buses = pd.Index(regions['name'], name="bus")
|
||||
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(regions))
|
||||
matrix = vstack(list(progressbar(pool.imap(calculate_potential, regions.index))))
|
||||
|
||||
potentials = config['capacity_per_sqkm'] * vlanduse._cutout_cell_areas(cutout)
|
||||
potmatrix = matrix * spdiag(potentials.ravel())
|
||||
if not config.get('keep_all_available_areas', False):
|
||||
potmatrix.data[potmatrix.data < 1.] = 0 # ignore weather cells where only less than 1 MW can be installed
|
||||
potmatrix.eliminate_zeros()
|
||||
|
||||
resource = config['resource']
|
||||
func = getattr(cutout, resource.pop('method'))
|
||||
resource = config['resource'] # pv panel config / wind turbine config
|
||||
correction_factor = config.get('correction_factor', 1.)
|
||||
if correction_factor != 1.:
|
||||
logger.warning('correction_factor is set as {}'.format(correction_factor))
|
||||
capacity_factor = correction_factor * func(capacity_factor=True, show_progress='Compute capacity factors: ', **resource).stack(spatial=('y', 'x')).values
|
||||
layoutmatrix = potmatrix * spdiag(capacity_factor)
|
||||
|
||||
profile, capacities = func(matrix=layoutmatrix, index=buses, per_unit=True,
|
||||
return_capacity=True, show_progress='Compute profiles: ',
|
||||
**resource)
|
||||
|
||||
capacity_per_sqkm = config['capacity_per_sqkm']
|
||||
p_nom_max_meth = config.get('potential', 'conservative')
|
||||
|
||||
if p_nom_max_meth == 'simple':
|
||||
p_nom_max = xr.DataArray(np.asarray(potmatrix.sum(axis=1)).squeeze(), [buses])
|
||||
elif p_nom_max_meth == 'conservative':
|
||||
# p_nom_max has to be calculated for each bus and is the minimal ratio
|
||||
# (min over all weather grid cells of the bus region) between the available
|
||||
# potential (potmatrix) and the used normalised layout (layoutmatrix /
|
||||
# capacities), so we would like to calculate i.e. potmatrix / (layoutmatrix /
|
||||
# capacities). Since layoutmatrix = potmatrix * capacity_factor, this
|
||||
# corresponds to capacities/max(capacity factor in the voronoi cell)
|
||||
p_nom_max = xr.DataArray([1./np.max(capacity_factor[inds]) if len(inds) else 0.
|
||||
for inds in np.split(potmatrix.indices, potmatrix.indptr[1:-1])], [buses]) * capacities
|
||||
if isinstance(config.get("corine", {}), list):
|
||||
config['corine'] = {'grid_codes': config['corine']}
|
||||
|
||||
if correction_factor != 1.:
|
||||
logger.info(f'correction_factor is set as {correction_factor}')
|
||||
|
||||
|
||||
cutout = atlite.Cutout(paths['cutout'])
|
||||
regions = gpd.read_file(paths.regions).set_index('name').rename_axis('bus')
|
||||
buses = regions.index
|
||||
|
||||
excluder = atlite.ExclusionContainer(crs=3035, res=100)
|
||||
|
||||
if config['natura']:
|
||||
excluder.add_raster(paths.natura, nodata=0, allow_no_overlap=True)
|
||||
|
||||
corine = config.get("corine", {})
|
||||
if "grid_codes" in corine:
|
||||
codes = corine["grid_codes"]
|
||||
excluder.add_raster(paths.corine, codes=codes, invert=True, crs=3035)
|
||||
if corine.get("distance", 0.) > 0.:
|
||||
codes = corine["distance_grid_codes"]
|
||||
buffer = corine["distance"]
|
||||
excluder.add_raster(paths.corine, codes=codes, buffer=buffer, crs=3035)
|
||||
|
||||
if "max_depth" in config:
|
||||
func = lambda v: v <= -config['max_depth']
|
||||
excluder.add_raster(paths.gebco, codes=func, crs=4236, nodata=-1000)
|
||||
|
||||
if 'min_shore_distance' in config:
|
||||
buffer = config['min_shore_distance']
|
||||
excluder.add_geometry(paths.country_shapes, buffer=buffer)
|
||||
|
||||
if 'max_shore_distance' in config:
|
||||
buffer = config['max_shore_distance']
|
||||
excluder.add_geometry(paths.country_shapes, buffer=buffer, invert=True)
|
||||
|
||||
kwargs = dict(nprocesses=nprocesses, disable_progressbar=noprogress)
|
||||
if noprogress:
|
||||
logger.info('Calculate landuse availabilities...')
|
||||
start = time.time()
|
||||
availability = cutout.availabilitymatrix(regions, excluder, **kwargs)
|
||||
duration = time.time() - start
|
||||
logger.info(f'Completed availability calculation ({duration:2.2f}s)')
|
||||
else:
|
||||
raise AssertionError('Config key `potential` should be one of "simple" (default) or "conservative",'
|
||||
' not "{}"'.format(p_nom_max_meth))
|
||||
availability = cutout.availabilitymatrix(regions, excluder, **kwargs)
|
||||
|
||||
layout = xr.DataArray(np.asarray(potmatrix.sum(axis=0)).reshape(cutout.shape),
|
||||
[cutout.meta.indexes[ax] for ax in ['y', 'x']])
|
||||
area = cutout.grid.to_crs(3035).area / 1e6
|
||||
area = xr.DataArray(area.values.reshape(cutout.shape),
|
||||
[cutout.coords['y'], cutout.coords['x']])
|
||||
|
||||
# Determine weighted average distance from substation
|
||||
cell_coords = cutout.grid_coordinates()
|
||||
potential = capacity_per_sqkm * availability.sum('bus') * area
|
||||
func = getattr(cutout, resource.pop('method'))
|
||||
resource['dask_kwargs'] = {'num_workers': nprocesses}
|
||||
capacity_factor = correction_factor * func(capacity_factor=True, **resource)
|
||||
layout = capacity_factor * area * capacity_per_sqkm
|
||||
profile, capacities = func(matrix=availability.stack(spatial=['y','x']),
|
||||
layout=layout, index=buses,
|
||||
per_unit=True, return_capacity=True, **resource)
|
||||
|
||||
logger.info(f"Calculating maximal capacity per bus (method '{p_nom_max_meth}')")
|
||||
if p_nom_max_meth == 'simple':
|
||||
p_nom_max = capacity_per_sqkm * availability @ area
|
||||
elif p_nom_max_meth == 'conservative':
|
||||
max_cap_factor = capacity_factor.where(availability!=0).max(['x', 'y'])
|
||||
p_nom_max = capacities / max_cap_factor
|
||||
else:
|
||||
raise AssertionError('Config key `potential` should be one of "simple" '
|
||||
f'(default) or "conservative", not "{p_nom_max_meth}"')
|
||||
|
||||
|
||||
|
||||
logger.info('Calculate average distances.')
|
||||
layoutmatrix = (layout * availability).stack(spatial=['y','x'])
|
||||
|
||||
coords = cutout.grid[['x', 'y']]
|
||||
bus_coords = regions[['x', 'y']]
|
||||
|
||||
average_distance = []
|
||||
for i in regions.index:
|
||||
row = layoutmatrix[i]
|
||||
distances = haversine(regions.loc[i, ['x', 'y']], cell_coords[row.indices])[0]
|
||||
average_distance.append((distances * (row.data / row.data.sum())).sum())
|
||||
centre_of_mass = []
|
||||
for bus in buses:
|
||||
row = layoutmatrix.sel(bus=bus).data
|
||||
nz_b = row != 0
|
||||
row = row[nz_b]
|
||||
co = coords[nz_b]
|
||||
distances = haversine(bus_coords.loc[bus], co)
|
||||
average_distance.append((distances * (row / row.sum())).sum())
|
||||
centre_of_mass.append(co.values.T @ (row / row.sum()))
|
||||
|
||||
average_distance = xr.DataArray(average_distance, [buses])
|
||||
centre_of_mass = xr.DataArray(centre_of_mass, [buses, ('spatial', ['x', 'y'])])
|
||||
|
||||
|
||||
ds = xr.merge([(correction_factor * profile).rename('profile'),
|
||||
capacities.rename('weight'),
|
||||
p_nom_max.rename('p_nom_max'),
|
||||
layout.rename('potential'),
|
||||
average_distance.rename('average_distance')])
|
||||
capacities.rename('weight'),
|
||||
p_nom_max.rename('p_nom_max'),
|
||||
potential.rename('potential'),
|
||||
average_distance.rename('average_distance')])
|
||||
|
||||
|
||||
if snakemake.wildcards.technology.startswith("offwind"):
|
||||
import geopandas as gpd
|
||||
from shapely.geometry import LineString
|
||||
|
||||
offshore_shape = gpd.read_file(snakemake.input.offshore_shapes).unary_union
|
||||
logger.info('Calculate underwater fraction of connections.')
|
||||
offshore_shape = gpd.read_file(paths['offshore_shapes']).unary_union
|
||||
underwater_fraction = []
|
||||
for i in regions.index:
|
||||
row = layoutmatrix[i]
|
||||
centre_of_mass = (cell_coords[row.indices] * (row.data / row.data.sum())[:,np.newaxis]).sum(axis=0)
|
||||
line = LineString([centre_of_mass, regions.loc[i, ['x', 'y']]])
|
||||
underwater_fraction.append(line.intersection(offshore_shape).length / line.length)
|
||||
for bus in buses:
|
||||
p = centre_of_mass.sel(bus=bus).data
|
||||
line = LineString([p, regions.loc[bus, ['x', 'y']]])
|
||||
frac = line.intersection(offshore_shape).length/line.length
|
||||
underwater_fraction.append(frac)
|
||||
|
||||
ds['underwater_fraction'] = xr.DataArray(underwater_fraction, [buses])
|
||||
|
||||
# select only buses with some capacity and minimal capacity factor
|
||||
ds = ds.sel(bus=((ds['profile'].mean('time') > config.get('min_p_max_pu', 0.)) &
|
||||
(ds['p_nom_max'] > config.get('min_p_nom_max', 0.))))
|
||||
(ds['p_nom_max'] > config.get('min_p_nom_max', 0.))))
|
||||
|
||||
if 'clip_p_max_pu' in config:
|
||||
ds['profile'].values[ds['profile'].values < config['clip_p_max_pu']] = 0.
|
||||
min_p_max_pu = config['clip_p_max_pu']
|
||||
ds['profile'] = ds['profile'].where(ds['profile'] >= min_p_max_pu, 0)
|
||||
|
||||
ds.to_netcdf(snakemake.output.profile)
|
||||
|
||||
@ -1,75 +0,0 @@
|
||||
# SPDX-FileCopyrightText: 2019-2020 Fabian Hofmann (FIAS)
|
||||
#
|
||||
# SPDX-License-Identifier: GPL-3.0-or-later
|
||||
|
||||
"""
|
||||
.. image:: https://zenodo.org/badge/DOI/10.5281/zenodo.3517949.svg
|
||||
:target: https://doi.org/10.5281/zenodo.3517949
|
||||
|
||||
Cutouts are spatiotemporal subsets of the European weather data from the `ECMWF ERA5 <https://software.ecmwf.int/wiki/display/CKB/ERA5+data+documentation>`_ reanalysis dataset and the `CMSAF SARAH-2 <https://wui.cmsaf.eu/safira/action/viewDoiDetails?acronym=SARAH_V002>`_ solar surface radiation dataset for the year 2013 (3.9 GB).
|
||||
They have been prepared by and are for use with the `atlite <https://github.com/PyPSA/atlite>`_ tool. You can either generate them yourself using the ``build_cutouts`` rule or retrieve them directly from `zenodo <https://doi.org/10.5281/zenodo.3517949>`_ through the rule ``retrieve_cutout`` described here.
|
||||
|
||||
.. note::
|
||||
To download cutouts yourself from the `ECMWF ERA5 <https://software.ecmwf.int/wiki/display/CKB/ERA5+data+documentation>`_ you need to `set up the CDS API <https://cds.climate.copernicus.eu/api-how-to>`_.
|
||||
|
||||
The :ref:`tutorial` uses smaller `cutouts <https://zenodo.org/record/3518020/files/pypsa-eur-tutorial-cutouts.tar.xz>`_ than required for the full model (19 MB)
|
||||
|
||||
.. image:: https://zenodo.org/badge/DOI/10.5281/zenodo.3518020.svg
|
||||
:target: https://doi.org/10.5281/zenodo.3518020
|
||||
|
||||
|
||||
**Relevant Settings**
|
||||
|
||||
.. code:: yaml
|
||||
|
||||
tutorial:
|
||||
enable:
|
||||
build_cutout:
|
||||
|
||||
.. seealso::
|
||||
Documentation of the configuration file ``config.yaml`` at
|
||||
:ref:`toplevel_cf`
|
||||
|
||||
**Outputs**
|
||||
|
||||
- ``cutouts/{cutout}``: weather data from either the `ERA5 <https://www.ecmwf.int/en/forecasts/datasets/reanalysis-datasets/era5>`_ reanalysis weather dataset or `SARAH-2 <https://wui.cmsaf.eu/safira/action/viewProduktSearch>`_ satellite-based historic weather data.
|
||||
|
||||
.. seealso::
|
||||
For details see :mod:`build_cutout` and read the `atlite documentation <https://atlite.readthedocs.io>`_.
|
||||
|
||||
"""
|
||||
|
||||
import logging
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
from pathlib import Path
|
||||
import tarfile
|
||||
from _helpers import progress_retrieve, configure_logging
|
||||
|
||||
if __name__ == "__main__":
|
||||
if 'snakemake' not in globals():
|
||||
from _helpers import mock_snakemake
|
||||
snakemake = mock_snakemake('retrieve_cutout')
|
||||
rootpath = '..'
|
||||
else:
|
||||
rootpath = '.'
|
||||
|
||||
configure_logging(snakemake) # TODO Make logging compatible with progressbar (see PR #102)
|
||||
|
||||
if snakemake.config['tutorial']:
|
||||
url = "https://zenodo.org/record/3518020/files/pypsa-eur-tutorial-cutouts.tar.xz"
|
||||
else:
|
||||
url = "https://zenodo.org/record/3517949/files/pypsa-eur-cutouts.tar.xz"
|
||||
|
||||
# Save location
|
||||
tarball_fn = Path(f"{rootpath}/cutouts.tar.xz")
|
||||
|
||||
logger.info(f"Downloading cutouts from '{url}'.")
|
||||
progress_retrieve(url, tarball_fn)
|
||||
|
||||
logger.info(f"Extracting cutouts.")
|
||||
tarfile.open(tarball_fn).extractall(path=rootpath)
|
||||
|
||||
tarball_fn.unlink()
|
||||
|
||||
logger.info(f"Cutouts available in '{Path(tarball_fn.stem).stem}'.")
|
||||
@ -1,49 +0,0 @@
|
||||
# Copyright 2019-2020 Fabian Hofmann (FIAS)
|
||||
# SPDX-FileCopyrightText: : 2017-2020 The PyPSA-Eur Authors
|
||||
#
|
||||
# SPDX-License-Identifier: GPL-3.0-or-later
|
||||
|
||||
"""
|
||||
.. image:: https://zenodo.org/badge/DOI/10.5281/zenodo.3518215.svg
|
||||
:target: https://doi.org/10.5281/zenodo.3518215
|
||||
|
||||
This rule, as a substitute for :mod:`build_natura_raster`, downloads an already rasterized version (`natura.tiff <https://zenodo.org/record/3518215/files/natura.tiff>`_) of `Natura 2000 <https://en.wikipedia.org/wiki/Natura_2000>`_ natural protection areas to reduce computation times. The file is placed into the ``resources`` sub-directory.
|
||||
|
||||
**Relevant Settings**
|
||||
|
||||
.. code:: yaml
|
||||
|
||||
enable:
|
||||
build_natura_raster:
|
||||
|
||||
.. seealso::
|
||||
Documentation of the configuration file ``config.yaml`` at
|
||||
:ref:`toplevel_cf`
|
||||
|
||||
**Outputs**
|
||||
|
||||
- ``resources/natura.tiff``: Rasterized version of `Natura 2000 <https://en.wikipedia.org/wiki/Natura_2000>`_ natural protection areas to reduce computation times.
|
||||
|
||||
.. seealso::
|
||||
For details see :mod:`build_natura_raster`.
|
||||
|
||||
"""
|
||||
|
||||
import logging
|
||||
|
||||
from _helpers import progress_retrieve, configure_logging
|
||||
|
||||
logger = logging.getLogger(__name__)
|
||||
|
||||
if __name__ == "__main__":
|
||||
if 'snakemake' not in globals():
|
||||
from _helpers import mock_snakemake
|
||||
snakemake = mock_snakemake('retrieve_natura_raster')
|
||||
configure_logging(snakemake) # TODO Make logging compatible with progressbar (see PR #102)
|
||||
|
||||
url = "https://zenodo.org/record/3518215/files/natura.tiff"
|
||||
|
||||
logger.info(f"Downloading natura raster from '{url}'.")
|
||||
progress_retrieve(url, snakemake.output[0])
|
||||
|
||||
logger.info(f"Natura raster available as '{snakemake.output[0]}'.")
|
||||
@ -54,16 +54,15 @@ electricity:
|
||||
atlite:
|
||||
nprocesses: 4
|
||||
cutouts:
|
||||
europe-2013-era5:
|
||||
europe-2013-era5-tutorial:
|
||||
module: era5
|
||||
xs: [4., 15.]
|
||||
ys: [56., 46.]
|
||||
months: [3, 3]
|
||||
years: [2013, 2013]
|
||||
x: [4., 15.]
|
||||
y: [46., 56.]
|
||||
time: ["2013-03", "2013-03"]
|
||||
|
||||
renewable:
|
||||
onwind:
|
||||
cutout: europe-2013-era5
|
||||
cutout: europe-2013-era5-tutorial
|
||||
resource:
|
||||
method: wind
|
||||
turbine: Vestas_V112_3MW
|
||||
@ -80,7 +79,7 @@ renewable:
|
||||
potential: simple # or conservative
|
||||
clip_p_max_pu: 1.e-2
|
||||
offwind-ac:
|
||||
cutout: europe-2013-era5
|
||||
cutout: europe-2013-era5-tutorial
|
||||
resource:
|
||||
method: wind
|
||||
turbine: NREL_ReferenceTurbine_5MW_offshore
|
||||
@ -92,7 +91,7 @@ renewable:
|
||||
potential: simple # or conservative
|
||||
clip_p_max_pu: 1.e-2
|
||||
offwind-dc:
|
||||
cutout: europe-2013-era5
|
||||
cutout: europe-2013-era5-tutorial
|
||||
resource:
|
||||
method: wind
|
||||
turbine: NREL_ReferenceTurbine_5MW_offshore
|
||||
@ -105,7 +104,7 @@ renewable:
|
||||
potential: simple # or conservative
|
||||
clip_p_max_pu: 1.e-2
|
||||
solar:
|
||||
cutout: europe-2013-era5
|
||||
cutout: europe-2013-era5-tutorial
|
||||
resource:
|
||||
method: pv
|
||||
panel: CSi
|
||||
@ -148,9 +147,9 @@ transformers:
|
||||
|
||||
load:
|
||||
url: https://data.open-power-system-data.org/time_series/2019-06-05/time_series_60min_singleindex.csv
|
||||
power_statistics: True # only for files from <2019; set false in order to get ENTSOE transparency data
|
||||
power_statistics: True # only for files from <2019; set false in order to get ENTSOE transparency data
|
||||
interpolate_limit: 3 # data gaps up until this size are interpolated linearly
|
||||
time_shift_for_large_gaps: 1w # data gaps up until this size are copied by copying from
|
||||
time_shift_for_large_gaps: 1w # data gaps up until this size are copied by copying from
|
||||
manual_adjustments: true # false
|
||||
scaling_factor: 1.0
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user