diff --git a/Snakefile b/Snakefile index 7678a401f..b1bc6e334 100644 --- a/Snakefile +++ b/Snakefile @@ -188,7 +188,7 @@ rule build_renewable_profiles: 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() + if any(key in ["max_depth", "min_depth"] for key in config["renewable"][w.technology].keys()) else []), country_shapes='resources/country_shapes.geojson', offshore_shapes='resources/offshore_shapes.geojson', @@ -226,6 +226,7 @@ rule add_electricity: geth_hydro_capacities='data/geth2015_hydro_capacities.csv', load='resources/load.csv', nuts3_shapes='resources/nuts3_shapes.geojson', + gebco='data/bundle/GEBCO_2014_2D.nc', **{f"profile_{tech}": f"resources/profile_{tech}.nc" for tech in config['renewable']} output: "networks/elec.nc" diff --git a/config.default.yaml b/config.default.yaml index d2bf61598..769e9a2e8 100755 --- a/config.default.yaml +++ b/config.default.yaml @@ -36,6 +36,7 @@ enable: build_natura_raster: false retrieve_natura_raster: true custom_busmap: false + split_offshore_regions: True #splits big offshore regions into smaller regions electricity: voltages: [220., 300., 380.] @@ -110,19 +111,37 @@ renewable: cutout: europe-2013-era5 resource: method: wind - turbine: NREL_ReferenceTurbine_5MW_offshore - capacity_per_sqkm: 2 + turbine: NREL_ReferenceTurbine_2020ATB_12MW_offshore + capacity_per_sqkm: 2.5 correction_factor: 0.8855 # proxy for wake losses # from 10.1016/j.energy.2018.08.153 # until done more rigorously in #153 corine: [44, 255] natura: true - max_depth: 50 + max_depth: 60 max_shore_distance: 30000 potential: simple # or conservative clip_p_max_pu: 1.e-2 + calculate_topology_cost: true offwind-dc: + cutout: europe-2013-era5 + resource: + method: wind + turbine: NREL_ReferenceTurbine_2020ATB_12MW_offshore + capacity_per_sqkm: 3.5 # from DEA + correction_factor: 0.8855 + # proxy for wake losses + # from 10.1016/j.energy.2018.08.153 + # until done more rigorously in #153 + corine: [44, 255] + natura: true + max_depth: 60 + min_shore_distance: 30000 + potential: simple # or conservative + clip_p_max_pu: 1.e-2 + calculate_topology_cost: true + offwind-float: cutout: europe-2013-era5 resource: method: wind @@ -135,8 +154,7 @@ renewable: # until done more rigorously in #153 corine: [44, 255] natura: true - max_depth: 50 - min_shore_distance: 30000 + min_depth: 60 potential: simple # or conservative clip_p_max_pu: 1.e-2 solar: @@ -203,7 +221,8 @@ costs: marginal_cost: # EUR/MWh solar: 0.01 onwind: 0.015 - offwind: 0.015 + offwind-ac: 0.015 + offwind-dc: 0.015 hydro: 0. H2: 0. electrolysis: 0. diff --git a/data/costs.csv b/data/costs.csv index 8953eb8a8..e84795f43 100644 --- a/data/costs.csv +++ b/data/costs.csv @@ -1,7 +1,9 @@ technology,year,parameter,value,unit,source solar-rooftop,2030,discount rate,0.04,per unit,standard for decentral onwind,2030,lifetime,30,years,DEA https://ens.dk/en/our-services/projections-and-models/technology-data -offwind,2030,lifetime,30,years,DEA https://ens.dk/en/our-services/projections-and-models/technology-data +offwind-ac,2030,lifetime,30,years,DEA https://ens.dk/en/our-services/projections-and-models/technology-data +offwind-dc,2030,lifetime,30,years,DEA https://ens.dk/en/our-services/projections-and-models/technology-data +offwind-float,2030,lifetime,20,years,C. Maienza 2020 A life cycle cost model for floating offshore wind farms solar,2030,lifetime,25,years,IEA2010 solar-rooftop,2030,lifetime,25,years,IEA2010 solar-utility,2030,lifetime,25,years,IEA2010 @@ -17,13 +19,18 @@ geothermal,2030,lifetime,40,years,IEA2010 biomass,2030,lifetime,30,years,ECF2010 in DIW DataDoc http://hdl.handle.net/10419/80348 oil,2030,lifetime,30,years,ECF2010 in DIW DataDoc http://hdl.handle.net/10419/80348 onwind,2030,investment,1040,EUR/kWel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data -offwind,2030,investment,1640,EUR/kWel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data +offwind-ac,2030,investment,1640,EUR/kWel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data offwind-ac-station,2030,investment,250,EUR/kWel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data offwind-ac-connection-submarine,2030,investment,2685,EUR/MW/km,DEA https://ens.dk/en/our-services/projections-and-models/technology-data offwind-ac-connection-underground,2030,investment,1342,EUR/MW/km,DEA https://ens.dk/en/our-services/projections-and-models/technology-data +offwind-dc,2030,investment,1640,EUR/kWel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data offwind-dc-station,2030,investment,400,EUR/kWel,Haertel 2017; assuming one onshore and one offshore node + 13% learning reduction offwind-dc-connection-submarine,2030,investment,2000,EUR/MW/km,DTU report based on Fig 34 of https://ec.europa.eu/energy/sites/ener/files/documents/2014_nsog_report.pdf offwind-dc-connection-underground,2030,investment,1000,EUR/MW/km,Haertel 2017; average + 13% learning reduction +offwind-float,2030,investment,2973,EUR/kWel,"C. Maienza 2020 A life cycle cost model for floating offshore wind farms;""Spar Buyo without cost for cabelin year not specified""" +offwind-float-station,2030,investment,400,EUR/kWel,Haertel 2017; assuming one onshore and one offshore node + 13% learning reduction +offwind-float-connection-submarine,2030,investment,2000,EUR/MW/km,DTU report based on Fig 34 of https://ec.europa.eu/energy/sites/ener/files/documents/2014_nsog_report.pdf +offwind-float-connection-underground,2030,investment,1000,EUR/MW/km,"Haertel 2017; average + 13% learning reduction""" solar,2030,investment,600,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348 biomass,2030,investment,2209,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348 geothermal,2030,investment,3392,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348 @@ -39,7 +46,9 @@ nuclear,2030,investment,6000,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80 CCGT,2030,investment,800,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348 oil,2030,investment,400,EUR/kWel,DIW DataDoc http://hdl.handle.net/10419/80348 onwind,2030,FOM,2.450549,%/year,DEA https://ens.dk/en/our-services/projections-and-models/technology-data -offwind,2030,FOM,2.304878,%/year,DEA https://ens.dk/en/our-services/projections-and-models/technology-data +offwind-ac,2030,FOM,2.304878,%/year,DEA https://ens.dk/en/our-services/projections-and-models/technology-data +offwind-dc,2030,FOM,2.304878,%/year,DEA https://ens.dk/en/our-services/projections-and-models/technology-data +offwind-float,2030,FOM,1.31,%/year,C. Maienza 2020 A life cycle cost model for floating offshore wind farms ;Assumend lifetime 20a solar,2030,FOM,4.166667,%/year,DIW DataDoc http://hdl.handle.net/10419/80348 solar-rooftop,2030,FOM,2,%/year,ETIP PV solar-utility,2030,FOM,3,%/year,ETIP PV @@ -54,7 +63,9 @@ ror,2030,FOM,2,%/year,DIW DataDoc http://hdl.handle.net/10419/80348 CCGT,2030,FOM,2.5,%/year,DIW DataDoc http://hdl.handle.net/10419/80348 OCGT,2030,FOM,3.75,%/year,DIW DataDoc http://hdl.handle.net/10419/80348 onwind,2030,VOM,2.3,EUR/MWhel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data -offwind,2030,VOM,2.7,EUR/MWhel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data +offwind-ac,2030,VOM,2.7,EUR/MWhel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data +offwind-dc,2030,VOM,2.7,EUR/MWhel,DEA https://ens.dk/en/our-services/projections-and-models/technology-data +offwind-float,2030,VOM,0,EUR/MWhel,no information; included in FOM solar,2030,VOM,0.01,EUR/MWhel,RES costs made up to fix curtailment order coal,2030,VOM,6,EUR/MWhel,DIW DataDoc http://hdl.handle.net/10419/80348 PC (Advanced/SuperC) lignite,2030,VOM,7,EUR/MWhel,DIW DataDoc http://hdl.handle.net/10419/80348 @@ -192,4 +203,4 @@ HVDC submarine,2030,lifetime,40,years,Hagspiel HVDC submarine,2030,FOM,2,%/year,Hagspiel HVDC inverter pair,2030,investment,150000,EUR/MW,Hagspiel HVDC inverter pair,2030,lifetime,40,years,Hagspiel -HVDC inverter pair,2030,FOM,2,%/year,Hagspiel +HVDC inverter pair,2030,FOM,2,%/year,Hagspiel \ No newline at end of file diff --git a/scripts/_helpers.py b/scripts/_helpers.py index 410e05afc..04afeceaa 100644 --- a/scripts/_helpers.py +++ b/scripts/_helpers.py @@ -3,6 +3,7 @@ # SPDX-License-Identifier: MIT import pandas as pd +import numpy as np from pathlib import Path @@ -199,6 +200,46 @@ def aggregate_costs(n, flatten=False, opts=None, existing_only=False): return costs +def calculate_offwind_cost(WD, MW=12, no=54, D=236, HH=138, SP=343, DT=8): + """ + Helper calculating offshore wind capex considering the average water depth of the region. + + Parameters + ---------- + WD : xarray + Average water depth of the different regions + MW : float + Power of the wind turbine in MW + no: int + Average number of wind turbines in farm + D: int + Rotor diameter of wind turbine in meters + HH: int + Hub height of wind turbine in meters + SP: int + Specific power of wind turbine in W/m2 + DT: int + Distance between the wind turbine in number of rotor diameters + + Returns + ------- + capex: xarray + Capex of the wind turbine in the different regions + """ + RA=(D/2)**2*np.pi + IA=DT*D + wind_turbine_invest=(-0.6*SP+750+(0.53*HH*RA+5500)/(1000*MW))*1.1 + wind_turbine_install= 300*MW**(-0.6) + foundation_invest=(8*np.abs(WD)+30)*(1+(0.003*(350-np.min([400,SP])))) + foundation_install=2.5*np.abs(WD)+600*MW**(-0.6) + array_cable=IA*500/MW/1000 + turbine_transport=50 + insurance=100 + finance_cost=100 + continences=50 + capex=np.sum([wind_turbine_invest,wind_turbine_install, foundation_invest, foundation_install, array_cable, turbine_transport, insurance, finance_cost, continences])*1000 # in €/MW + return capex + def progress_retrieve(url, file): import urllib from progressbar import ProgressBar @@ -240,7 +281,7 @@ def mock_snakemake(rulename, **wildcards): if os.path.exists(p): snakefile = p break - workflow = sm.Workflow(snakefile, overwrite_configfiles=[]) + workflow = sm.Workflow(snakefile, overwrite_configfiles=[], rerun_triggers=[]) workflow.include(snakefile) workflow.global_resources = {} rule = workflow.get_rule(rulename) diff --git a/scripts/add_electricity.py b/scripts/add_electricity.py index ceef23904..ab9521a0f 100755 --- a/scripts/add_electricity.py +++ b/scripts/add_electricity.py @@ -263,6 +263,7 @@ def update_transmission_costs(n, costs, length_factor=1.0): def attach_wind_and_solar(n, costs, input_profiles, technologies, line_length_factor=1): + from _helpers import calculate_offwind_cost # TODO: rename tech -> carrier, technologies -> carriers for tech in technologies: @@ -275,29 +276,45 @@ def attach_wind_and_solar(n, costs, input_profiles, technologies, line_length_fa suptech = tech.split('-', 2)[0] if suptech == 'offwind': underwater_fraction = ds['underwater_fraction'].to_pandas() - connection_cost = (line_length_factor * - ds['average_distance'].to_pandas() * - (underwater_fraction * + cable_cost = (line_length_factor * + ds['average_distance'].to_pandas() * + (underwater_fraction * costs.at[tech + '-connection-submarine', 'capital_cost'] + (1. - underwater_fraction) * costs.at[tech + '-connection-underground', 'capital_cost'])) - capital_cost = (costs.at['offwind', 'capital_cost'] + - costs.at[tech + '-station', 'capital_cost'] + - connection_cost) + grid_connection_cost=costs.at[tech + '-station', 'capital_cost'] + cable_cost + calculate_topology_cost=technologies[tech].get("calculate_topology_cost", False) + if calculate_topology_cost: + turbine_cost=calculate_offwind_cost(ds["water_depth"].to_pandas())*(calculate_annuity(costs.at[tech, 'lifetime'], costs.at[tech, "discount rate"]) + costs.at[tech, "FOM"]/100.) * Nyears + else: + turbine_cost=costs.at[tech, 'capital_cost'] + capital_cost = (turbine_cost + grid_connection_cost) + logger.info("Added connection cost of {:0.0f}-{:0.0f} Eur/MW/a to {}" - .format(connection_cost.min(), connection_cost.max(), tech)) + .format(cable_cost.min(), cable_cost.max(), tech)) + n.madd("Generator", ds.indexes['bus'], ' ' + tech, + bus=ds.indexes['bus'].str.split("_").str[0], + carrier=tech, + p_nom_extendable=True, + p_nom_max=ds['p_nom_max'].to_pandas(), + weight=ds['weight'].to_pandas(), + marginal_cost=costs.at[tech, 'marginal_cost'], + capital_cost=capital_cost, + grid_connection_cost=grid_connection_cost, + turbine_cost=turbine_cost, + efficiency=costs.at[tech, 'efficiency'], + p_max_pu=ds['profile'].transpose('time', 'bus').to_pandas()) else: capital_cost = costs.at[tech, 'capital_cost'] - - n.madd("Generator", ds.indexes['bus'], ' ' + tech, - bus=ds.indexes['bus'], + n.madd("Generator", ds.indexes['bus'], ' ' + tech, + bus=ds.indexes['bus'].str.split("_").str[0], carrier=tech, p_nom_extendable=True, p_nom_max=ds['p_nom_max'].to_pandas(), weight=ds['weight'].to_pandas(), - marginal_cost=costs.at[suptech, 'marginal_cost'], + marginal_cost=costs.at[tech, 'marginal_cost'], capital_cost=capital_cost, - efficiency=costs.at[suptech, 'efficiency'], + efficiency=costs.at[tech, 'efficiency'], p_max_pu=ds['profile'].transpose('time', 'bus').to_pandas()) diff --git a/scripts/build_bus_regions.py b/scripts/build_bus_regions.py index d91d0575b..d5491d030 100644 --- a/scripts/build_bus_regions.py +++ b/scripts/build_bus_regions.py @@ -47,13 +47,101 @@ import pypsa import os import pandas as pd +import numpy as np import geopandas as gpd +from pyproj import Geod +from shapely import affinity +from shapely.geometry import Point +from sklearn.cluster import KMeans from vresutils.graph import voronoi_partition_pts logger = logging.getLogger(__name__) +def calculate_area(shape, ellipsoid="WGS84"): + geod = Geod(ellps=ellipsoid) + return abs(geod.geometry_area_perimeter(shape)[0])/1e6 + +def cluster_points(n_clusters, point_list): + ''' + Clusters the inner points of a region into n_clusters. + + Parameters + ---------- + n_clusters : + Number of clusters + point_list : + List of inner points. + + Returns + ------- + Returns list of cluster centers. + ''' + kmeans = KMeans(n_clusters=n_clusters, random_state=0).fit(point_list) + return kmeans.cluster_centers_ + + +def fill_shape_with_points(shape, oversize_factor, num=10): + ''' + Fills the shape of the offshore region with points. This is needed for splitting the regions into smaller regions. + + Parameters + ---------- + shape : + Shape of the region. + oversize_factor : int + Factor by which the original region is oversized. + num : int, optional + Number of points added in the x and y direction. + + Returns + ------- + inner_points : + Returns a list of points lying inside the shape. + ''' + + inner_points=list() + x_min, y_min, x_max, y_max= shape.bounds + iteration=0 + while True: + for x in np.linspace(x_min, x_max, num=num): + for y in np.linspace(y_min, y_max, num=num): + if Point(x,y).within(shape): + inner_points.append((x,y)) + if len(inner_points) > oversize_factor: + break + else: + #perturb bounds that not the same points are added again + num+=1 + x_min += abs(x_max-x_min)*0.01 + x_max -= abs(x_max-x_min)*0.01 + y_min += abs(y_max-y_min)*0.01 + y_max -= abs(y_max-y_min)*0.01 + return inner_points + +def build_voronoi_cells(shape, points): + ''' + Builds Voronoi cells from given points in the given shape. + + Parameters + ---------- + shape : + Shape where to build the cells. + points : + List of points. + + Returns + ------- + split region + Geopandas DataFrame containing the split regions. + ''' + split_region=gpd.GeoDataFrame({ + 'x': points[:,0], + 'y': points[:,0], + 'geometry': voronoi_partition_pts(points, shape), + }) + return split_region def save_to_geojson(s, fn): if os.path.exists(fn): os.unlink(fn) @@ -101,6 +189,22 @@ def save_to_geojson(s, fn): 'country': country }) offshore_regions_c = offshore_regions_c.loc[offshore_regions_c.area > 1e-2] + split_offshore_regions=snakemake.config["enable"].get('split_offshore_regions', False) + offshore_regions_c.drop_duplicates(subset="geometry", inplace=True) # some regions are duplicated + if not offshore_regions_c.empty and split_offshore_regions: + threshold=15000 #km2 threshold at which regions are splitted + region_oversize=offshore_regions_c.geometry.map(lambda x: calculate_area(x)/threshold) + for bus, region in offshore_regions_c[region_oversize>1].iterrows(): + shape=region.geometry + oversize_factor=region_oversize.loc[bus] + inner_points=fill_shape_with_points(shape, oversize_factor) + cluster_centers=cluster_points(int(np.ceil(oversize_factor)), inner_points) + inner_regions=build_voronoi_cells(shape, cluster_centers) + inner_regions.set_index(pd.Index([f"{bus}_{i}" for i in inner_regions.index], name="Bus"), inplace=True) + inner_regions['name']=inner_regions.index + inner_regions['country']=country + offshore_regions_c=pd.concat([offshore_regions_c.drop(bus),inner_regions]) + offshore_regions_c["area"]=offshore_regions_c.geometry.apply(lambda x: calculate_area(x)) offshore_regions.append(offshore_regions_c) save_to_geojson(pd.concat(onshore_regions, ignore_index=True), snakemake.output.regions_onshore) diff --git a/scripts/build_renewable_profiles.py b/scripts/build_renewable_profiles.py index a2b2eda6d..af7916a2b 100644 --- a/scripts/build_renewable_profiles.py +++ b/scripts/build_renewable_profiles.py @@ -31,6 +31,7 @@ distance: natura: max_depth: + min_depth: max_shore_distance: min_shore_distance: capacity_per_sqkm: @@ -198,7 +199,7 @@ if __name__ == '__main__': if 'snakemake' not in globals(): from _helpers import mock_snakemake - snakemake = mock_snakemake('build_renewable_profiles', technology='solar') + snakemake = mock_snakemake('build_renewable_profiles', technology='offwind-ac') configure_logging(snakemake) pgb.streams.wrap_stderr() @@ -240,7 +241,11 @@ # use named function np.greater with partially frozen argument instead # and exclude areas where: -max_depth > grid cell depth func = functools.partial(np.greater,-config['max_depth']) - excluder.add_raster(snakemake.input.gebco, codes=func, crs=4236, nodata=-1000) + excluder.add_raster(snakemake.input.gebco, codes=func, crs=4236, nodata=-1000, allow_no_overlap=True) + + if "min_depth" in config: + func = functools.partial(np.greater,-config['min_depth']) + excluder.add_raster(snakemake.input.gebco, codes=func, crs=4236, nodata=-1000, invert=True, allow_no_overlap=True) if 'min_shore_distance' in config: buffer = config['min_shore_distance'] @@ -312,8 +317,15 @@ potential.rename('potential'), average_distance.rename('average_distance')]) - - if snakemake.wildcards.technology.startswith("offwind"): + tech=snakemake.wildcards.technology + if tech.startswith("offwind"): + if tech.endswith!= "float": + with xr.open_dataset(snakemake.input.gebco) as gebco: + from rasterio.warp import Resampling + gebco=gebco.rename({"lon":"x", "lat":"y"}) + water_depth=atlite.gis.regrid(gebco.elevation,cutout.data.x, cutout.data.y,resampling=Resampling.average) + water_depth=(water_depth*availability.where(availability!=0)).mean(["x", "y"]) + ds['water_depth'] = xr.DataArray(water_depth, [buses]).fillna(0).clip(max=0) logger.info('Calculate underwater fraction of connections.') offshore_shape = gpd.read_file(snakemake.input['offshore_shapes']).unary_union underwater_fraction = [] @@ -322,7 +334,6 @@ 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