fix Nyears scaling

scripts/build_transport_demand.py

@@ -83,7 +83,7 @@ def build_transport_demand(traffic_fn, airtemp_fn, nodes, nodal_transport_data):
     )
 
     transport = (
-        (transport_shape.multiply(energy_totals_transport) * 1e6 * Nyears)
+        (transport_shape.multiply(energy_totals_transport) * 1e6 * nyears)
         .divide(efficiency_gain * ice_correction)
         .multiply(1 + dd_EV)
     )
@@ -181,7 +181,7 @@ if __name__ == "__main__":
 
     snapshots = pd.date_range(freq="h", **snakemake.config["snapshots"], tz="UTC")
 
-    Nyears = 1
+    nyears = len(snapshots) / 8760
 
     nodal_transport_data = build_nodal_transport_data(
         snakemake.input.transport_data, pop_layout

scripts/make_summary.py

@@ -699,7 +699,7 @@ if __name__ == "__main__":
         for planning_horizon in snakemake.config["scenario"]["planning_horizons"]
     }
 
-    Nyears = 1
+    Nyears = len(pd.date_range(freq="h", **snakemake.config["snapshots"])) / 8760
 
     costs_db = prepare_costs(
         snakemake.input.costs,

scripts/plot_network.py

@@ -23,7 +23,7 @@ from make_summary import assign_carriers
 from plot_summary import preferred_order, rename_techs
 from pypsa.plot import add_legend_circles, add_legend_lines, add_legend_patches
 
-plt.style.use(["ggplot", "matplotlibrc"])
+# plt.style.use(["ggplot", "matplotlibrc"])
 
 
 def rename_techs_tyndp(tech):
@@ -394,8 +394,7 @@ def plot_h2_map(network, regions):
         link_widths=link_widths_retro,
         branch_components=["Link"],
         ax=ax,
-        color_geomap=False,
+        **map_opts,
-        boundaries=map_opts["boundaries"],
     )
 
     regions.plot(
@@ -922,11 +921,11 @@ if __name__ == "__main__":
         snakemake = mock_snakemake(
             "plot_network",
             simpl="",
-            clusters="181",
-            ll="vopt",
             opts="",
-            sector_opts="Co2L0-730H-T-H-B-I-A-solar+p3-linemaxext10",
+            clusters="5",
-            planning_horizons="2050",
+            ll="v1.5",
+            sector_opts="CO2L0-1H-T-H-B-I-A-solar+p3-dist1",
+            planning_horizons="2030",
         )
 
     logging.basicConfig(level=snakemake.config["logging"]["level"])
@@ -938,6 +937,9 @@ if __name__ == "__main__":
 
     map_opts = snakemake.config["plotting"]["map"]
 
+    if map_opts["boundaries"] is None:
+        map_opts["boundaries"] = regions.total_bounds[[0, 2, 1, 3]] + [-1, 1, -1, 1]
+
     plot_map(
         n,
         components=["generators", "links", "stores", "storage_units"],

scripts/prepare_sector_network.py

@@ -676,7 +676,7 @@ def add_dac(n, costs):
     )
 
 
-def add_co2limit(n, Nyears=1.0, limit=0.0):
+def add_co2limit(n, nyears=1.0, limit=0.0):
     logger.info(f"Adding CO2 budget limit as per unit of 1990 levels of {limit}")
 
     countries = snakemake.config["countries"]
@@ -688,7 +688,7 @@ def add_co2limit(n, Nyears=1.0, limit=0.0):
 
     co2_limit = co2_totals.loc[countries, sectors].sum().sum()
 
-    co2_limit *= limit * Nyears
+    co2_limit *= limit * nyears
 
     n.add(
         "GlobalConstraint",
@@ -732,7 +732,7 @@ def cycling_shift(df, steps=1):
     return df
 
 
-def prepare_costs(cost_file, config, Nyears):
+def prepare_costs(cost_file, config, nyears):
     # set all asset costs and other parameters
     costs = pd.read_csv(cost_file, index_col=[0, 1]).sort_index()
 
@@ -750,7 +750,7 @@ def prepare_costs(cost_file, config, Nyears):
         return annuity(v["lifetime"], v["discount rate"]) + v["FOM"] / 100
 
     costs["fixed"] = [
-        annuity_factor(v) * v["investment"] * Nyears for i, v in costs.iterrows()
+        annuity_factor(v) * v["investment"] * nyears for i, v in costs.iterrows()
     ]
 
     return costs
@@ -1557,8 +1557,7 @@ def add_land_transport(n, costs):
         co2 = (
             ice_share
             / ice_efficiency
-            * transport[nodes].sum().sum()
+            * transport[nodes].sum(axis=1)
-            / 8760
             * costs.at["oil", "CO2 intensity"]
         )
 
@@ -2333,11 +2332,13 @@ def add_industry(n, costs):
     logger.info("Add industrial demand")
 
     nodes = pop_layout.index
+    nsnapshots = n.snapshot_weightings.generators.sum()
+    nyears = nsnapshots / 8760
 
     # 1e6 to convert TWh to MWh
     industrial_demand = (
         pd.read_csv(snakemake.input.industrial_demand, index_col=0) * 1e6
-    )
+    ) * nyears
 
     n.madd(
         "Bus",
@@ -2352,10 +2353,10 @@ def add_industry(n, costs):
             industrial_demand.loc[spatial.biomass.locations, "solid biomass"].rename(
                 index=lambda x: x + " solid biomass for industry"
             )
-            / 8760
+            / nsnapshots
         )
     else:
-        p_set = industrial_demand["solid biomass"].sum() / 8760
+        p_set = industrial_demand["solid biomass"].sum() / nsnapshots
 
     n.madd(
         "Load",
@@ -2402,7 +2403,7 @@ def add_industry(n, costs):
         unit="MWh_LHV",
     )
 
-    gas_demand = industrial_demand.loc[nodes, "methane"] / 8760.0
+    gas_demand = industrial_demand.loc[nodes, "methane"] / nsnapshots
 
     if options["gas_network"]:
         spatial_gas_demand = gas_demand.rename(index=lambda x: x + " gas for industry")
@@ -2454,7 +2455,7 @@ def add_industry(n, costs):
         suffix=" H2 for industry",
         bus=nodes + " H2",
         carrier="H2 for industry",
-        p_set=industrial_demand.loc[nodes, "hydrogen"] / 8760,
+        p_set=industrial_demand.loc[nodes, "hydrogen"] / nsnapshots,
     )
 
     shipping_hydrogen_share = get(options["shipping_hydrogen_share"], investment_year)
