revision gas infrastructure representation

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"
 
+
+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",
+            gas_network="data/gas_network/gas_network_dataset.csv"
-        country_shapes=pypsaeur("resources/country_shapes.geojson"),
+        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=10000
+        resources: mem_mb=4000
-    script: "scripts/build_gas_network.py"
+        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

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'

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)

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'):
+def string2list(string, with_none=True):
-    """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):
     """Convert string format to a list."""
-    p = re.compile('(?<!\\\\)\'')
-    string = p.sub('\"', string)
 
-    if with_None:
+    if with_none:
         p2 = re.compile('None')
         string = p2.sub('\"None\"', string)
+    else:
+        p = re.compile('(?<!\\\\)\'')
+        string = p.sub('\"', string)
 
     return json.loads(string)
 
 
-def load_gas_network(df_path):
-    """Load and format gas network data."""
-
-    df = pd.read_csv(df_path, sep=',')
-
-    df.long = df.long.apply(string2list)
-    df.lat = df.lat.apply(string2list)
-    df.node_id = df.node_id.apply(string2list)
-
-    # pipes which can be used in both directions
-    both_direct_df = df[df.is_bothDirection == 1].reset_index(drop=True)
-    both_direct_df.node_id = both_direct_df.node_id.apply(lambda x: [x[1], x[0]])
-    both_direct_df.long = both_direct_df.long.apply(lambda x: [x[1], x[0]])
-    both_direct_df.lat = both_direct_df.lat.apply(lambda x: [x[1], x[0]])
-
-    df_singledirect = pd.concat([df, both_direct_df]).reset_index(drop=True)
-    df_singledirect.drop('is_bothDirection', axis=1)
-
-    # create shapely geometry points
-    df['point1'] = df.apply(lambda x: Point((x['long'][0], x['lat'][0])), axis=1)
-    df['point2'] = df.apply(lambda x: Point((x['long'][1], x['lat'][1])), axis=1)
-    df['point1_name'] = df.node_id.str[0]
-    df['point2_name'] = df.node_id.str[1]
-
-    part1 = df[['point1', 'point1_name']]
-    part2 = df[['point2', 'point2_name']]
-    part1.columns = ['geometry', 'name']
-    part2.columns = ['geometry', 'name']
-    points = [part1, part2]
-    points = concat_gdf(points)
-    points = points.drop_duplicates()
-    points.reset_index(drop=True, inplace=True)
-
-    return df, points
-
-
-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')
-    bus_regions = bus_regions.reset_index()
-
-    return bus_regions
-
-
-def points2buses(input_points, bus_regions):
-    """Map gas network points to network buses depending on bus region."""
-
-    points = input_points.copy()
-    points['bus'] = None
-    buses_list = set(bus_regions.name)
-    for bus in buses_list:
-        mask = bus_regions[bus_regions.name == bus]
-        index = gpd.clip(points, mask).index
-        points.loc[index, 'bus'] = bus
-
-    return points
-
-
-def build_gas_network_topology(df, points2buses):
-    """Create gas network between pypsa buses.
-
-    Parameters
-    ----------
-    df : pd.DataFrame
-        gas network data
-    points2buses_map : pd.DataFrame
-        mapping of gas network points to pypsa buses
-
-    Returns
-    -------
-    gas_connections : pd.DataFrame
-        gas network connecting pypsa buses
-    """
-
-    tmp_df = points2buses[['bus', 'name']]
-
-    tmp_df.columns = ['buses_start', 'name']
-    gas_connections = df.merge(tmp_df, left_on='point1_name', right_on='name')
-    
-    tmp_df.columns = ['buses_destination', 'name']
-    gas_connections = gas_connections.merge(tmp_df, left_on='point2_name', right_on='name')
-    
-    # drop all pipes connecting the same bus
-    gas_connections = gas_connections[gas_connections.buses_start != gas_connections.buses_destination]
-    gas_connections.reset_index(drop=True, inplace=True)
-    gas_connections.drop(['point1', 'point2'], axis=1, inplace=True)
-
-    return gas_connections
-
-
-def check_missing(nodes, gas_connections):
-    """Check which nodes are not connected to the gas network."""
-
-    start_buses = gas_connections.buses_start.dropna().unique()
-    end_buses = gas_connections.buses_destination.dropna().unique()
-
-    missing_start = nodes[[bus not in start_buses for bus in nodes]]
-    missing_end = nodes[[bus not in end_buses for bus in nodes]]
-
-    logger.info(f"- The following buses are missing in gas network data as a start bus:"
-                f"\n {', '.join(map(str, missing_start))} \n"
-                f"- The following buses are missing in gas network data as an end bus:"
-                f"\n {', '.join(map(str, missing_end))} \n"
-                f"- The following buses are missing completely:"
-                f"\n {', '.join(map(str, missing_start.intersection(missing_end)))}")
-
-
-def clean_dataset(nodes, gas_connections):
-    """Convert units and save only necessary data."""
-
-    check_missing(nodes, gas_connections)
-
-    determine_pipe_capacity(gas_connections)
-    
-    cols = [
-        'is_bothDirection',
-        'capacity_recalculated',
-        'buses_start',
-        'buses_destination',
-        'id',
-        'length_km'
-    ]
-    clean_pipes = gas_connections[cols].dropna()
-
-
-    # convert GW -> 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)
+    m0 = (1500 - 0) / (500 - 0)
-    m1 = (11.25 - 5) / (900 - 600)
+    m1 = (5000 - 1500) / (600 - 500)
-    m2 = (21.7 - 11.25) / (1200 - 900)
+    m2 = (11250 - 5000) / (900 - 600)
+    m3 = (21700 - 11250) / (1200 - 900)
 
     # intercept
-    a0 = -16
+    a0 = 0
-    a1 = -7.5
+    a1 = -16000
-    a2 = -20.1
+    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):
+def find_terminal_points(df):
-    """Check pipe capacity depending on diameter and pressure."""
     
-    gas_network["capacity_recalculated"] = gas_network.diameter_mm.apply(diameter2capacity)
+    latlon = []
 
-    # if pipe capacity smaller than 1.5 GW take original pipe capacity
+    for attr in ["lat", "long"]:
-    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"])
     
-    # for pipes without diameter assume 500 mm diameter
+        s = df[attr].apply(string2list)
-    gas_network["capacity_recalculated"].fillna(1.5, inplace=True)
 
-    # for nord stream take orginal data
+        s = s.apply(lambda x: [x[0], x[-1]])
-    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"]
+        latlon.append(pd.DataFrame(s.to_list(),
+            columns=[f"{attr}0", f"{attr}1"]
+        ))
+    
+    latlon = pd.concat(latlon, axis=1)
+    
+    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',
+        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')
 
     logging.basicConfig(level=snakemake.config['logging_level'])
 
-    # import gas network data
+    gas_network = process_gas_network_data(snakemake.input.gas_network)
-    gas_network, points = load_gas_network(snakemake.input.gas_network)
 
-    # get clustered bus regions
+    gas_network.to_csv(snakemake.output.cleaned_gas_network)
-    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)

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)

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)

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)
 

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):
+def add_storage_and_grids(n, costs):
-    print("adding electricity and hydrogen storage")
+    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],
+    lines_sel = n.lines[attrs]
-                            "links": n.links.loc[n.links.carrier == "DC", 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:

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}'.")

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