diff --git a/.gitignore b/.gitignore index b4734ab2..559dde47 100644 --- a/.gitignore +++ b/.gitignore @@ -19,6 +19,7 @@ gurobi.log /data /data/links_p_nom.csv /cutouts +/dask-worker-space doc/_build diff --git a/doc/release_notes.rst b/doc/release_notes.rst index b514528d..2941678c 100644 --- a/doc/release_notes.rst +++ b/doc/release_notes.rst @@ -67,6 +67,7 @@ Upcoming Release * Cache data and cutouts folders. This cache will be updated weekly. * Add rule to automatically retrieve Natura2000 natural protection areas. Switch of file format to GPKG. +* The unused argument ``simple_hvdc_costs`` in :mod:`add_electricity` was removed. * Iterative solving with impedance updates is skipped if there are no expandable lines. @@ -74,6 +75,10 @@ Upcoming Release * Use updated SARAH-2 and ERA5 cutouts with slightly wider scope to east and additional variables. +* Fix crs bug. Change crs 4236 to 4326. + +* Update rasterio version to correctly calculate exclusion raster + PyPSA-Eur 0.4.0 (22th September 2021) ===================================== diff --git a/doc/tutorial.rst b/doc/tutorial.rst index 512683a5..1c247f8e 100644 --- a/doc/tutorial.rst +++ b/doc/tutorial.rst @@ -43,7 +43,7 @@ For more information on the data dependencies of PyPSA-Eur, continue reading :re How to customise PyPSA-Eur? =========================== -The model can be adapted to only include selected countries (e.g. Germany) instead of all European countries to limit the spatial scope. +The model can be adapted to only include selected countries (e.g. Belgium) instead of all European countries to limit the spatial scope. .. literalinclude:: ../config.tutorial.yaml :language: yaml @@ -53,41 +53,43 @@ Likewise, the example's temporal scope can be restricted (e.g. to a single month .. literalinclude:: ../config.tutorial.yaml :language: yaml - :lines: 24-27 + :start-at: snapshots: + :end-before: enable: It is also possible to allow less or more carbon-dioxide emissions. Here, we limit the emissions of Germany 100 Megatonnes per year. .. literalinclude:: ../config.tutorial.yaml :language: yaml - :lines: 38,40 + :lines: 40,42 PyPSA-Eur also includes a database of existing conventional powerplants. We can select which types of powerplants we like to be included with fixed capacities: .. literalinclude:: ../config.tutorial.yaml :language: yaml - :lines: 38,54 + :lines: 40,56 To accurately model the temporal and spatial availability of renewables such as wind and solar energy, we rely on historical weather data. It is advisable to adapt the required range of coordinates to the selection of countries. .. literalinclude:: ../config.tutorial.yaml :language: yaml - :lines: 56-63 + :start-at: atlite: + :end-before: renewable: We can also decide which weather data source should be used to calculate potentials and capacity factor time-series for each carrier. For example, we may want to use the ERA-5 dataset for solar and not the default SARAH-2 dataset. .. literalinclude:: ../config.tutorial.yaml :language: yaml - :lines: 65,108-109 + :lines: 67,110,111 Finally, it is possible to pick a solver. For instance, this tutorial uses the open-source solvers CBC and Ipopt and does not rely on the commercial solvers Gurobi or CPLEX (for which free academic licenses are available). .. literalinclude:: ../config.tutorial.yaml :language: yaml - :lines: 171,181-182 + :lines: 173,183,184 .. note:: @@ -126,11 +128,6 @@ orders ``snakemake`` to run the script ``solve_network`` that produces the solve .. until https://github.com/snakemake/snakemake/issues/46 closed -.. warning:: - On Windows the previous command may currently cause a ``MissingRuleException`` due to problems with output files in subfolders. - This is an `open issue `_ at `snakemake `_. - Windows users should add the option ``--keep-target-files`` to the command or instead run ``snakemake -j 1 solve_all_networks``. - This triggers a workflow of multiple preceding jobs that depend on each rule's inputs and outputs: .. graphviz:: @@ -271,7 +268,8 @@ the wildcards given in ``scenario`` in the configuration file ``config.yaml`` ar .. literalinclude:: ../config.tutorial.yaml :language: yaml - :lines: 14-18 + :start-at: scenario: + :end-before: countries: In this example we would not only solve a 6-node model of Germany but also a 2-node model. diff --git a/envs/environment.yaml b/envs/environment.yaml index 31d037b1..fd6dff49 100644 --- a/envs/environment.yaml +++ b/envs/environment.yaml @@ -11,7 +11,7 @@ dependencies: - pip - pypsa>=0.18.1 - - atlite>=0.2.5 + - atlite>=0.2.6 - dask # Dependencies of the workflow itself @@ -24,7 +24,7 @@ dependencies: - yaml - pytables - lxml - - powerplantmatching>=0.4.8 + - powerplantmatching>=0.5.3 - numpy - pandas - geopandas @@ -37,7 +37,7 @@ dependencies: - pyomo - matplotlib - proj - - fiona <= 1.18.20 # Till issue https://github.com/Toblerity/Fiona/issues/1085 is not solved + - fiona<=1.18.20 # Till issue https://github.com/Toblerity/Fiona/issues/1085 is not solved # Keep in conda environment when calling ipython - ipython @@ -45,8 +45,8 @@ dependencies: # GIS dependencies: - cartopy - descartes - - rasterio - fiona # explicit for Windows + - rasterio<=1.2.9 # 1.2.10 creates error https://github.com/PyPSA/atlite/issues/238 # PyPSA-Eur-Sec Dependencies - geopy diff --git a/scripts/_helpers.py b/scripts/_helpers.py index 410e05af..6e47c053 100644 --- a/scripts/_helpers.py +++ b/scripts/_helpers.py @@ -231,6 +231,7 @@ def mock_snakemake(rulename, **wildcards): import os from pypsa.descriptors import Dict from snakemake.script import Snakemake + from packaging.version import Version, parse script_dir = Path(__file__).parent.resolve() assert Path.cwd().resolve() == script_dir, \ @@ -240,7 +241,8 @@ def mock_snakemake(rulename, **wildcards): if os.path.exists(p): snakefile = p break - workflow = sm.Workflow(snakefile, overwrite_configfiles=[]) + kwargs = dict(rerun_triggers=[]) if parse(sm.__version__) > Version("7.7.0") else {} + workflow = sm.Workflow(snakefile, overwrite_configfiles=[], **kwargs) workflow.include(snakefile) workflow.global_resources = {} rule = workflow.get_rule(rulename) diff --git a/scripts/add_electricity.py b/scripts/add_electricity.py index 7dffe60f..ceef2390 100755 --- a/scripts/add_electricity.py +++ b/scripts/add_electricity.py @@ -94,7 +94,6 @@ import geopandas as gpd import powerplantmatching as pm from powerplantmatching.export import map_country_bus -from vresutils.costdata import annuity from vresutils import transfer as vtransfer idx = pd.IndexSlice @@ -105,6 +104,18 @@ logger = logging.getLogger(__name__) def normed(s): return s/s.sum() +def calculate_annuity(n, r): + """Calculate the annuity factor for an asset with lifetime n years and + discount rate of r, e.g. annuity(20, 0.05) * 20 = 1.6""" + + if isinstance(r, pd.Series): + return pd.Series(1/n, index=r.index).where(r == 0, r/(1. - 1./(1.+r)**n)) + elif r > 0: + return r / (1. - 1./(1.+r)**n) + else: + return 1 / n + + def _add_missing_carriers_from_costs(n, costs, carriers): missing_carriers = pd.Index(carriers).difference(n.carriers.index) if missing_carriers.empty: return @@ -138,7 +149,7 @@ def load_costs(tech_costs, config, elec_config, Nyears=1.): "investment" : 0, "lifetime" : 25}) - costs["capital_cost"] = ((annuity(costs["lifetime"], costs["discount rate"]) + + costs["capital_cost"] = ((calculate_annuity(costs["lifetime"], costs["discount rate"]) + costs["FOM"]/100.) * costs["investment"] * Nyears) @@ -227,7 +238,7 @@ def attach_load(n, regions, load, nuts3_shapes, countries, scaling=1.): n.madd("Load", substation_lv_i, bus=substation_lv_i, p_set=load) -def update_transmission_costs(n, costs, length_factor=1.0, simple_hvdc_costs=False): +def update_transmission_costs(n, costs, length_factor=1.0): # TODO: line length factor of lines is applied to lines and links. # Separate the function to distinguish. @@ -242,16 +253,12 @@ def update_transmission_costs(n, costs, length_factor=1.0, simple_hvdc_costs=Fal # may be missing. Therefore we have to return here. if n.links.loc[dc_b].empty: return - if simple_hvdc_costs: - costs = (n.links.loc[dc_b, 'length'] * length_factor * - costs.at['HVDC overhead', 'capital_cost']) - else: - costs = (n.links.loc[dc_b, 'length'] * length_factor * - ((1. - n.links.loc[dc_b, 'underwater_fraction']) * - costs.at['HVDC overhead', 'capital_cost'] + - n.links.loc[dc_b, 'underwater_fraction'] * - costs.at['HVDC submarine', 'capital_cost']) + - costs.at['HVDC inverter pair', 'capital_cost']) + costs = (n.links.loc[dc_b, 'length'] * length_factor * + ((1. - n.links.loc[dc_b, 'underwater_fraction']) * + costs.at['HVDC overhead', 'capital_cost'] + + n.links.loc[dc_b, 'underwater_fraction'] * + costs.at['HVDC submarine', 'capital_cost']) + + costs.at['HVDC inverter pair', 'capital_cost']) n.links.loc[dc_b, 'capital_cost'] = costs diff --git a/scripts/build_bus_regions.py b/scripts/build_bus_regions.py index d91d0575..8003d370 100644 --- a/scripts/build_bus_regions.py +++ b/scripts/build_bus_regions.py @@ -47,9 +47,10 @@ from _helpers import configure_logging import pypsa import os import pandas as pd +import numpy as np import geopandas as gpd - -from vresutils.graph import voronoi_partition_pts +from shapely.geometry import Polygon +from scipy.spatial import Voronoi logger = logging.getLogger(__name__) @@ -61,6 +62,53 @@ def save_to_geojson(s, fn): s.to_file(fn, driver='GeoJSON', schema=schema) +def voronoi_partition_pts(points, outline): + """ + Compute the polygons of a voronoi partition of `points` within the + polygon `outline`. Taken from + https://github.com/FRESNA/vresutils/blob/master/vresutils/graph.py + Attributes + ---------- + points : Nx2 - ndarray[dtype=float] + outline : Polygon + Returns + ------- + polygons : N - ndarray[dtype=Polygon|MultiPolygon] + """ + + points = np.asarray(points) + + if len(points) == 1: + polygons = [outline] + else: + xmin, ymin = np.amin(points, axis=0) + xmax, ymax = np.amax(points, axis=0) + xspan = xmax - xmin + yspan = ymax - ymin + + # to avoid any network positions outside all Voronoi cells, append + # the corners of a rectangle framing these points + vor = Voronoi(np.vstack((points, + [[xmin-3.*xspan, ymin-3.*yspan], + [xmin-3.*xspan, ymax+3.*yspan], + [xmax+3.*xspan, ymin-3.*yspan], + [xmax+3.*xspan, ymax+3.*yspan]]))) + + polygons = [] + for i in range(len(points)): + poly = Polygon(vor.vertices[vor.regions[vor.point_region[i]]]) + + if not poly.is_valid: + poly = poly.buffer(0) + + poly = poly.intersection(outline) + + polygons.append(poly) + + + return np.array(polygons, dtype=object) + + if __name__ == "__main__": if 'snakemake' not in globals(): from _helpers import mock_snakemake diff --git a/scripts/build_load_data.py b/scripts/build_load_data.py index 10921782..55270e49 100755 --- a/scripts/build_load_data.py +++ b/scripts/build_load_data.py @@ -116,14 +116,19 @@ def nan_statistics(df): keys=['total', 'consecutive', 'max_total_per_month'], axis=1) -def copy_timeslice(load, cntry, start, stop, delta): +def copy_timeslice(load, cntry, start, stop, delta, fn_load=None): start = pd.Timestamp(start) stop = pd.Timestamp(stop) - if start-delta in load.index and stop in load.index and cntry in load: - load.loc[start:stop, cntry] = load.loc[start-delta:stop-delta, cntry].values + if (start in load.index and stop in load.index): + if start-delta in load.index and stop-delta in load.index and cntry in load: + load.loc[start:stop, cntry] = load.loc[start-delta:stop-delta, cntry].values + elif fn_load is not None: + duration = pd.date_range(freq='h', start=start-delta, end=stop-delta) + load_raw = load_timeseries(fn_load, duration, [cntry], powerstatistics) + load.loc[start:stop, cntry] = load_raw.loc[start-delta:stop-delta, cntry].values -def manual_adjustment(load, powerstatistics): +def manual_adjustment(load, fn_load, powerstatistics): """ Adjust gaps manual for load data from OPSD time-series package. @@ -150,6 +155,8 @@ def manual_adjustment(load, powerstatistics): powerstatistics: bool Whether argument load comprises the electricity consumption data of the ENTSOE power statistics or of the ENTSOE transparency map + load_fn: str + File name or url location (file format .csv) Returns ------- @@ -175,7 +182,11 @@ def manual_adjustment(load, powerstatistics): copy_timeslice(load, 'CH', '2010-11-04 04:00', '2010-11-04 22:00', Delta(days=1)) copy_timeslice(load, 'NO', '2010-12-09 11:00', '2010-12-09 18:00', Delta(days=1)) # whole january missing - copy_timeslice(load, 'GB', '2009-12-31 23:00', '2010-01-31 23:00', Delta(days=-364)) + copy_timeslice(load, 'GB', '2010-01-01 00:00', '2010-01-31 23:00', Delta(days=-365), fn_load) + # 1.1. at midnight gets special treatment + copy_timeslice(load, 'IE', '2016-01-01 00:00', '2016-01-01 01:00', Delta(days=-366), fn_load) + copy_timeslice(load, 'PT', '2016-01-01 00:00', '2016-01-01 01:00', Delta(days=-366), fn_load) + copy_timeslice(load, 'GB', '2016-01-01 00:00', '2016-01-01 01:00', Delta(days=-366), fn_load) else: if 'ME' in load: @@ -206,7 +217,7 @@ if __name__ == "__main__": load = load_timeseries(snakemake.input[0], years, countries, powerstatistics) if snakemake.config['load']['manual_adjustments']: - load = manual_adjustment(load, powerstatistics) + load = manual_adjustment(load, snakemake.input[0], powerstatistics) logger.info(f"Linearly interpolate gaps of size {interpolate_limit} and less.") load = load.interpolate(method='linear', limit=interpolate_limit) diff --git a/scripts/build_renewable_profiles.py b/scripts/build_renewable_profiles.py index b6ac6c65..a49df1f5 100644 --- a/scripts/build_renewable_profiles.py +++ b/scripts/build_renewable_profiles.py @@ -189,6 +189,7 @@ import logging from pypsa.geo import haversine from shapely.geometry import LineString import time +from dask.distributed import Client, LocalCluster from _helpers import configure_logging @@ -203,7 +204,7 @@ if __name__ == '__main__': pgb.streams.wrap_stderr() nprocesses = int(snakemake.threads) - noprogress = not snakemake.config['atlite'].get('show_progress', True) + noprogress = not snakemake.config['atlite'].get('show_progress', False) config = snakemake.config['renewable'][snakemake.wildcards.technology] resource = config['resource'] # pv panel config / wind turbine config correction_factor = config.get('correction_factor', 1.) @@ -216,7 +217,9 @@ if __name__ == '__main__': if correction_factor != 1.: logger.info(f'correction_factor is set as {correction_factor}') - + cluster = LocalCluster(n_workers=nprocesses, threads_per_worker=1) + client = Client(cluster, asynchronous=True) + cutout = atlite.Cutout(snakemake.input['cutout']) regions = gpd.read_file(snakemake.input.regions).set_index('name').rename_axis('bus') buses = regions.index @@ -242,7 +245,7 @@ if __name__ == '__main__': # 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=4326, nodata=-1000) if 'min_shore_distance' in config: buffer = config['min_shore_distance'] @@ -268,7 +271,7 @@ if __name__ == '__main__': potential = capacity_per_sqkm * availability.sum('bus') * area func = getattr(cutout, resource.pop('method')) - resource['dask_kwargs'] = {'num_workers': nprocesses} + resource['dask_kwargs'] = {"scheduler": client} 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']), diff --git a/scripts/cluster_network.py b/scripts/cluster_network.py index 642db4da..1d5608e2 100644 --- a/scripts/cluster_network.py +++ b/scripts/cluster_network.py @@ -281,7 +281,14 @@ def clustering_for_n_clusters(n, n_clusters, custom_busmap=False, aggregate_carr aggregate_generators_carriers=aggregate_carriers, aggregate_one_ports=["Load", "StorageUnit"], line_length_factor=line_length_factor, - generator_strategies={'p_nom_max': p_nom_max_strategy, 'p_nom_min': pd.Series.sum}, + generator_strategies={'p_nom_max': p_nom_max_strategy, + 'p_nom_min': pd.Series.sum, + 'p_min_pu': pd.Series.mean, + 'marginal_cost': pd.Series.mean, + 'committable': np.any, + 'ramp_limit_up': pd.Series.max, + 'ramp_limit_down': pd.Series.max, + }, scale_link_capital_costs=False) if not n.links.empty: diff --git a/scripts/make_summary.py b/scripts/make_summary.py index 854e9463..c070d33f 100644 --- a/scripts/make_summary.py +++ b/scripts/make_summary.py @@ -171,6 +171,9 @@ def calculate_capacity(n,label,capacity): if 'p_nom_opt' in c.df.columns: c_capacities = abs(c.df.p_nom_opt.multiply(c.df.sign)).groupby(c.df.carrier).sum() capacity = include_in_summary(capacity, [c.list_name], label, c_capacities) + elif 'e_nom_opt' in c.df.columns: + c_capacities = abs(c.df.e_nom_opt.multiply(c.df.sign)).groupby(c.df.carrier).sum() + capacity = include_in_summary(capacity, [c.list_name], label, c_capacities) for c in n.iterate_components(n.passive_branch_components): c_capacities = c.df['s_nom_opt'].groupby(c.df.carrier).sum() @@ -185,11 +188,11 @@ def calculate_capacity(n,label,capacity): def calculate_supply(n, label, supply): """calculate the max dispatch of each component at the buses where the loads are attached""" - load_types = n.loads.carrier.value_counts().index + load_types = n.buses.carrier.unique() for i in load_types: - buses = n.loads.bus[n.loads.carrier == i].values + buses = n.buses.query("carrier == @i").index bus_map = pd.Series(False,index=n.buses.index) @@ -232,11 +235,11 @@ def calculate_supply(n, label, supply): def calculate_supply_energy(n, label, supply_energy): """calculate the total dispatch of each component at the buses where the loads are attached""" - load_types = n.loads.carrier.value_counts().index + load_types = n.buses.carrier.unique() for i in load_types: - buses = n.loads.bus[n.loads.carrier == i].values + buses = n.buses.query("carrier == @i").index bus_map = pd.Series(False,index=n.buses.index) @@ -404,7 +407,7 @@ def make_summaries(networks_dict, paths, config, country='all'): Nyears = n.snapshot_weightings.objective.sum() / 8760. costs = load_costs(paths[0], config['costs'], config['electricity'], Nyears) - update_transmission_costs(n, costs, simple_hvdc_costs=False) + update_transmission_costs(n, costs) assign_carriers(n) diff --git a/scripts/plot_p_nom_max.py b/scripts/plot_p_nom_max.py index ea66d612..e79ad274 100644 --- a/scripts/plot_p_nom_max.py +++ b/scripts/plot_p_nom_max.py @@ -19,7 +19,7 @@ Description """ import logging -from _helpers import configure_logging, retrieve_snakemake_keys +from _helpers import configure_logging import pypsa import pandas as pd @@ -53,13 +53,11 @@ if __name__ == "__main__": clusts= '5,full', country= 'all') configure_logging(snakemake) - paths, config, wildcards, logs, out = retrieve_snakemake_keys(snakemake) - plot_kwds = dict(drawstyle="steps-post") - clusters = wildcards.clusts.split(',') - techs = wildcards.techs.split(',') - country = wildcards.country + clusters = snakemake.wildcards.clusts.split(',') + techs = snakemake.wildcards.techs.split(',') + country = snakemake.wildcards.country if country == 'all': country = None else: @@ -68,7 +66,7 @@ if __name__ == "__main__": fig, axes = plt.subplots(1, len(techs)) for j, cluster in enumerate(clusters): - net = pypsa.Network(paths[j]) + net = pypsa.Network(snakemake.input[j]) for i, tech in enumerate(techs): cum_p_nom_max(net, tech, country).plot(x="p_max_pu", y="cum_p_nom_max", @@ -81,4 +79,4 @@ if __name__ == "__main__": plt.legend(title="Cluster level") - fig.savefig(out[0], transparent=True, bbox_inches='tight') + fig.savefig(snakemake.output[0], transparent=True, bbox_inches='tight') diff --git a/scripts/plot_summary.py b/scripts/plot_summary.py index 48f064b0..bc2bd30c 100644 --- a/scripts/plot_summary.py +++ b/scripts/plot_summary.py @@ -21,7 +21,7 @@ Description import os import logging -from _helpers import configure_logging, retrieve_snakemake_keys +from _helpers import configure_logging import pandas as pd import matplotlib.pyplot as plt @@ -170,12 +170,12 @@ if __name__ == "__main__": attr='', ext='png', country='all') configure_logging(snakemake) - paths, config, wildcards, logs, out = retrieve_snakemake_keys(snakemake) + config = snakemake.config - summary = wildcards.summary + summary = snakemake.wildcards.summary try: func = globals()[f"plot_{summary}"] except KeyError: raise RuntimeError(f"plotting function for {summary} has not been defined") - func(os.path.join(paths[0], f"{summary}.csv"), config, out[0]) + func(os.path.join(snakemake.input[0], f"{summary}.csv"), config, snakemake.output[0]) diff --git a/scripts/prepare_links_p_nom.py b/scripts/prepare_links_p_nom.py index 6bd4bca4..b83089d6 100644 --- a/scripts/prepare_links_p_nom.py +++ b/scripts/prepare_links_p_nom.py @@ -37,7 +37,7 @@ Description """ import logging -from _helpers import configure_logging, retrieve_snakemake_keys +from _helpers import configure_logging import pandas as pd @@ -63,8 +63,6 @@ if __name__ == "__main__": snakemake = mock_snakemake('prepare_links_p_nom', simpl='', network='elec') configure_logging(snakemake) - paths, config, wildcards, logs, out = retrieve_snakemake_keys(snakemake) - links_p_nom = pd.read_html('https://en.wikipedia.org/wiki/List_of_HVDC_projects', header=0, match="SwePol")[0] mw = "Power (MW)" @@ -76,4 +74,4 @@ if __name__ == "__main__": links_p_nom['x1'], links_p_nom['y1'] = extract_coordinates(links_p_nom['Converterstation 1']) links_p_nom['x2'], links_p_nom['y2'] = extract_coordinates(links_p_nom['Converterstation 2']) - links_p_nom.dropna(subset=['x1', 'y1', 'x2', 'y2']).to_csv(out[0], index=False) + links_p_nom.dropna(subset=['x1', 'y1', 'x2', 'y2']).to_csv(snakemake.output[0], index=False) diff --git a/scripts/prepare_network.py b/scripts/prepare_network.py index f984ace6..206e220b 100755 --- a/scripts/prepare_network.py +++ b/scripts/prepare_network.py @@ -104,7 +104,7 @@ def set_transmission_limit(n, ll_type, factor, costs, Nyears=1): ref = (lines_s_nom @ n.lines[col] + n.links.loc[links_dc_b, "p_nom"] @ n.links.loc[links_dc_b, col]) - update_transmission_costs(n, costs, simple_hvdc_costs=False) + update_transmission_costs(n, costs) if factor == 'opt' or float(factor) > 1.0: n.lines['s_nom_min'] = lines_s_nom