@@ -2470,11 +2471,11 @@ def add_industry(n, costs):
     domestic_navigation = pop_weighted_energy_totals.loc[
         nodes, "total domestic navigation"
     ].squeeze()
-    international_navigation = pd.read_csv(
+    international_navigation = (
-        snakemake.input.shipping_demand, index_col=0
+        pd.read_csv(snakemake.input.shipping_demand, index_col=0).squeeze() * nyears
-    ).squeeze()
+    )
     all_navigation = domestic_navigation + international_navigation
-    p_set = all_navigation * 1e6 / 8760
+    p_set = all_navigation * 1e6 / nsnapshots
 
     if shipping_hydrogen_share:
         oil_efficiency = options.get(
@@ -2681,7 +2682,7 @@ def add_industry(n, costs):
     )
 
     demand_factor = options.get("HVC_demand_factor", 1)
-    p_set = demand_factor * industrial_demand.loc[nodes, "naphtha"].sum() / 8760
+    p_set = demand_factor * industrial_demand.loc[nodes, "naphtha"].sum() / nsnapshots
     if demand_factor != 1:
         logger.warning(f"Changing HVC demand by {demand_factor*100-100:+.2f}%.")
 
@@ -2699,7 +2700,7 @@ def add_industry(n, costs):
         demand_factor
         * pop_weighted_energy_totals.loc[nodes, all_aviation].sum(axis=1).sum()
         * 1e6
-        / 8760
+        / nsnapshots
     )
     if demand_factor != 1:
         logger.warning(f"Changing aviation demand by {demand_factor*100-100:+.2f}%.")
