resolve merge conflicts

.gitignore

@@ -19,6 +19,7 @@ gurobi.log
 /data
 /data/links_p_nom.csv
 /cutouts
+/dask-worker-space
 
 doc/_build
 

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)
 =====================================

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 <https://github.com/snakemake/snakemake/issues/46>`_ at `snakemake <https://snakemake.readthedocs.io/>`_.
-    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.
 

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

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)

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 = (n.links.loc[dc_b, 'length'] * length_factor *
+            ((1. - n.links.loc[dc_b, 'underwater_fraction']) *
-                 costs.at['HVDC overhead', 'capital_cost'])
+            costs.at['HVDC overhead', 'capital_cost'] +
-    else:
+            n.links.loc[dc_b, 'underwater_fraction'] *
-        costs = (n.links.loc[dc_b, 'length'] * length_factor *
+            costs.at['HVDC submarine', 'capital_cost']) +
-                ((1. - n.links.loc[dc_b, 'underwater_fraction']) *
+            costs.at['HVDC inverter pair', 'capital_cost'])
-                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
 
 

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 shapely.geometry import Polygon
-from vresutils.graph import voronoi_partition_pts
+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

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:
+    if (start in load.index and stop in load.index):
-        load.loc[start:stop, cntry] = load.loc[start-delta:stop-delta, cntry].values
+        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)

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,6 +217,8 @@ 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')
@@ -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']),

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:

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)
 

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(',')
+    clusters = snakemake.wildcards.clusts.split(',')
-    techs = wildcards.techs.split(',')
+    techs = snakemake.wildcards.techs.split(',')
-    country = wildcards.country
+    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')

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])

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)

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