diff --git a/Snakefile b/Snakefile index ba99f54f..3632e498 100644 --- a/Snakefile +++ b/Snakefile @@ -81,8 +81,14 @@ rule build_simplified_population_layouts: if config["sector"]["gas_network"]: + datafiles = [ + "IGGIELGN_LNGs.geojson", + "IGGIELGN_BorderPoints.geojson", + "IGGIELGN_Productions.geojson", + ] + rule retrieve_gas_infrastructure_data: - output: "data/gas_network/scigrid-gas/data/IGGIELGN_LNGs.csv" + output: expand("data/gas_network/scigrid-gas/data/{files}", files=datafiles) script: 'scripts/retrieve_gas_infrastructure_data.py' rule build_gas_network: @@ -93,20 +99,20 @@ if config["sector"]["gas_network"]: resources: mem_mb=4000 script: "scripts/build_gas_network.py" - rule build_gas_import_locations: + rule build_gas_input_locations: input: - lng="data/gas_network/scigrid-gas/data/IGGIELGN_LNGs.geojson" - entry="data/gas_network/scigrid-gas/data/IGGIELGN_BorderPoints.geojson" - production="data/gas_network/scigrid-gas/data/IGGIELGN_Productions.geojson" + lng="data/gas_network/scigrid-gas/data/IGGIELGN_LNGs.geojson", + entry="data/gas_network/scigrid-gas/data/IGGIELGN_BorderPoints.geojson", + production="data/gas_network/scigrid-gas/data/IGGIELGN_Productions.geojson", regions_onshore=pypsaeur("resources/regions_onshore_elec_s{simpl}_{clusters}.geojson"), output: - gas_input_nodes="resources/gas_input_nodes_s{simpl}_{clusters}.csv" + gas_input_nodes="resources/gas_input_locations_s{simpl}_{clusters}.csv" resources: mem_mb=2000, - script: "scripts/build_gas_import_locations.py" + script: "scripts/build_gas_input_locations.py" rule cluster_gas_network: input: - cleaned_gas_network="data/gas_network/gas_network_dataset.csv", + cleaned_gas_network="resources/gas_network.csv", regions_onshore=pypsaeur("resources/regions_onshore_elec_s{simpl}_{clusters}.geojson"), regions_offshore=pypsaeur("resources/regions_offshore_elec_s{simpl}_{clusters}.geojson") output: @@ -114,7 +120,7 @@ if config["sector"]["gas_network"]: resources: mem_mb=4000 script: "scripts/cluster_gas_network.py" - gas_infrastructure = {**rules.cluster_gas_network.output, **rules.build_gas_import_locations.output} + gas_infrastructure = {**rules.cluster_gas_network.output, **rules.build_gas_input_locations.output} else: gas_infrastructure = {} diff --git a/config.default.yaml b/config.default.yaml index 61a70fac..59a27753 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -243,12 +243,14 @@ sector: electricity_distribution_grid: false electricity_distribution_grid_cost_factor: 1.0 #multiplies cost in data/costs.csv electricity_grid_connection: true # only applies to onshore wind and utility PV + H2_network: true gas_network: true H2_retrofit: true # if set to True existing gas pipes can be retrofitted to H2 pipes # according to hydrogen backbone strategy (April, 2020) p.15 # https://gasforclimate2050.eu/wp-content/uploads/2020/07/2020_European-Hydrogen-Backbone_Report.pdf # 60% of original natural gas capacity could be used in cost-optimal case as H2 capacity H2_retrofit_capacity_per_CH4: 0.6 # ratio for H2 capacity per original CH4 capacity of retrofitted pipelines + gas_network_connectivity_upgrade: 3 # https://networkx.org/documentation/stable/reference/algorithms/generated/networkx.algorithms.connectivity.edge_augmentation.k_edge_augmentation.html#networkx.algorithms.connectivity.edge_augmentation.k_edge_augmentation gas_distribution_grid: true gas_distribution_grid_cost_factor: 1.0 #multiplies cost in data/costs.csv biomass_transport: false # biomass transport between nodes @@ -449,6 +451,7 @@ plotting: gas boilers: '#db6a25' gas boiler marginal: '#db6a25' gas: '#e05b09' + fossil gas: '#e05b09' natural gas: '#e05b09' CCGT: '#a85522' CCGT marginal: '#a85522' @@ -457,6 +460,7 @@ plotting: gas for industry: '#853403' gas for industry CC: '#692e0a' gas pipeline: '#ebbca0' + gas pipeline new: '#a87c62' # oil oil: '#c9c9c9' oil boiler: '#adadad' @@ -546,6 +550,7 @@ plotting: H2 storage: '#bf13a0' land transport fuel cell: '#6b3161' H2 pipeline: '#f081dc' + H2 pipeline retrofitted: '#ba99b5' H2 Fuel Cell: '#c251ae' H2 Electrolysis: '#ff29d9' # syngas diff --git a/scripts/build_biomass_potentials.py b/scripts/build_biomass_potentials.py index 68d87808..2b6812f3 100644 --- a/scripts/build_biomass_potentials.py +++ b/scripts/build_biomass_potentials.py @@ -175,7 +175,7 @@ def convert_nuts2_to_regions(bio_nuts2, regions): # calculate area of nuts2 regions bio_nuts2["area_nuts2"] = area(bio_nuts2) - overlay = gpd.overlay(regions, bio_nuts2) + overlay = gpd.overlay(regions, bio_nuts2, keep_geom_type=True) # calculate share of nuts2 area inside region overlay["share"] = area(overlay) / overlay["area_nuts2"] diff --git a/scripts/build_gas_import_locations.py b/scripts/build_gas_input_locations.py similarity index 82% rename from scripts/build_gas_import_locations.py rename to scripts/build_gas_input_locations.py index 8ee6d577..b19c1efd 100644 --- a/scripts/build_gas_import_locations.py +++ b/scripts/build_gas_input_locations.py @@ -1,5 +1,5 @@ """ -Build import locations for fossil gas from entry-points and LNG terminals. +Build import locations for fossil gas from entry-points, LNG terminals and production sites. """ import logging @@ -16,11 +16,8 @@ def read_scigrid_gas(fn): return df -def build_gas_input_locations(lng_fn, entry_fn, prod_fn): +def build_gas_input_locations(lng_fn, entry_fn, prod_fn, countries): - countries = snakemake.config["countries"] - countries[countries.index('GB')] = 'UK' - # LNG terminals lng = read_scigrid_gas(lng_fn) @@ -37,7 +34,8 @@ def build_gas_input_locations(lng_fn, entry_fn, prod_fn): prod = read_scigrid_gas(prod_fn) prod = prod.loc[ (prod.geometry.y > 35) & - (prod.geometry.x < 30) + (prod.geometry.x < 30) & + (prod.country_code != "DE") ] return gpd.GeoDataFrame( @@ -60,10 +58,13 @@ if __name__ == "__main__": onshore_regions = gpd.read_file(snakemake.input.regions_onshore).set_index('name') + countries = onshore_regions.index.str[:2].unique().str.replace("GB", "UK") + gas_input_locations = build_gas_input_locations( snakemake.input.lng, snakemake.input.entry, - snakemake.input.production + snakemake.input.production, + countries ) # recommended to use projected CRS rather than geographic CRS diff --git a/scripts/cluster_gas_network.py b/scripts/cluster_gas_network.py index 00dd165d..9f192a92 100755 --- a/scripts/cluster_gas_network.py +++ b/scripts/cluster_gas_network.py @@ -76,8 +76,9 @@ def aggregate_parallel_pipes(df): "diameter_mm": "mean", "length": 'mean', 'tags': ' '.join, + "p_min_pu": 'min', } - df = df.groupby(df.index).agg(strategies) + return df.groupby(df.index).agg(strategies) if __name__ == "__main__": @@ -105,6 +106,6 @@ if __name__ == "__main__": gas_network = build_clustered_gas_network(df, bus_regions) reindex_pipes(gas_network) - aggregate_parallel_pipes(gas_network) + gas_network = aggregate_parallel_pipes(gas_network) gas_network.to_csv(snakemake.output.clustered_gas_network) \ No newline at end of file diff --git a/scripts/plot_network.py b/scripts/plot_network.py index 2a6073d7..0c1049c0 100644 --- a/scripts/plot_network.py +++ b/scripts/plot_network.py @@ -236,8 +236,6 @@ def plot_h2_map(network): linewidth_factor = 1e4 # MW below which not drawn line_lower_threshold = 1e3 - bus_color = "m" - link_color = "c" # Drop non-electric buses so they don't clutter the plot n.buses.drop(n.buses.index[n.buses.carrier != "AC"], inplace=True) @@ -251,38 +249,68 @@ def plot_h2_map(network): n.links.drop(n.links.index[~n.links.carrier.str.contains("H2 pipeline")], inplace=True) - link_widths = n.links.p_nom_opt / linewidth_factor - link_widths[n.links.p_nom_opt < line_lower_threshold] = 0. - link_color = n.links.carrier.map({"H2 pipeline":"red", - "H2 pipeline retrofitted": "blue"}) + h2_new = n.links.loc[n.links.carrier=="H2 pipeline", "p_nom_opt"] + h2_retro = n.links.loc[n.links.carrier=='H2 pipeline retrofitted'] + + positive_order = h2_retro.bus0 < h2_retro.bus1 + h2_retro_p = h2_retro[positive_order] + swap_buses = {"bus0": "bus1", "bus1": "bus0"} + h2_retro_n = h2_retro[~positive_order].rename(columns=swap_buses) + h2_retro = pd.concat([h2_retro_p, h2_retro_n]) + + h2_retro.index = h2_retro.apply( + lambda x: f"H2 pipeline {x.bus0.replace(' H2', '')} -> {x.bus1.replace(' H2', '')}", + axis=1 + ) + + h2_retro = h2_retro["p_nom_opt"] + + link_widths_total = (h2_new + h2_retro) / linewidth_factor + link_widths_total = link_widths_total.groupby(level=0).sum().reindex(n.links.index).fillna(0.) + link_widths_total[n.links.p_nom_opt < line_lower_threshold] = 0. + + retro = n.links.p_nom_opt.where(n.links.carrier=='H2 pipeline retrofitted', other=0.) + link_widths_retro = retro / linewidth_factor + link_widths_retro[n.links.p_nom_opt < line_lower_threshold] = 0. n.links.bus0 = n.links.bus0.str.replace(" H2", "") n.links.bus1 = n.links.bus1.str.replace(" H2", "") - print(link_widths.sort_values()) - - print(n.links[["bus0", "bus1"]]) - fig, ax = plt.subplots( figsize=(7, 6), subplot_kw={"projection": ccrs.PlateCarree()} ) + + bus_colors = { + "H2 Electrolysis": "m", + "H2 Fuel Cell": "slateblue" + } n.plot( bus_sizes=bus_sizes, - bus_colors={"H2 Electrolysis": bus_color, - "H2 Fuel Cell": "slateblue"}, - link_colors=link_color, - link_widths=link_widths, + bus_colors=bus_colors, + link_colors='#afc6c7', + link_widths=link_widths_total, branch_components=["Link"], - ax=ax, **map_opts + ax=ax, + **map_opts + ) + + n.plot( + geomap=False, + bus_sizes=0, + link_colors='#72d3d6', + link_widths=link_widths_retro, + branch_components=["Link"], + ax=ax, + **map_opts ) handles = make_legend_circles_for( [50000, 10000], scale=bus_size_factor, - facecolor=bus_color + facecolor='k' ) labels = ["{} GW".format(s) for s in (50, 10)] @@ -330,261 +358,127 @@ def plot_ch4_map(network): n = network.copy() - supply_energy = get_nodal_balance().droplevel([0,1]).sort_index() - if "gas pipeline" not in n.links.carrier.unique(): return assign_location(n) - bus_size_factor = 1e7 + bus_size_factor = 4e7 linewidth_factor = 1e4 # MW below which not drawn - line_lower_threshold = 5e3 - bus_color = "maroon" - link_color = "lightcoral" + line_lower_threshold = 500 # Drop non-electric buses so they don't clutter the plot n.buses.drop(n.buses.index[n.buses.carrier != "AC"], inplace=True) elec = n.generators[n.generators.carrier=="gas"].index - methanation_i = n.links[n.links.carrier.isin(["helmeth", "Sabatier"])].index - - bus_sizes = n.generators_t.p.loc[:,elec].mul(n.snapshot_weightings, axis=0).sum().groupby(n.generators.loc[elec,"bus"]).sum() / bus_size_factor + + bus_sizes = n.generators_t.p.loc[:,elec].mul(n.snapshot_weightings.generators, axis=0).sum().groupby(n.generators.loc[elec,"bus"]).sum() / bus_size_factor bus_sizes.rename(index=lambda x: x.replace(" gas", ""), inplace=True) bus_sizes = bus_sizes.reindex(n.buses.index).fillna(0) - bus_sizes.index = pd.MultiIndex.from_product( - [bus_sizes.index, ["fossil gas"]]) + bus_sizes.index = pd.MultiIndex.from_product([bus_sizes.index, ["fossil gas"]]) - methanation = abs(n.links_t.p1.loc[:,methanation_i].mul(n.snapshot_weightings, axis=0)).sum().groupby(n.links.loc[methanation_i,"bus1"]).sum() / bus_size_factor + methanation_i = n.links[n.links.carrier.isin(["helmeth", "Sabatier"])].index + methanation = abs(n.links_t.p1.loc[:,methanation_i].mul(n.snapshot_weightings.generators, axis=0)).sum().groupby(n.links.loc[methanation_i,"bus1"]).sum() / bus_size_factor methanation = methanation.groupby(methanation.index).sum().rename(index=lambda x: x.replace(" gas", "")) # make a fake MultiIndex so that area is correct for legend - methanation.index = pd.MultiIndex.from_product( - [methanation.index, ["methanation"]]) + methanation.index = pd.MultiIndex.from_product([methanation.index, ["methanation"]]) biogas_i = n.stores[n.stores.carrier=="biogas"].index - biogas = n.stores_t.p.loc[:,biogas_i].mul(n.snapshot_weightings, axis=0).sum().groupby(n.stores.loc[biogas_i,"bus"]).sum() / bus_size_factor + biogas = n.stores_t.p.loc[:,biogas_i].mul(n.snapshot_weightings.generators, axis=0).sum().groupby(n.stores.loc[biogas_i,"bus"]).sum() / bus_size_factor biogas = biogas.groupby(biogas.index).sum().rename(index=lambda x: x.replace(" biogas", "")) # make a fake MultiIndex so that area is correct for legend - biogas.index = pd.MultiIndex.from_product( - [biogas.index, ["biogas"]]) + biogas.index = pd.MultiIndex.from_product([biogas.index, ["biogas"]]) bus_sizes = pd.concat([bus_sizes, methanation, biogas]) bus_sizes.sort_index(inplace=True) - n.links.drop(n.links.index[n.links.carrier != "gas pipeline"], inplace=True) + to_remove = n.links.index[~n.links.carrier.str.contains("gas pipeline")] + n.links.drop(to_remove, inplace=True) - link_widths = n.links.p_nom_opt / linewidth_factor + link_widths = n.links.p_nom_opt / linewidth_factor link_widths[n.links.p_nom_opt < line_lower_threshold] = 0. + link_widths_orig = n.links.p_nom / linewidth_factor + link_widths_orig[n.links.p_nom < line_lower_threshold] = 0. + + link_color = n.links.carrier.map({"gas pipeline": "lightcoral", + "gas pipeline new": "red"}) + n.links.bus0 = n.links.bus0.str.replace(" gas", "") n.links.bus1 = n.links.bus1.str.replace(" gas", "") - print(link_widths.sort_values()) + bus_colors = { + "fossil gas": 'maroon', + "methanation": "steelblue", + "biogas": "seagreen" + } - print(n.links[["bus0", "bus1"]]) + fig, ax = plt.subplots(figsize=(7,6), subplot_kw={"projection": ccrs.PlateCarree()}) - fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()}) + n.plot( + bus_sizes=bus_sizes, + bus_colors=bus_colors, + link_colors='lightgrey', + link_widths=link_widths_orig, + branch_components=["Link"], + ax=ax, + **map_opts + ) - fig.set_size_inches(7, 6) - - n.plot(bus_sizes=bus_sizes, - bus_colors={"fossil gas": bus_color, - "methanation": "steelblue", - "biogas": "seagreen"}, - link_colors=link_color, - link_widths=link_widths, - branch_components=["Link"], - ax=ax, boundaries=(-10, 30, 34, 70)) + n.plot( + geomap=False, + ax=ax, + bus_sizes=0., + link_colors=link_color, + link_widths=link_widths, + branch_components=["Link"], + **map_opts + ) handles = make_legend_circles_for( - [200, 1000], scale=bus_size_factor, facecolor=bus_color) + [200000, 1000000], + scale=bus_size_factor, + facecolor='k' + ) labels = ["{} MW".format(s) for s in (200, 1000)] - l2 = ax.legend(handles, labels, - loc="upper left", bbox_to_anchor=(0.01, 1.01), - labelspacing=1.0, - framealpha=1., - title='Biomass potential', - handler_map=make_handler_map_to_scale_circles_as_in(ax)) + + l2 = ax.legend( + handles, labels, + loc="upper left", + bbox_to_anchor=(-0.03, 1.01), + labelspacing=1.0, + frameon=False, + title='gas generation', + handler_map=make_handler_map_to_scale_circles_as_in(ax) + ) + ax.add_artist(l2) handles = [] labels = [] for s in (50, 10): - handles.append(plt.Line2D([0], [0], color=link_color, - linewidth=s * 1e3 / linewidth_factor)) + handles.append(plt.Line2D([0], [0], color="k", linewidth=s * 1e3 / linewidth_factor)) labels.append("{} GW".format(s)) - l1_1 = ax.legend(handles, labels, - loc="upper left", bbox_to_anchor=(0.30, 1.01), - framealpha=1, - labelspacing=0.8, handletextpad=1.5, - title='CH4 pipeline capacity') + + l1_1 = ax.legend( + handles, labels, + loc="upper left", + bbox_to_anchor=(0.28, 1.01), + frameon=False, + labelspacing=0.8, + handletextpad=1.5, + title='gas pipeline capacity' + ) + ax.add_artist(l1_1) - fig.savefig(snakemake.output.map.replace("-costs-all","-ch4_network"), transparent=True, - bbox_inches="tight") + fig.savefig( + snakemake.output.map.replace("-costs-all","-ch4_network"), + bbox_inches="tight" + ) - ################################################## - supply_energy.drop("gas pipeline", level=1, inplace=True) - supply_energy = supply_energy[abs(supply_energy)>5] - supply_energy.rename(index=lambda x: x.replace(" gas",""), level=0, inplace=True) - - - demand = supply_energy[supply_energy<0].groupby(level=[0,1]).sum() - supply = supply_energy[supply_energy>0].groupby(level=[0,1]).sum() - - #### DEMAND ####################################### - bus_size_factor = 2e7 - bus_sizes = abs(demand) / bus_size_factor - - fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()}) - - fig.set_size_inches(7, 6) - - n.plot(bus_sizes=bus_sizes, - bus_colors={"CHP": "r", - "OCGT": "wheat", - "SMR": "darkkhaki", - "SMR CC": "tan", - "gas boiler": "orange", - "gas for industry": "grey", - 'gas for industry CC': "lightgrey"}, - link_colors=link_color, - link_widths=link_widths, - branch_components=["Link"], - ax=ax, boundaries=(-10, 30, 34, 70)) - - handles = make_legend_circles_for( - [10e6, 20e6], scale=bus_size_factor, facecolor=bus_color) - labels = ["{} TWh".format(s) for s in (10, 20)] - l2 = ax.legend(handles, labels, - loc="upper left", bbox_to_anchor=(0.01, 1.01), - labelspacing=1.0, - framealpha=1., - title='CH4 demand', - handler_map=make_handler_map_to_scale_circles_as_in(ax)) - ax.add_artist(l2) - - handles = [] - labels = [] - - for s in (50, 10): - handles.append(plt.Line2D([0], [0], color=link_color, - linewidth=s * 1e3 / linewidth_factor)) - labels.append("{} GW".format(s)) - l1_1 = ax.legend(handles, labels, - loc="upper left", bbox_to_anchor=(0.30, 1.01), - framealpha=1, - labelspacing=0.8, handletextpad=1.5, - title='CH4 pipeline capacity') - ax.add_artist(l1_1) - - fig.savefig(snakemake.output.map.replace("-costs-all","-ch4_demand"), transparent=True, - bbox_inches="tight") - - - #### SUPPLY ####################################### - bus_size_factor = 2e7 - bus_sizes = supply / bus_size_factor - - fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()}) - - fig.set_size_inches(7, 6) - - n.plot(bus_sizes=bus_sizes, - bus_colors={"gas": "maroon", - "methanation": "steelblue", - "helmeth": "slateblue", - "biogas": "seagreen"}, - link_colors=link_color, - link_widths=link_widths, - branch_components=["Link"], - ax=ax, boundaries=(-10, 30, 34, 70)) - - handles = make_legend_circles_for( - [10e6, 20e6], scale=bus_size_factor, facecolor="black") - labels = ["{} TWh".format(s) for s in (10, 20)] - l2 = ax.legend(handles, labels, - loc="upper left", bbox_to_anchor=(0.01, 1.01), - labelspacing=1.0, - framealpha=1., - title='CH4 supply', - handler_map=make_handler_map_to_scale_circles_as_in(ax)) - ax.add_artist(l2) - - handles = [] - labels = [] - - for s in (50, 10): - handles.append(plt.Line2D([0], [0], color=link_color, - linewidth=s * 1e3 / linewidth_factor)) - labels.append("{} GW".format(s)) - l1_1 = ax.legend(handles, labels, - loc="upper left", bbox_to_anchor=(0.30, 1.01), - framealpha=1, - labelspacing=0.8, handletextpad=1.5, - title='CH4 pipeline capacity') - ax.add_artist(l1_1) - - fig.savefig(snakemake.output.map.replace("-costs-all","-ch4_supply"), transparent=True, - bbox_inches="tight") - - ########################################################################### - net = pd.concat([demand.groupby(level=0).sum().rename("demand"), - supply.groupby(level=0).sum().rename("supply")], axis=1).sum(axis=1) - mask = net>0 - net_importer = net.mask(mask).rename("net_importer") - net_exporter = net.mask(~mask).rename("net_exporter") - - bus_size_factor = 2e7 - linewidth_factor = 1e-1 - bus_sizes = pd.concat([abs(net_importer), net_exporter], axis=1).fillna(0).stack() / bus_size_factor - - link_widths = abs(n.links_t.p0).max().loc[n.links.index] / n.links.p_nom_opt - link_widths /= linewidth_factor - - - fig, ax = plt.subplots(subplot_kw={"projection": ccrs.PlateCarree()}) - - fig.set_size_inches(7, 6) - - n.plot(bus_sizes=bus_sizes, - bus_colors={"net_importer": "r", - "net_exporter": "darkgreen", - }, - link_colors="lightgrey", - link_widths=link_widths, - branch_components=["Link"], - ax=ax, boundaries=(-10, 30, 34, 70)) - - handles = make_legend_circles_for( - [10e6, 20e6], scale=bus_size_factor, facecolor="black") - labels = ["{} TWh".format(s) for s in (10, 20)] - l2 = ax.legend(handles, labels, - loc="upper left", bbox_to_anchor=(0.01, 1.01), - labelspacing=1.0, - framealpha=1., - title='Net Import/Export', - handler_map=make_handler_map_to_scale_circles_as_in(ax)) - ax.add_artist(l2) - - handles = [] - labels = [] - - for s in (0.5, 1): - handles.append(plt.Line2D([0], [0], color="lightgrey", - linewidth=s / linewidth_factor)) - labels.append("{} per unit".format(s)) - l1_1 = ax.legend(handles, labels, - loc="upper left", bbox_to_anchor=(0.30, 1.01), - framealpha=1, - labelspacing=0.8, handletextpad=1.5, - title='maximum used CH4 pipeline capacity') - ax.add_artist(l1_1) - - fig.savefig(snakemake.output.map.replace("-costs-all","-ch4_net_balance"), transparent=True, - bbox_inches="tight") def plot_map_without(network): @@ -785,51 +679,6 @@ def plot_series(network, carrier="AC", name="test"): transparent=True) -def get_nodal_balance(carrier="gas"): - - bus_map = (n.buses.carrier == carrier) - bus_map.at[""] = False - supply_energy = pd.Series(dtype="float64") - - for c in n.iterate_components(n.one_port_components): - - items = c.df.index[c.df.bus.map(bus_map).fillna(False)] - - if len(items) == 0: - continue - - s = round(c.pnl.p.multiply(n.snapshot_weightings,axis=0).sum().multiply(c.df['sign']).loc[items] - .groupby([c.df.bus, c.df.carrier]).sum()) - s = pd.concat([s], keys=[c.list_name]) - s = pd.concat([s], keys=[carrier]) - - supply_energy = supply_energy.reindex(s.index.union(supply_energy.index)) - supply_energy.loc[s.index] = s - - - for c in n.iterate_components(n.branch_components): - - for end in [col[3:] for col in c.df.columns if col[:3] == "bus"]: - - items = c.df.index[c.df["bus" + str(end)].map(bus_map,na_action=False)] - - if len(items) == 0: - continue - - s = ((-1)*c.pnl["p"+end][items].multiply(n.snapshot_weightings,axis=0).sum() - .groupby([c.df.loc[items,'bus{}'.format(end)], c.df.loc[items,'carrier']]).sum()) - s.index = s.index - s = pd.concat([s], keys=[c.list_name]) - s = pd.concat([s], keys=[carrier]) - - supply_energy = supply_energy.reindex(s.index.union(supply_energy.index)) - - supply_energy.loc[s.index] = s - - supply_energy = supply_energy.rename(index=lambda x: rename_techs(x), level=3) - return supply_energy - - if __name__ == "__main__": if 'snakemake' not in globals(): from helper import mock_snakemake diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index df8bc081..1a58d01e 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -8,6 +8,7 @@ import pytz import pandas as pd import numpy as np import xarray as xr +import networkx as nx from itertools import product from scipy.stats import beta @@ -16,6 +17,10 @@ from vresutils.costdata import annuity from build_energy_totals import build_eea_co2, build_eurostat_co2, build_co2_totals from helper import override_component_attrs +from networkx.algorithms.connectivity.edge_augmentation import k_edge_augmentation +from networkx.algorithms import complement +from pypsa.geo import haversine_pts + import logging logger = logging.getLogger(__name__) @@ -131,40 +136,6 @@ def get(item, investment_year=None): return item -def create_network_topology(n, prefix, connector=" -> "): - """ - Create a network topology like the power transmission network. - - Parameters - ---------- - n : pypsa.Network - prefix : str - connector : str - - Returns - ------- - pd.DataFrame with columns bus0, bus1 and length - """ - - ln_attrs = ["bus0", "bus1", "length"] - lk_attrs = ["bus0", "bus1", "length", "underwater_fraction"] - - candidates = pd.concat([ - n.lines[ln_attrs], - n.links.loc[n.links.carrier == "DC", lk_attrs] - ]).fillna(0) - - positive_order = candidates.bus0 < candidates.bus1 - candidates_p = candidates[positive_order] - swap_buses = {"bus0": "bus1", "bus1": "bus0"} - candidates_n = candidates[~positive_order].rename(columns=swap_buses) - candidates = pd.concat([candidates_p, candidates_n]) - - topo = candidates.groupby(["bus0", "bus1"], as_index=False).mean() - topo.index = topo.apply(lambda c: prefix + c.bus0 + connector + c.bus1, axis=1) - return topo - - def co2_emissions_year(countries, opts, year): """ Calculate CO2 emissions in one specific year (e.g. 1990 or 2018). @@ -252,14 +223,21 @@ def add_lifetime_wind_solar(n, costs): n.generators.loc[gen_i, "lifetime"] = costs.at[carrier, 'lifetime'] -def create_network_topology(n, prefix, connector=" -> ", bidirectional=True): +def haversine(p): + coord0 = n.buses.loc[p.bus0, ['x', 'y']].values + coord1 = n.buses.loc[p.bus1, ['x', 'y']].values + return 1.5 * haversine_pts(coord0, coord1) + + +def create_network_topology(n, prefix, carriers=["DC"], connector=" -> ", bidirectional=True): """ - Create a network topology like the power transmission network. + Create a network topology from transmission lines and link carrier selection. Parameters ---------- n : pypsa.Network prefix : str + carriers : list-like connector : str bidirectional : bool, default True True: one link for each connection @@ -267,7 +245,7 @@ def create_network_topology(n, prefix, connector=" -> ", bidirectional=True): Returns ------- - pd.DataFrame with columns bus0, bus1 and length + pd.DataFrame with columns bus0, bus1, length, underwater_fraction """ ln_attrs = ["bus0", "bus1", "length"] @@ -275,9 +253,13 @@ def create_network_topology(n, prefix, connector=" -> ", bidirectional=True): candidates = pd.concat([ n.lines[ln_attrs], - n.links.loc[n.links.carrier == "DC", lk_attrs] + n.links.loc[n.links.carrier.isin(carriers), lk_attrs] ]).fillna(0) + # base network topology purely on location not carrier + candidates["bus0"] = candidates.bus0.map(n.buses.location) + candidates["bus1"] = candidates.bus1.map(n.buses.location) + positive_order = candidates.bus0 < candidates.bus1 candidates_p = candidates[positive_order] swap_buses = {"bus0": "bus1", "bus1": "bus0"} @@ -1110,25 +1092,16 @@ def add_storage_and_grids(n, costs): logger.info("Add gas network") - cols = [ - "bus0", - "bus1", - "p_min_pu", - "p_nom", - "tags", - "length" - "build_year" - ] fn = snakemake.input.clustered_gas_network - gas_pipes = pd.read_csv(fn, usecols=cols, index_col=0) + gas_pipes = pd.read_csv(fn, index_col=0) if options["H2_retrofit"]: - gas_pipes["p_nom_max"] = gas_pipes.gas_pipes.p_nom + gas_pipes["p_nom_max"] = gas_pipes.p_nom gas_pipes["p_nom_min"] = 0. gas_pipes["capital_cost"] = 0. else: gas_pipes["p_nom_max"] = np.inf - gas_pipes["p_nom_min"] = gas_pipes.gas_pipes.p_nom + gas_pipes["p_nom_min"] = gas_pipes.p_nom gas_pipes["capital_cost"] = gas_pipes.length * costs.at['CH4 (g) pipeline', 'fixed'] n.madd("Link", @@ -1144,11 +1117,12 @@ def add_storage_and_grids(n, costs): capital_cost=gas_pipes.capital_cost, tags=gas_pipes.tags, carrier="gas pipeline", - lifetime=50 + lifetime=costs.at['CH4 (g) pipeline', 'lifetime'] ) # remove fossil generators where there is neither # production, LNG terminal, nor entry-point beyond system scope + fn = snakemake.input.gas_input_nodes gas_input_nodes = pd.read_csv(fn, index_col=0, squeeze=True).values remove_i = n.generators[ @@ -1157,8 +1131,41 @@ def add_storage_and_grids(n, costs): ].index n.generators.drop(remove_i, inplace=True) - # TODO candidate gas network topology + # add candidates for new gas pipelines to achieve full connectivity + G = nx.Graph() + + gas_buses = n.buses.loc[n.buses.carrier=='gas', 'location'] + G.add_nodes_from(np.unique(gas_buses.values)) + + sel = gas_pipes.p_nom > 1500 + attrs = ["bus0", "bus1", "length"] + G.add_weighted_edges_from(gas_pipes.loc[sel, attrs].values) + + # find all complement edges + complement_edges = pd.DataFrame(complement(G).edges, columns=["bus0", "bus1"]) + complement_edges["length"] = complement_edges.apply(haversine, axis=1) + + # apply k_edge_augmentation weighted by length of complement edges + k_edge = options.get("gas_network_connectivity_upgrade", 3) + augmentation = k_edge_augmentation(G, k_edge, avail=complement_edges.values) + new_gas_pipes = pd.DataFrame(augmentation, columns=["bus0", "bus1"]) + new_gas_pipes["length"] = new_gas_pipes.apply(haversine, axis=1) + + new_gas_pipes.index = new_gas_pipes.apply( + lambda x: f"gas pipeline new {x.bus0} <-> {x.bus1}", axis=1) + + n.madd("Link", + new_gas_pipes.index, + bus0=new_gas_pipes.bus0 + " gas", + bus1=new_gas_pipes.bus1 + " gas", + p_min_pu=-1, # new gas pipes are bidirectional + p_nom_extendable=True, + length=new_gas_pipes.length, + capital_cost=new_gas_pipes.length * costs.at['CH4 (g) pipeline', 'fixed'], + carrier="gas pipeline new", + lifetime=costs.at['CH4 (g) pipeline', 'lifetime'] + ) # retroftting existing CH4 pipes to H2 pipes if options["gas_network"] and options["H2_retrofit"]: @@ -1172,49 +1179,32 @@ def add_storage_and_grids(n, costs): h2_pipes.index, bus0=h2_pipes.bus0 + " H2", bus1=h2_pipes.bus1 + " H2", - p_min_pu=-1., # allow that all H2 pipelines can be used in other direction - p_nom_max=h2_pipes.pipe_capacity_MW * options["H2_retrofit_capacity_per_CH4"], + p_min_pu=-1., # allow that all H2 retrofit pipelines can be used in both directions + p_nom_max=h2_pipes.p_nom * options["H2_retrofit_capacity_per_CH4"], p_nom_extendable=True, - length=h2_pipes.length_km, - capital_cost=costs.at['H2 (g) pipeline repurposed', 'fixed'] * h2_pipes.length_km, - type=gas_pipes.num_parallel, - tags=h2_pipes.id, + length=h2_pipes.length, + capital_cost=costs.at['H2 (g) pipeline repurposed', 'fixed'] * h2_pipes.length, + tags=h2_pipes.tags, carrier="H2 pipeline retrofitted", - lifetime=50 + lifetime=costs.at['H2 (g) pipeline repurposed', 'lifetime'] ) - attrs = ["bus0", "bus1", "length"] - h2_links = pd.DataFrame(columns=attrs) + if options.get("H2_network", True): - lines_sel = n.lines[attrs] - links_sel = n.links.loc[n.links.carrier.isin(["DC", "gas pipeline"]), attrs] + h2_pipes = create_network_topology(n, "H2 pipeline ", carriers=["DC", "gas pipeline"]) - candidates = pd.concat({ - "lines": lines_sel, - "links": links_sel, - }) - - for candidate in candidates.index: - buses = [candidates.at[candidate, "bus0"], candidates.at[candidate, "bus1"]] - buses.sort() - name = f"H2 pipeline {buses[0]} -> {buses[1]}" - if name not in h2_links.index: - h2_links.at[name, "bus0"] = buses[0] - h2_links.at[name, "bus1"] = buses[1] - h2_links.at[name, "length"] = candidates.at[candidate, "length"] - - # TODO Add efficiency losses - n.madd("Link", - h2_links.index, - bus0=h2_links.bus0.values + " H2", - bus1=h2_links.bus1.values + " H2", - p_min_pu=-1, - p_nom_extendable=True, - length=h2_links.length.values, - capital_cost=costs.at['H2 (g) pipeline', 'fixed'] * h2_links.length.values, - carrier="H2 pipeline", - lifetime=costs.at['H2 (g) pipeline', 'lifetime'] - ) + # TODO Add efficiency losses + n.madd("Link", + h2_pipes.index, + bus0=h2_pipes.bus0.values + " H2", + bus1=h2_pipes.bus1.values + " H2", + p_min_pu=-1, + p_nom_extendable=True, + length=h2_pipes.length.values, + capital_cost=costs.at['H2 (g) pipeline', 'fixed'] * h2_pipes.length.values, + carrier="H2 pipeline", + lifetime=costs.at['H2 (g) pipeline', 'lifetime'] + ) n.add("Carrier", "battery")