@@ -2718,7 +2719,8 @@ def add_industry(n, costs):
     co2_release = ["naphtha for industry", "kerosene for aviation"]
     co2 = (
         n.loads.loc[co2_release, "p_set"].sum() * costs.at["oil", "CO2 intensity"]
-        - industrial_demand.loc[nodes, "process emission from feedstock"].sum() / 8760
+        - industrial_demand.loc[nodes, "process emission from feedstock"].sum()
+        / nsnapshots
     )
 
     n.add(
@@ -2741,12 +2743,13 @@ def add_industry(n, costs):
             for node in nodes
         ],
         carrier="low-temperature heat for industry",
-        p_set=industrial_demand.loc[nodes, "low-temperature heat"] / 8760,
+        p_set=industrial_demand.loc[nodes, "low-temperature heat"] / nsnapshots,
     )
 
     # remove today's industrial electricity demand by scaling down total electricity demand
     for ct in n.buses.country.dropna().unique():
         # TODO map onto n.bus.country
+
         loads_i = n.loads.index[
             (n.loads.index.str[:2] == ct) & (n.loads.carrier == "electricity")
         ]
@@ -2765,7 +2768,7 @@ def add_industry(n, costs):
         suffix=" industry electricity",
         bus=nodes,
         carrier="industry electricity",
-        p_set=industrial_demand.loc[nodes, "electricity"] / 8760,
+        p_set=industrial_demand.loc[nodes, "electricity"] / nsnapshots,
     )
 
     n.madd(
@@ -2782,10 +2785,10 @@ def add_industry(n, costs):
             -industrial_demand.loc[nodes, sel]
             .sum(axis=1)
             .rename(index=lambda x: x + " process emissions")
-            / 8760
+            / nsnapshots
         )
     else:
-        p_set = -industrial_demand.loc[nodes, sel].sum(axis=1).sum() / 8760
+        p_set = -industrial_demand.loc[nodes, sel].sum(axis=1).sum() / nsnapshots
 
     # this should be process emissions fossil+feedstock
     # then need load on atmosphere for feedstock emissions that are currently going to atmosphere via Link Fischer-Tropsch demand
@@ -2829,10 +2832,10 @@ def add_industry(n, costs):
                 industrial_demand.loc[spatial.ammonia.locations, "ammonia"].rename(
                     index=lambda x: x + " NH3"
                 )
-                / 8760
+                / nsnapshots
             )
         else:
-            p_set = industrial_demand["ammonia"].sum() / 8760
+            p_set = industrial_demand["ammonia"].sum() / nsnapshots
 
         n.madd(
             "Load",
@@ -2884,6 +2887,7 @@ def add_agriculture(n, costs):
     logger.info("Add agriculture, forestry and fishing sector.")
 
     nodes = pop_layout.index
+    nsnapshots = n.snapshot_weightings.generators.sum()
 
     # electricity
 
@@ -2895,7 +2899,7 @@ def add_agriculture(n, costs):
         carrier="agriculture electricity",
         p_set=pop_weighted_energy_totals.loc[nodes, "total agriculture electricity"]
         * 1e6
-        / 8760,
+        / nsnapshots,
     )
 
     # heat
@@ -2908,7 +2912,7 @@ def add_agriculture(n, costs):
         carrier="agriculture heat",
         p_set=pop_weighted_energy_totals.loc[nodes, "total agriculture heat"]
         * 1e6
-        / 8760,
+        / nsnapshots,
     )
 
     # machinery
@@ -2944,7 +2948,7 @@ def add_agriculture(n, costs):
             / efficiency_gain
             * machinery_nodal_energy
             * 1e6
-            / 8760,
+            / nsnapshots,
         )
 
     if oil_share > 0:
@@ -2953,14 +2957,14 @@ def add_agriculture(n, costs):
             ["agriculture machinery oil"],
             bus=spatial.oil.nodes,
             carrier="agriculture machinery oil",
-            p_set=oil_share * machinery_nodal_energy.sum() * 1e6 / 8760,
+            p_set=oil_share * machinery_nodal_energy.sum() * 1e6 / nsnapshots,
         )
 
         co2 = (
             oil_share
             * machinery_nodal_energy.sum()
             * 1e6
-            / 8760
+            / nsnapshots
             * costs.at["oil", "CO2 intensity"]
         )
 
