Merge branch 'master' into voronoi

doc/release_notes.rst

@@ -62,12 +62,19 @@ Upcoming Release
 
 * New network topology extracted from the ENTSO-E interactive map.
 
+* 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.
 
 * Switch from Germany to Belgium for continuous integration and tutorial to save resources.
 
 * 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,7 +45,7 @@ dependencies:
   # GIS dependencies:
   - cartopy
   - descartes
-  - rasterio
+  - rasterio<=1.2.9  # 1.2.10 creates error https://github.com/PyPSA/atlite/issues/238
 
   # PyPSA-Eur-Sec Dependencies
   - geopy

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_renewable_profiles.py

@@ -240,7 +240,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']

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/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