diff --git a/Snakefile b/Snakefile index 860f6cee..0c9d22a9 100644 --- a/Snakefile +++ b/Snakefile @@ -89,38 +89,31 @@ rule build_bus_regions: rule build_cutout: output: "cutouts/{cutout}" - resources: mem=5000 + resources: mem=config['atlite'].get('nprocesses', 4) * 1000 threads: config['atlite'].get('nprocesses', 4) benchmark: "benchmarks/build_cutout_{cutout}" # group: 'feedin_preparation' script: "scripts/build_cutout.py" -def memory_build_renewable_potentials(wildcards): - corine_config = config["renewable"][wildcards.technology]["corine"] - return 12000 if corine_config.get("distance") is None else 24000 - -rule build_renewable_potentials: - input: - cutout=lambda wildcards: "cutouts/" + config["renewable"][wildcards.technology]['cutout'], - corine="data/bundle/corine/g250_clc06_V18_5.tif", - natura="data/bundle/natura/Natura2000_end2015.shp" - output: "resources/potentials_{technology}.nc" - resources: mem=memory_build_renewable_potentials - benchmark: "benchmarks/build_renewable_potentials_{technology}" - # group: 'feedin_preparation' - script: "scripts/build_renewable_potentials.py" +rule build_natura_raster: + input: "data/bundle/natura/Natura2000_end2015.shp" + output: "resources/natura.tiff" + script: "scripts/build_natura_raster.py" rule build_renewable_profiles: input: base_network="networks/base.nc", - potentials="resources/potentials_{technology}.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', regions=lambda wildcards: ("resources/regions_onshore.geojson" if wildcards.technology in ('onwind', 'solar') else "resources/regions_offshore.geojson"), cutout=lambda wildcards: "cutouts/" + config["renewable"][wildcards.technology]['cutout'] - output: - profile="resources/profile_{technology}.nc", - resources: mem=5000 + output: profile="resources/profile_{technology}.nc", + resources: mem=config['atlite'].get('nprocesses', 2) * 5000 + threads: config['atlite'].get('nprocesses', 2) benchmark: "benchmarks/build_renewable_profiles_{technology}" # group: 'feedin_preparation' script: "scripts/build_renewable_profiles.py" diff --git a/config.yaml b/config.yaml index 10ade504..e39ab16a 100644 --- a/config.yaml +++ b/config.yaml @@ -58,7 +58,7 @@ renewable: method: wind turbine: Vestas_V112_3MW # ScholzPhd Tab 4.3.1: 10MW/km^2 - capacity_per_sqm: 3 + capacity_per_sqkm: 3 # correction_factor: 0.93 corine: #The selection of CORINE Land Cover [1] types that are allowed for wind and solar are based on [2] p.42 / p.28 @@ -70,19 +70,20 @@ renewable: 24, 25, 26, 27, 28, 29, 31, 32] distance: 1000 distance_grid_codes: [1, 2, 3, 4, 5, 6] - natura: true + natura: true + potential: conservative # or heuristic offwind: cutout: europe-2013-era5 resource: method: wind turbine: NREL_ReferenceTurbine_5MW_offshore # ScholzPhd Tab 4.3.1: 10MW/km^2 - capacity_per_sqm: 3 - height_cutoff: 50 + capacity_per_sqkm: 3 # correction_factor: 0.93 - corine: - grid_codes: [44, 255] - natura: true + corine: [44, 255] + natura: true + max_depth: 50 + potential: conservative # or heuristic solar: cutout: europe-2013-sarah resource: @@ -92,12 +93,12 @@ renewable: slope: 35. azimuth: 180. # ScholzPhd Tab 4.3.1: 170 MW/km^2 - capacity_per_sqm: 1.7 + capacity_per_sqkm: 1.7 correction_factor: 0.877 - corine: - grid_codes: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, - 18, 19, 20, 26, 31, 32] - natura: true + corine: [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, + 14, 15, 16, 17, 18, 19, 20, 26, 31, 32] + natura: true + potential: conservative # or heuristic hydro: cutout: europe-2013-era5 carriers: [ror, PHS, hydro] diff --git a/environment.yaml b/environment.yaml index 12e652d2..f7a092d5 100644 --- a/environment.yaml +++ b/environment.yaml @@ -40,6 +40,8 @@ dependencies: - pypsa>=0.13 - vresutils>=0.2.5 - git+https://github.com/FRESNA/atlite.git + - git+https://github.com/FZJ-IEK3-VSA/glaes.git + - git+https://github.com/FZJ-IEK3-VSA/geokit.git #- git+https://github.com/FRESNA/powerplantmatching.git #- https://software.ecmwf.int/wiki/download/attachments/56664858/ecmwf-api-client-python.tgz diff --git a/scripts/build_natura_raster.py b/scripts/build_natura_raster.py new file mode 100644 index 00000000..b3a01f98 --- /dev/null +++ b/scripts/build_natura_raster.py @@ -0,0 +1,18 @@ +import atlite +from osgeo import gdal +import geokit as gk + +def determine_cutout_xXyY(cutout_name): + cutout = atlite.Cutout(cutout_name, cutout_dir="../cutouts") + x, X, y, Y = cutout.extent + dx = (X - x) / (cutout.shape[1] - 1) + dy = (Y - y) / (cutout.shape[0] - 1) + return [x - dx/2., X + dx/2., y - dy/2., Y + dy/2.] + +cutout_names = np.unique([res['cutout'] for res in config['renewable'].values()]) +xs, Xs, ys, Ys = zip(*(determine_cutout_xyXY(cutout) for cutout in cutout_names)) +xXyY = min(xs), max(Xs), min(ys), max(Ys) + +natura = gk.vector.loadVector(snakemake.input[0]) +extent = gk.Extent.from_xXyY(xXyY).castTo(3035).fit(100) +extent.rasterize(natura, pixelWidth=100, pixelHeight=100, output=snakemake.output[0]) diff --git a/scripts/build_renewable_profiles.py b/scripts/build_renewable_profiles.py index 1fe6391e..7abd4202 100644 --- a/scripts/build_renewable_profiles.py +++ b/scripts/build_renewable_profiles.py @@ -5,51 +5,158 @@ import atlite import numpy as np import xarray as xr import pandas as pd -import geopandas as gpd +from multiprocessing import Pool +import glaes as gl +import geokit as gk +from osgeo import gdal +from scipy.sparse import csr_matrix, vstack + +from vresutils import landuse as vlanduse +from vresutils.array import spdiag + +import progressbar as pgb import logging logger = logging.getLogger(__name__) -logging.basicConfig(level=snakemake.config['logging_level']) -config = snakemake.config['renewable'][snakemake.wildcards.technology] +bounds = dx = dy = gebco = clc = natura = None +def init_globals(n_bounds, n_dx, n_dy): + # global in each process of the multiprocessing.Pool + global bounds, dx, dy, gebco, clc, natura -time = pd.date_range(freq='m', **snakemake.config['snapshots']) -params = dict(years=slice(*time.year[[0, -1]]), months=slice(*time.month[[0, -1]])) + bounds = n_bounds + dx = n_dx + dy = n_dy -regions = gpd.read_file(snakemake.input.regions).set_index('name') -regions.index.name = 'bus' + gebco = gk.raster.loadRaster(snakemake.input.gebco) + gebco.SetProjection(gk.srs.loadSRS(4326).ExportToWkt()) -cutout = atlite.Cutout(config['cutout'], - cutout_dir=os.path.dirname(snakemake.input.cutout), - **params) + clc = gk.raster.loadRaster(snakemake.input.corine) + clc.SetProjection(gk.srs.loadSRS(3035).ExportToWkt()) -# Potentials -potentials = xr.open_dataarray(snakemake.input.potentials) + natura = gk.raster.loadRaster(snakemake.input.natura) -# Indicatormatrix -indicatormatrix = cutout.indicatormatrix(regions.geometry) +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 -resource = config['resource'] -func = getattr(cutout, resource.pop('method')) -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, **resource) -layout = capacity_factor * potentials +def calculate_potential(gid): + feature = gk.vector.extractFeature(snakemake.input.regions, where=gid) + ec = gl.ExclusionCalculator(feature.geom, srs=4326, pixelRes=1e-3) # about 100m in EU, it also works in LAEA 3035 -profile, capacities = func(matrix=indicatormatrix, index=regions.index, - layout=layout, per_unit=True, return_capacity=True, - **resource) + 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"]) -relativepotentials = (potentials / layout).stack(spatial=('y', 'x')).values -p_nom_max = xr.DataArray([np.nanmin(relativepotentials[row.nonzero()[1]]) - if row.getnnz() > 0 else 0 - for row in indicatormatrix.tocsr()], - [capacities.coords['bus']]) * capacities + if config.get("natura", False): + ec.excludeRasterType(natura, value=1) + if "max_depth" in config: + ec.excludeRasterType(gebco, (None, -config["max_depth"])) -ds = xr.merge([(correction_factor * profile).rename('profile'), - capacities.rename('weight'), - p_nom_max.rename('p_nom_max'), - layout.rename('potential')]) -(ds.sel(bus=ds['profile'].mean('time') > config.get('min_p_max_pu', 0.)) - .to_netcdf(snakemake.output.profile)) + # TODO compute a distance field as a raster beforehand + if 'max_shore_distance' in config: + ec.excludeVectorType(snakemake.input.country_shapes, buffer=config['max_shore_distance'], invert=True) + if 'min_shore_distance' in config: + ec.excludeVectorType(snakemake.input.country_shapes, buffer=config['min_shore_distance']) + + 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__': + 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 + + with Pool(initializer=init_globals, initargs=(bounds, dx, dy), + maxtasksperchild=20, processes=snakemake.config['atlite'].get('nprocesses', 2)) as pool: + features = gk.vector.extractFeatures(snakemake.input.regions, onlyAttr=True) #.iloc[:10] + buses = pd.Index(features['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(features)) + matrix = vstack(list(progressbar(pool.imap(calculate_potential, features.index)))) + + potentials = config['capacity_per_sqkm'] * vlanduse._cutout_cell_areas(cutout) + potmatrix = matrix * spdiag(potentials.ravel()) + potmatrix.data[potmatrix.data < 1.] = 0 # ignore weather cells where only less than 1 MW can be installed + potmatrix.eliminate_zeros() + + with pgb.ProgressBar(prefix='Compute capacity factors: ', max_value=1) as progressbar: + progressbar.update(0) + resource = config['resource'] + func = getattr(cutout, resource.pop('method')) + 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, **resource).stack(spatial=('y', 'x')).values + layoutmatrix = potmatrix * spdiag(capacity_factor) + progressbar.update(1) + + + with pgb.ProgressBar(prefix='Compute profiles: ', max_value=1) as progressbar: + progressbar.update(0) + profile, capacities = func(matrix=layoutmatrix, index=buses, per_unit=True, + return_capacity=True, **resource) + progressbar.update(1) + + + p_nom_max_meth = config.get('potential', 'conservative') + + if 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 + elif p_nom_max_meth == 'heuristic': + p_nom_max = 0.8 * xr.DataArray(np.asarray(potmatrix.sum(axis=1)).squeeze(), [buses]) + else: + raise AssertionError('Config key `potential` should be one of "conservative" (default) or "heuristic",' + ' not "{}"'.format(p_nom_max_meth)) + + layout = xr.DataArray(np.asarray(potmatrix.sum(axis=0)).reshape(cutout.shape), + [cutout.meta.indexes[ax] for ax in ['y', 'x']]) + + ds = xr.merge([(correction_factor * profile).rename('profile'), + capacities.rename('weight'), + p_nom_max.rename('p_nom_max'), + layout.rename('potential')]) + (ds.sel(bus=(ds['profile'].mean('time') > config.get('min_p_max_pu', 0.)) & (ds['p_nom_max'] > 0.)) + .to_netcdf(snakemake.output.profile))