@@ -3229,19 +3233,19 @@ def set_temporal_aggregation(n, opts, solver_name):
     return n
 
 
-# %%
 if __name__ == "__main__":
     if "snakemake" not in globals():
         from _helpers import mock_snakemake
 
         snakemake = mock_snakemake(
             "prepare_sector_network",
+            configfiles="test/config.overnight.yaml",
             simpl="",
             opts="",
-            clusters="37",
+            clusters="5",
             ll="v1.5",
-            sector_opts="cb40ex0-365H-T-H-B-I-A-solar+p3-dist1",
+            sector_opts="CO2L0-24H-T-H-B-I-A-solar+p3-dist1",
-            planning_horizons="2020",
+            planning_horizons="2030",
         )
 
     logging.basicConfig(level=snakemake.config["logging"]["level"])
@@ -3258,16 +3262,17 @@ if __name__ == "__main__":
     n = pypsa.Network(snakemake.input.network, override_component_attrs=overrides)
 
     pop_layout = pd.read_csv(snakemake.input.clustered_pop_layout, index_col=0)
-    Nyears = n.snapshot_weightings.generators.sum() / 8760
+    nsnapshots = n.snapshot_weightings.objective.sum()
+    nyears = nsnapshots / 8760
 
     costs = prepare_costs(
         snakemake.input.costs,
         snakemake.config["costs"],
-        Nyears,
+        nyears,
     )
 
-    pop_weighted_energy_totals = pd.read_csv(
+    pop_weighted_energy_totals = (
-        snakemake.input.pop_weighted_energy_totals, index_col=0
+        pd.read_csv(snakemake.input.pop_weighted_energy_totals, index_col=0) * nyears
     )
 
     patch_electricity_network(n)
@@ -3369,7 +3374,7 @@ if __name__ == "__main__":
         limit = float(limit.replace("p", ".").replace("m", "-"))
         break
     logger.info(f"Add CO2 limit from {limit_type}")
-    add_co2limit(n, Nyears, limit)
+    add_co2limit(n, nyears, limit)
 
     for o in opts:
         if not o[:10] == "linemaxext":

scripts/solve_network.py

@@ -31,7 +31,7 @@ def add_land_use_constraint(n, config):
         _add_land_use_constraint(n, config)
 
 
-def _add_land_use_constraint(n):
+def _add_land_use_constraint(n, config):
     # warning: this will miss existing offwind which is not classed AC-DC and has carrier 'offwind'
 
     for carrier in ["solar", "onwind", "offwind-ac", "offwind-dc"]:

test/config.myopic.yaml

@@ -58,3 +58,15 @@ solving:
     name: glpk
     options: glpk-default
   mem: 4000
+
+plotting:
+  map:
+    boundaries:
+  eu_node_location:
+    x: -5.5
+    y: 46.
+  costs_max: 1000
+  costs_threshold: 0.0000001
+  energy_max:
+  energy_min:
+  energy_threshold: 0.000001

test/config.overnight.yaml

@@ -16,7 +16,7 @@ scenario:
   clusters:
   - 5
   sector_opts:
-  - CO2L0-24H-T-H-B-I-A-solar+p3-dist1
+  - CO2L0-1H-T-H-B-I-A-solar+p3-dist1
   planning_horizons:
   - 2030
 
@@ -59,3 +59,15 @@ solving:
     name: glpk
     options: glpk-default
   mem: 4000
+
+plotting:
+  map:
+    boundaries:
+  eu_node_location:
+    x: -5.5
+    y: 46.
+  costs_max: 1000
+  costs_threshold: 0.0000001
+  energy_max:
+  energy_min:
+  energy_threshold: 0.000001

test/config.test1.yaml

@@ -65,3 +65,16 @@ solving:
   solver:
     name: glpk
     options: "glpk-default"
+
+
+plotting:
+  map:
+    boundaries:
+  eu_node_location:
+    x: -5.5
+    y: 46.
+  costs_max: 1000
+  costs_threshold: 0.0000001
+  energy_max:
+  energy_min:
+  energy_threshold: 0.000001