From 6a00d5bfca4f803e203d1ecd07e377a68fc9c917 Mon Sep 17 00:00:00 2001 From: Fabian Neumann Date: Wed, 3 Nov 2021 20:34:43 +0100 Subject: [PATCH] revision gas infrastructure representation --- Snakefile | 54 +++- config.default.yaml | 1 - scripts/build_gas_import_locations.py | 76 +++++ scripts/build_gas_network.py | 289 ++++++-------------- scripts/cluster_gas_network.py | 110 ++++++++ scripts/helper.py | 12 + scripts/plot_network.py | 6 +- scripts/prepare_sector_network.py | 199 +++++++------- scripts/retrieve_gas_infrastructure_data.py | 36 +++ scripts/solve_network.py | 4 +- 10 files changed, 465 insertions(+), 322 deletions(-) create mode 100644 scripts/build_gas_import_locations.py mode change 100755 => 100644 scripts/build_gas_network.py create mode 100755 scripts/cluster_gas_network.py create mode 100644 scripts/retrieve_gas_infrastructure_data.py diff --git a/Snakefile b/Snakefile index 1bec5683..ba99f54f 100644 --- a/Snakefile +++ b/Snakefile @@ -78,16 +78,46 @@ rule build_simplified_population_layouts: benchmark: "benchmarks/build_clustered_population_layouts/s{simpl}" script: "scripts/build_clustered_population_layouts.py" -rule build_gas_network: - input: - gas_network="data/gas_network/gas_network_dataset.csv", - country_shapes=pypsaeur("resources/country_shapes.geojson"), - regions_onshore=pypsaeur("resources/regions_onshore_elec_s{simpl}_{clusters}.geojson"), - regions_offshore=pypsaeur("resources/regions_offshore_elec_s{simpl}_{clusters}.geojson") - output: - clustered_gas_network="resources/gas_network_elec_s{simpl}_{clusters}.csv" - resources: mem_mb=10000 - script: "scripts/build_gas_network.py" + +if config["sector"]["gas_network"]: + + rule retrieve_gas_infrastructure_data: + output: "data/gas_network/scigrid-gas/data/IGGIELGN_LNGs.csv" + script: 'scripts/retrieve_gas_infrastructure_data.py' + + rule build_gas_network: + input: + gas_network="data/gas_network/gas_network_dataset.csv" + output: + cleaned_gas_network="resources/gas_network.csv" + resources: mem_mb=4000 + script: "scripts/build_gas_network.py" + + rule build_gas_import_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" + regions_onshore=pypsaeur("resources/regions_onshore_elec_s{simpl}_{clusters}.geojson"), + output: + gas_input_nodes="resources/gas_input_nodes_s{simpl}_{clusters}.csv" + resources: mem_mb=2000, + script: "scripts/build_gas_import_locations.py" + + rule cluster_gas_network: + input: + cleaned_gas_network="data/gas_network/gas_network_dataset.csv", + regions_onshore=pypsaeur("resources/regions_onshore_elec_s{simpl}_{clusters}.geojson"), + regions_offshore=pypsaeur("resources/regions_offshore_elec_s{simpl}_{clusters}.geojson") + output: + clustered_gas_network="resources/gas_network_elec_s{simpl}_{clusters}.csv" + resources: mem_mb=4000 + script: "scripts/cluster_gas_network.py" + + gas_infrastructure = {**rules.cluster_gas_network.output, **rules.build_gas_import_locations.output} +else: + gas_infrastructure = {} + rule build_heat_demands: input: @@ -354,7 +384,6 @@ rule prepare_sector_network: energy_totals_name='resources/energy_totals.csv', co2_totals_name='resources/co2_totals.csv', transport_name='resources/transport_data.csv', - clustered_gas_network="resources/gas_network_elec_s{simpl}_{clusters}.csv", traffic_data_KFZ="data/emobility/KFZ__count", traffic_data_Pkw="data/emobility/Pkw__count", biomass_potentials='resources/biomass_potentials_s{simpl}_{clusters}.csv', @@ -387,7 +416,8 @@ rule prepare_sector_network: solar_thermal_urban="resources/solar_thermal_urban_elec_s{simpl}_{clusters}.nc", solar_thermal_rural="resources/solar_thermal_rural_elec_s{simpl}_{clusters}.nc", **build_retro_cost_output, - **build_biomass_transport_costs_output + **build_biomass_transport_costs_output, + **gas_infrastructure output: RDIR + '/prenetworks/elec_s{simpl}_{clusters}_lv{lv}_{opts}_{sector_opts}_{planning_horizons}.nc' threads: 1 resources: mem_mb=2000 diff --git a/config.default.yaml b/config.default.yaml index 1254e926..61a70fac 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -457,7 +457,6 @@ plotting: gas for industry: '#853403' gas for industry CC: '#692e0a' gas pipeline: '#ebbca0' - Gas pipeline: '#ebbca0' # oil oil: '#c9c9c9' oil boiler: '#adadad' diff --git a/scripts/build_gas_import_locations.py b/scripts/build_gas_import_locations.py new file mode 100644 index 00000000..8ee6d577 --- /dev/null +++ b/scripts/build_gas_import_locations.py @@ -0,0 +1,76 @@ +""" +Build import locations for fossil gas from entry-points and LNG terminals. +""" + +import logging +logger = logging.getLogger(__name__) + +import pandas as pd +import geopandas as gpd + + +def read_scigrid_gas(fn): + df = gpd.read_file(fn) + df = pd.concat([df, df.param.apply(pd.Series)], axis=1) + df.drop(["param", "uncertainty", "method"], axis=1, inplace=True) + return df + + +def build_gas_input_locations(lng_fn, entry_fn, prod_fn): + + countries = snakemake.config["countries"] + countries[countries.index('GB')] = 'UK' + + # LNG terminals + lng = read_scigrid_gas(lng_fn) + + # Entry points from outside the model scope + entry = read_scigrid_gas(entry_fn) + entry["from_country"] = entry.from_country.str.rstrip() + entry = entry.loc[ + ~(entry.from_country.isin(countries) & entry.to_country.isin(countries)) & # only take non-EU entries + ~entry.name.str.contains("Tegelen") | # malformed datapoint + (entry.from_country == "NO") # entries from NO to GB + ] + + # production sites inside the model scope + prod = read_scigrid_gas(prod_fn) + prod = prod.loc[ + (prod.geometry.y > 35) & + (prod.geometry.x < 30) + ] + + return gpd.GeoDataFrame( + geometry=pd.concat([prod.geometry, entry.geometry, lng.geometry]).reset_index(drop=True), + crs=4326 + ) + + +if __name__ == "__main__": + + if 'snakemake' not in globals(): + from helper import mock_snakemake + snakemake = mock_snakemake( + 'build_gas_import_locations', + simpl='', + clusters='37', + ) + + logging.basicConfig(level=snakemake.config['logging_level']) + + onshore_regions = gpd.read_file(snakemake.input.regions_onshore).set_index('name') + + gas_input_locations = build_gas_input_locations( + snakemake.input.lng, + snakemake.input.entry, + snakemake.input.production + ) + + # recommended to use projected CRS rather than geographic CRS + gas_input_nodes = gpd.sjoin_nearest( + gas_input_locations.to_crs(3035), + onshore_regions.to_crs(3035), + how='left' + ).index_right.unique() + + pd.Series(gas_input_nodes, name='gas_input_nodes').to_csv(snakemake.output.gas_input_nodes) \ No newline at end of file diff --git a/scripts/build_gas_network.py b/scripts/build_gas_network.py old mode 100755 new mode 100644 index 0f213416..cf4a86a1 --- a/scripts/build_gas_network.py +++ b/scripts/build_gas_network.py @@ -1,5 +1,5 @@ """ -Builds clustered natural gas network based on data from: +Preprocess gas network based on data from: [1] the SciGRID Gas project (https://www.gas.scigrid.de/) @@ -15,174 +15,25 @@ import re import json import pandas as pd -import geopandas as gpd -import numpy as np - from shapely.geometry import Point +from pypsa.geo import haversine_pts -def concat_gdf(gdf_list, crs='EPSG:4326'): - """Convert to gepandas dataframe with given Coordinate Reference System (crs).""" - return gpd.GeoDataFrame(pd.concat(gdf_list),crs=crs) - - -def string2list(string, with_None=True): +def string2list(string, with_none=True): """Convert string format to a list.""" - p = re.compile('(? MW - clean_pipes.loc[:, 'capacity_recalculated'] *= 1e3 - - # rename columns - to_rename = { - 'capacity_recalculated': 'pipe_capacity_MW', - 'buses_start': 'bus0', - 'buses_destination': 'bus1' - } - clean_pipes.rename(columns=to_rename, inplace=True) - - return clean_pipes - - def diameter2capacity(pipe_diameter_mm): - """Calculate pipe capacity based on diameter. + """Calculate pipe capacity in MW based on diameter in mm. 20 inch (500 mm) 50 bar -> 1.5 GW CH4 pipe capacity (LHV) 24 inch (600 mm) 50 bar -> 5 GW CH4 pipe capacity (LHV) @@ -193,69 +44,107 @@ def diameter2capacity(pipe_diameter_mm): """ # slopes definitions - m0 = (5 - 1.5) / (600 - 500) - m1 = (11.25 - 5) / (900 - 600) - m2 = (21.7 - 11.25) / (1200 - 900) + m0 = (1500 - 0) / (500 - 0) + m1 = (5000 - 1500) / (600 - 500) + m2 = (11250 - 5000) / (900 - 600) + m3 = (21700 - 11250) / (1200 - 900) # intercept - a0 = -16 - a1 = -7.5 - a2 = -20.1 + a0 = 0 + a1 = -16000 + a2 = -7500 + a3 = -20100 if pipe_diameter_mm < 500: - return np.nan - elif pipe_diameter_mm < 600: return a0 + m0 * pipe_diameter_mm - elif pipe_diameter_mm < 900: + elif pipe_diameter_mm < 600: return a1 + m1 * pipe_diameter_mm - else: + elif pipe_diameter_mm < 900: return a2 + m2 * pipe_diameter_mm + else: + return a3 + m3 * pipe_diameter_mm -def determine_pipe_capacity(gas_network): - """Check pipe capacity depending on diameter and pressure.""" +def find_terminal_points(df): + + latlon = [] - gas_network["capacity_recalculated"] = gas_network.diameter_mm.apply(diameter2capacity) + for attr in ["lat", "long"]: - # if pipe capacity smaller than 1.5 GW take original pipe capacity - low_cap = gas_network.Capacity_GWh_h < 1.5 - gas_network.loc[low_cap, "capacity_recalculated"] = gas_network.loc[low_cap, "capacity_recalculated"].fillna(gas_network.loc[low_cap, "Capacity_GWh_h"]) + s = df[attr].apply(string2list) + + s = s.apply(lambda x: [x[0], x[-1]]) + + latlon.append(pd.DataFrame(s.to_list(), + columns=[f"{attr}0", f"{attr}1"] + )) - # for pipes without diameter assume 500 mm diameter - gas_network["capacity_recalculated"].fillna(1.5, inplace=True) + latlon = pd.concat(latlon, axis=1) - # for nord stream take orginal data - nord_stream = gas_network[gas_network.max_pressure_bar==220].index - gas_network.loc[nord_stream, "capacity_recalculated"] = gas_network.loc[nord_stream, "Capacity_GWh_h"] + points = latlon.apply( + lambda x: { + "point0": Point(x.long0, x.lat0), + "point1": Point(x.long1, x.lat1) + }, + axis=1, + result_type='expand' + ) + + return pd.concat([df, points], axis=1) + + +def process_gas_network_data(fn): + + df = pd.read_csv(fn, sep=',') + + df = find_terminal_points(df) + + to_drop = ["name", "source_id", "country_code", "node_id", + "long", "lat", "lat_mean", "long_mean", "num_compressor"] + df.drop(to_drop, axis=1, inplace=True) + + to_rename = { + "is_bothDirection": "bidirectional", + "start_year": "build_year", + "length_km": "length", + "Capacity_GWh_h": "p_nom_data", + "id": "tags", + } + df.rename(columns=to_rename, inplace=True) + + df.bidirectional = df.bidirectional.astype(bool) + + # convert from GWh/h to MW + df.p_nom_data *= 1e3 + + # for pipes with missing diameter, assume 500 mm + df.loc[df.diameter_mm.isna(), "diameter_mm"] = 500. + + # for nord stream and small pipelines take original capacity data + # otherwise inferred values from pipe diameter + df["p_nom"] = df.diameter_mm.map(diameter2capacity) + df.p_nom.update( + df.p_nom_data.where((df.diameter_mm < 500) | (df.max_pressure_bar == 220)) + ) + + df["length_haversine"] = df.apply( + lambda p: 1.5 * haversine_pts([p.point0.x, p.point1.y], [p.point1.x, p.point1.y]), + axis=1 + ) + + df.length.update(df.length_haversine.where(df.length.isna())) + + return df if __name__ == "__main__": if 'snakemake' not in globals(): from helper import mock_snakemake - snakemake = mock_snakemake('build_gas_network', - network='elec', simpl='', clusters='37', - lv='1.0', opts='', planning_horizons='2020', - sector_opts='168H-T-H-B-I') + snakemake = mock_snakemake('build_gas_network') logging.basicConfig(level=snakemake.config['logging_level']) - # import gas network data - gas_network, points = load_gas_network(snakemake.input.gas_network) + gas_network = process_gas_network_data(snakemake.input.gas_network) - # get clustered bus regions - bus_regions = load_bus_regions( - snakemake.input.regions_onshore, - snakemake.input.regions_offshore - ) - nodes = pd.Index(bus_regions.name.unique()) - - # map gas network points to network buses - points2buses_map = points2buses(points, bus_regions) - - # create gas network between pypsa nodes - gas_connections = build_gas_network_topology(gas_network, points2buses_map) - - gas_connections = clean_dataset(nodes, gas_connections) - - gas_connections.to_csv(snakemake.output.clustered_gas_network) + gas_network.to_csv(snakemake.output.cleaned_gas_network) \ No newline at end of file diff --git a/scripts/cluster_gas_network.py b/scripts/cluster_gas_network.py new file mode 100755 index 00000000..00dd165d --- /dev/null +++ b/scripts/cluster_gas_network.py @@ -0,0 +1,110 @@ +"""Cluster gas network.""" + +import logging +logger = logging.getLogger(__name__) + +import pandas as pd +import geopandas as gpd + +from shapely import wkt + + +def concat_gdf(gdf_list, crs='EPSG:4326'): + """Concatenate multiple geopandas dataframes with common coordinate reference system (crs).""" + return gpd.GeoDataFrame(pd.concat(gdf_list), crs=crs) + + +def load_bus_regions(onshore_path, offshore_path): + """Load pypsa-eur on- and offshore regions and concat.""" + + bus_regions_offshore = gpd.read_file(offshore_path) + bus_regions_onshore = gpd.read_file(onshore_path) + bus_regions = concat_gdf([bus_regions_offshore, bus_regions_onshore]) + bus_regions = bus_regions.dissolve(by='name', aggfunc='sum') + + return bus_regions + + +def build_clustered_gas_network(df, bus_regions): + + for i in [0,1]: + + gdf = gpd.GeoDataFrame(geometry=df[f"point{i}"], crs="EPSG:4326") + + bus_mapping = gpd.sjoin(gdf, bus_regions, how="left", op="within").index_right + bus_mapping = bus_mapping.groupby(bus_mapping.index).first() + + df[f"bus{i}"] = bus_mapping + + df.drop(["point0", "point1"], axis=1, inplace=True) + + # drop pipes where not both buses are inside regions + df = df.loc[~df.bus0.isna() & ~df.bus1.isna()] + + # drop pipes within one region + df = df.loc[df.bus1 != df.bus0] + + # create new numbered index + df.reset_index(drop=True, inplace=True) + + return df + + +def reindex_pipes(df): + + def make_index(x): + connector = " <-> " if x.bidirectional else " -> " + return "gas pipeline " + x.bus0 + connector + x.bus1 + + df.index = df.apply(make_index, axis=1) + + df["p_min_pu"] = df.bidirectional.apply(lambda bi: -1 if bi else 0) + df.drop("bidirectional", axis=1, inplace=True) + + df.sort_index(axis=1, inplace=True) + + +def aggregate_parallel_pipes(df): + + strategies = { + 'bus0': 'first', + 'bus1': 'first', + "p_nom": 'sum', + "p_nom_data": 'sum', + "max_pressure_bar": "mean", + "build_year": "mean", + "diameter_mm": "mean", + "length": 'mean', + 'tags': ' '.join, + } + df = df.groupby(df.index).agg(strategies) + + +if __name__ == "__main__": + + if 'snakemake' not in globals(): + from helper import mock_snakemake + snakemake = mock_snakemake( + 'cluster_gas_network', + simpl='', + clusters='37' + ) + + logging.basicConfig(level=snakemake.config['logging_level']) + + fn = snakemake.input.cleaned_gas_network + df = pd.read_csv(fn, index_col=0) + for col in ["point0", "point1"]: + df[col] = df[col].apply(wkt.loads) + + bus_regions = load_bus_regions( + snakemake.input.regions_onshore, + snakemake.input.regions_offshore + ) + + gas_network = build_clustered_gas_network(df, bus_regions) + + reindex_pipes(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/helper.py b/scripts/helper.py index 21f4e8c0..b176ccee 100644 --- a/scripts/helper.py +++ b/scripts/helper.py @@ -89,3 +89,15 @@ def mock_snakemake(rulename, **wildcards): os.chdir(script_dir) return snakemake + +# from pypsa-eur/_helpers.py +def progress_retrieve(url, file): + import urllib + from progressbar import ProgressBar + + pbar = ProgressBar(0, 100) + + def dlProgress(count, blockSize, totalSize): + pbar.update( int(count * blockSize * 100 / totalSize) ) + + urllib.request.urlretrieve(url, file, reporthook=dlProgress) \ No newline at end of file diff --git a/scripts/plot_network.py b/scripts/plot_network.py index a78b6551..2a6073d7 100644 --- a/scripts/plot_network.py +++ b/scripts/plot_network.py @@ -332,7 +332,7 @@ def plot_ch4_map(network): supply_energy = get_nodal_balance().droplevel([0,1]).sort_index() - if "Gas pipeline" not in n.links.carrier.unique(): + if "gas pipeline" not in n.links.carrier.unique(): return assign_location(n) @@ -372,7 +372,7 @@ def plot_ch4_map(network): 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) + n.links.drop(n.links.index[n.links.carrier != "gas pipeline"], inplace=True) link_widths = n.links.p_nom_opt / linewidth_factor link_widths[n.links.p_nom_opt < line_lower_threshold] = 0. @@ -426,7 +426,7 @@ def plot_ch4_map(network): bbox_inches="tight") ################################################## - supply_energy.drop("Gas pipeline", level=1, inplace=True) + 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) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 9cf9acf0..df8bc081 100644 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -1033,8 +1033,8 @@ def add_electricity_grid_connection(n, costs): n.generators.loc[gens, "capital_cost"] += costs.at['electricity grid connection', 'fixed'] -def add_storage(n, costs): - print("adding electricity and hydrogen storage") +def add_storage_and_grids(n, costs): + print("adding electricity and hydrogen storage as well as hydrogen and gas grids") nodes = pop_layout.index @@ -1106,11 +1106,93 @@ def add_storage(n, costs): capital_cost=h2_capital_cost ) + if options["gas_network"]: + + 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) + + if options["H2_retrofit"]: + gas_pipes["p_nom_max"] = gas_pipes.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["capital_cost"] = gas_pipes.length * costs.at['CH4 (g) pipeline', 'fixed'] + + n.madd("Link", + gas_pipes.index, + bus0=gas_pipes.bus0 + " gas", + bus1=gas_pipes.bus1 + " gas", + p_min_pu=gas_pipes.p_min_pu, + p_nom=gas_pipes.p_nom, + p_nom_extendable=True, + p_nom_max=gas_pipes.p_nom_max, + p_nom_min=gas_pipes.p_nom_min, + length=gas_pipes.length, + capital_cost=gas_pipes.capital_cost, + tags=gas_pipes.tags, + carrier="gas pipeline", + lifetime=50 + ) + + # 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[ + (n.generators.carrier=="gas") & + ~n.generators.bus.map(n.buses.location).isin(gas_input_nodes) + ].index + n.generators.drop(remove_i, inplace=True) + + # TODO candidate gas network topology + + + # retroftting existing CH4 pipes to H2 pipes + if options["gas_network"] and options["H2_retrofit"]: + + gas_pipe_i = n.links[n.links.carrier == "gas pipeline"].index + n.links.loc[gas_pipe_i, "p_nom_extendable"] = True + h2_pipes = gas_pipes.rename(index=lambda x: + x.replace("gas pipeline", "H2 pipeline retrofitted")) + + n.madd("Link", + 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_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, + carrier="H2 pipeline retrofitted", + lifetime=50 + ) + attrs = ["bus0", "bus1", "length"] h2_links = pd.DataFrame(columns=attrs) - candidates = pd.concat({"lines": n.lines[attrs], - "links": n.links.loc[n.links.carrier == "DC", attrs]}) + lines_sel = n.lines[attrs] + links_sel = n.links.loc[n.links.carrier.isin(["DC", "gas pipeline"]), attrs] + + candidates = pd.concat({ + "lines": lines_sel, + "links": links_sel, + }) for candidate in candidates.index: buses = [candidates.at[candidate, "bus0"], candidates.at[candidate, "bus1"]] @@ -1134,96 +1216,6 @@ def add_storage(n, costs): lifetime=costs.at['H2 (g) pipeline', 'lifetime'] ) - if options["gas_network"]: - - logger.info("Add gas network") - - cols = [ - "bus0", - "bus1", - "is_bothDirection", - "pipe_capacity_MW", - "id", - "length_km" - ] - gas_pipes = pd.read_csv(snakemake.input.clustered_gas_network, usecols=cols) - - def make_index(x): - connector = " <-> " if x.is_bothDirection else " -> " - return "Gas pipeline " + x.bus0 + connector + x.bus1 - - gas_pipes.index = gas_pipes.apply(make_index, axis=1) - - # group parallel pipes together - strategies = { - 'bus0': 'first', - 'bus1': 'first', - 'is_bothDirection': 'first', - "pipe_capacity_MW": 'sum', - "length_km": 'sum', - 'id': ' '.join, - } - gas_pipes = gas_pipes.groupby(gas_pipes.index).agg(strategies) - - gas_pipes["num_parallel"] = gas_pipes.index.value_counts() - gas_pipes["p_min_pu"] = gas_pipes.apply(lambda x: -1 if x.is_bothDirection else 0, axis=1) - - if options["H2_retrofit"]: - gas_pipes["p_nom_max"] = gas_pipes.gas_pipes.pipe_capacity_MW - 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.pipe_capacity_MW - gas_pipes["capital_cost"] = gas_pipes.length_km * costs.at['CH4 (g) pipeline', 'fixed'] - - n.madd("Link", - gas_pipes.index, - bus0=gas_pipes.bus0 + " gas", - bus1=gas_pipes.bus1 + " gas", - p_min_pu=gas_pipes.p_min_pu, - p_nom=gas_pipes.pipe_capacity_MW, - p_nom_extendable=True, - p_nom_max=gas_pipes.p_nom_max, - p_nom_min=gas_pipes.p_nom_min, - length=gas_pipes.length_km, - capital_cost=gas_pipes.capital_cost, - type=gas_pipes.num_parallel, - tags=gas_pipes.id, - carrier="Gas pipeline", - lifetime=50 - ) - - # remove fossil generators at all connected nodes - # TODO what should be assumed here? rather located at LNG terminals? - missing = nodes.difference(pd.concat([gas_pipes.bus0, gas_pipes.bus1]).unique()) - remove_i = n.generators[(n.generators.carrier=="gas") - & (~n.generators.bus.str.replace(" gas","").isin(missing))].index - n.generators.drop(remove_i, inplace=True) - - # retroftting existing CH4 pipes to H2 pipes - if options["gas_network"] and options["H2_retrofit"]: - - gas_pipe_i = n.links[n.links.carrier == "Gas pipeline"].index - n.links.loc[gas_pipe_i, "p_nom_extendable"] = True - h2_pipes = gas_pipes.rename(index=lambda x: - x.replace("Gas pipeline", "H2 pipeline retrofitted")) - - n.madd("Link", - 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_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, - carrier="H2 pipeline retrofitted", - lifetime=50 - ) - n.add("Carrier", "battery") n.madd("Bus", @@ -1963,14 +1955,6 @@ def add_industry(n, costs): # 1e6 to convert TWh to MWh industrial_demand = pd.read_csv(snakemake.input.industrial_demand, index_col=0) * 1e6 - methane_demand = industrial_demand.loc[nodes, "methane"].div(8760).rename(index=lambda x: x + " gas for industry") - - # need to aggregate methane demand if gas not nodally resolved - if not options["gas_network"]: - methane_demand = methane_demand.sum() - - solid_biomass_by_country = industrial_demand["solid biomass"].groupby(pop_layout.ct).sum() - n.madd("Bus", spatial.biomass.industry, location=spatial.biomass.locations, @@ -2018,11 +2002,18 @@ def add_industry(n, costs): location=spatial.gas.locations, carrier="gas for industry") + gas_demand = industrial_demand.loc[nodes, "methane"] / 8760. + + if options["gas_network"]: + spatial_gas_demand = gas_demand.rename(index=lambda x: x + " gas for industry") + else: + spatial_gas_demand = gas_demand.sum() + n.madd("Load", spatial.gas.industry, bus=spatial.gas.industry, carrier="gas for industry", - p_set=methane_demand + p_set=spatial_gas_demand ) n.madd("Link", @@ -2463,7 +2454,7 @@ if __name__ == "__main__": add_generation(n, costs) - add_storage(n, costs) + add_storage_and_grids(n, costs) # TODO merge with opts cost adjustment below for o in opts: diff --git a/scripts/retrieve_gas_infrastructure_data.py b/scripts/retrieve_gas_infrastructure_data.py new file mode 100644 index 00000000..4bae3e29 --- /dev/null +++ b/scripts/retrieve_gas_infrastructure_data.py @@ -0,0 +1,36 @@ +""" +Retrieve gas infrastructure data from https://zenodo.org/record/4767098/files/IGGIELGN.zip +""" + +import logging +from helper import progress_retrieve + +import zipfile +from pathlib import Path + +logger = logging.getLogger(__name__) + + +if __name__ == "__main__": + if 'snakemake' not in globals(): + from helper import mock_snakemake + snakemake = mock_snakemake('retrieve_gas_network_data') + rootpath = '..' + else: + rootpath = '.' + + url = "https://zenodo.org/record/4767098/files/IGGIELGN.zip" + + # Save locations + zip_fn = Path(f"{rootpath}/IGGIELGN.zip") + to_fn = Path(f"{rootpath}/data/gas_network/scigrid-gas") + + logger.info(f"Downloading databundle from '{url}'.") + progress_retrieve(url, zip_fn) + + logger.info(f"Extracting databundle.") + zipfile.ZipFile(zip_fn).extractall(to_fn) + + zip_fn.unlink() + + logger.info(f"Gas infrastructure data available in '{to_fn}'.") diff --git a/scripts/solve_network.py b/scripts/solve_network.py index 25666caf..c784c913 100644 --- a/scripts/solve_network.py +++ b/scripts/solve_network.py @@ -189,7 +189,7 @@ def add_chp_constraints(n): def add_pipe_retrofit_constraint(n): """Add constraint for retrofitting existing CH4 pipelines to H2 pipelines.""" - gas_pipes_i = n.links[n.links.carrier=="Gas pipeline"].index + gas_pipes_i = n.links[n.links.carrier=="gas pipeline"].index h2_retrofitted_i = n.links[n.links.carrier=='H2 pipeline retrofitted'].index if h2_retrofitted_i.empty or gas_pipes_i.empty: return @@ -201,7 +201,7 @@ def add_pipe_retrofit_constraint(n): CH4_per_H2 = 1 / n.config["sector"]["H2_retrofit_capacity_per_CH4"] lhs = linexpr( - (CH4_per_H2, link_p_nom.loc[h2_retrofitted_i].rename(index=lambda x: x.replace("H2 pipeline retrofitted", "Gas pipeline"))), + (CH4_per_H2, link_p_nom.loc[h2_retrofitted_i].rename(index=lambda x: x.replace("H2 pipeline retrofitted", "gas pipeline"))), (1, link_p_nom.loc[gas_pipes_i]) )