diff --git a/scripts/build_gas_network.py b/scripts/build_gas_network.py index a3690dbb..0f213416 100755 --- a/scripts/build_gas_network.py +++ b/scripts/build_gas_network.py @@ -1,56 +1,27 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Builds clustered natural gas network based on data from: + [1] the SciGRID Gas project (https://www.gas.scigrid.de/) + [2] ENTSOG capacity map (https://www.entsog.eu/sites/default/files/2019-10/Capacities%20for%20Transmission%20Capacity%20Map%20RTS008_NS%20-%20DWH_final.xlsx) - -Relevant Settings ------------------ -.. code:: yaml - sector: - - -Inputs ------- -gas network data from SciGRID gas and ENTSOG: - - gas_network="data/gas_network/gas_network_dataset.csv", - combined gas network data set from [1] and [2] - - 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") - -Outputs -------- -clustered gas network data for corresponding PyPSA-Eur-Sec network - - gas_network='resources/gas_network_{clusters}.csv' - -Description ------------ - """ -import geoplot -import geoplot.crs as gcrs -import matplotlib.pyplot as plt import logging logger = logging.getLogger(__name__) import re -import pandas as pd import json -from shapely.geometry import LineString,Point - +import pandas as pd import geopandas as gpd import numpy as np -#-----------------########################################################## -# helper functions # -#-----------------# -def concat_gdf(gdf_list, crs = 'EPSG:4326'): +from shapely.geometry import Point + + +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) @@ -67,11 +38,9 @@ def string2list(string, with_None=True): return json.loads(string) -#-----------------############################################################ -# main functions # -#-----------------# -def preprocessing(df_path): - """Load and format gas network data.""" +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) @@ -93,7 +62,6 @@ def preprocessing(df_path): 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'] @@ -106,90 +74,114 @@ def preprocessing(df_path): return df, points -def load_region(onshore_path, offshore_path): +def load_bus_regions(onshore_path, offshore_path): """Load pypsa-eur on- and offshore regions and concat.""" - buses_region_offshore = gpd.read_file(offshore_path) - buses_region_onshore = gpd.read_file(onshore_path) - buses_region = concat_gdf([buses_region_offshore, buses_region_onshore]) - buses_region = buses_region.dissolve(by='name', aggfunc='sum') - buses_region = buses_region.reset_index() - return buses_region + 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 create_points2buses_map(input_points, buses_region): +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(buses_region.name) + buses_list = set(bus_regions.name) for bus in buses_list: - mask = buses_region[buses_region.name == bus] + mask = bus_regions[bus_regions.name == bus] index = gpd.clip(points, mask).index points.loc[index, 'bus'] = bus return points -def create_cross_regions_network(df, points2buses_map): +def build_gas_network_topology(df, points2buses): """Create gas network between pypsa buses. - Input: - df : gas network data (pd.DataFrame) - points2buses_map : map gas network points to pypsa buses (pd.DataFrame) + Parameters + ---------- + df : pd.DataFrame + gas network data + points2buses_map : pd.DataFrame + mapping of gas network points to pypsa buses - Return: - cross_buses_gas_network : gas network connecting pypsa buses - (pd.DataFrame) + Returns + ------- + gas_connections : pd.DataFrame + gas network connecting pypsa buses """ - tmp_df = points2buses_map[['bus', 'name']] - tmp_df.columns = ['buses_start','name'] - cross_buses_gas_network = df.merge(tmp_df, left_on='point1_name', - right_on='name') + + 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'] - cross_buses_gas_network = cross_buses_gas_network.merge(tmp_df, - left_on='point2_name', - right_on='name') + gas_connections = gas_connections.merge(tmp_df, left_on='point2_name', right_on='name') + # drop all pipes connecting the same bus - cross_buses_gas_network = cross_buses_gas_network[cross_buses_gas_network.buses_start \ - != cross_buses_gas_network.buses_destination] - cross_buses_gas_network.reset_index(drop=True, inplace=True) - cross_buses_gas_network.drop(['point1','point2'], axis=1, inplace=True) + 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 cross_buses_gas_network + return gas_connections -def check_missing(nodes, cross_buses_gas_network): +def check_missing(nodes, gas_connections): """Check which nodes are not connected to the gas network.""" - missing0 = nodes[[bus not in cross_buses_gas_network.buses_start.dropna().unique() - for bus in nodes]] - missing1 = nodes[[bus not in cross_buses_gas_network.buses_destination.dropna().unique() - for bus in nodes]] - logger.info("\n - The following buses are missing in gas network data as a start bus: \n {} \n" - "- The following buses are missing in gas network data as an end bus: \n {} \n " - "- The following buses are missing completely: \n {}" - .format(', '.join(map(str, missing0)), - ', '.join(map(str, missing1)), - ', '.join(map(str, missing0.intersection(missing1))))) + + 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(cross_buses_gas_network): +def clean_dataset(nodes, gas_connections): """Convert units and save only necessary data.""" - inspect_pipe_capacity(cross_buses_gas_network) - cols = ['is_bothDirection', 'capacity_recalculated','buses_start', - 'buses_destination', 'id', 'length_km'] - clean_pipes = cross_buses_gas_network[cols].dropna() + + 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 - clean_pipes.rename(columns={'capacity_recalculated': 'pipe_capacity_MW', - 'buses_start': 'bus0', - 'buses_destination': 'bus1'}, inplace=True) + 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 recalculate_pipe_capacity(pipe_diameter_mm): +def diameter2capacity(pipe_diameter_mm): """Calculate pipe capacity based on diameter. 20 inch (500 mm) 50 bar -> 1.5 GW CH4 pipe capacity (LHV) @@ -197,97 +189,73 @@ def recalculate_pipe_capacity(pipe_diameter_mm): 36 inch (900 mm) 50 bar -> 11.25 GW CH4 pipe capacity (LHV) 48 inch (1200 mm) 80 bar -> 21.7 GW CH4 pipe capacity (LHV) - Based on p.15 of (https://gasforclimate2050.eu/wp-content/uploads/2020/07/2020_European-Hydrogen-Backbone_Report.pdf""" - # slope - m0 = (5-1.5) / (600-500) - m1 = (11.25-5)/(900-600) - m2 = (21.7-11.25)/(1200-900) + Based on p.15 of https://gasforclimate2050.eu/wp-content/uploads/2020/07/2020_European-Hydrogen-Backbone_Report.pdf + """ - if pipe_diameter_mm<500: + # slopes definitions + m0 = (5 - 1.5) / (600 - 500) + m1 = (11.25 - 5) / (900 - 600) + m2 = (21.7 - 11.25) / (1200 - 900) + + # intercept + a0 = -16 + a1 = -7.5 + a2 = -20.1 + + if pipe_diameter_mm < 500: return np.nan - if pipe_diameter_mm<600 and pipe_diameter_mm>=500: - return -16 + m0 * pipe_diameter_mm - if pipe_diameter_mm<900 and pipe_diameter_mm>=600: - return -7.5 + m1 * pipe_diameter_mm + elif pipe_diameter_mm < 600: + return a0 + m0 * pipe_diameter_mm + elif pipe_diameter_mm < 900: + return a1 + m1 * pipe_diameter_mm else: - return -20.1 + m2 * pipe_diameter_mm + return a2 + m2 * pipe_diameter_mm -def inspect_pipe_capacity(gas_network): + +def determine_pipe_capacity(gas_network): """Check pipe capacity depending on diameter and pressure.""" - gas_network["capacity_recalculated"] = gas_network.diameter_mm.apply(recalculate_pipe_capacity) - low_cap = gas_network.Capacity_GWh_h < 1.5 + + gas_network["capacity_recalculated"] = gas_network.diameter_mm.apply(diameter2capacity) + # if pipe capacity smaller than 1.5 GW take original pipe capacity - gas_network.loc[low_cap, "capacity_recalculated"] = gas_network.loc[low_cap, "capacity_recalculated"].fillna(gas_network.loc[low_cap,"Capacity_GWh_h"]) + 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 gas_network["capacity_recalculated"].fillna(1.5, inplace=True) - # nord stream take orginal data + + # 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"] -# ----------- VISULAISATION -------------------------------------------------- -def create_view_object(cbgn_no_duplicate,buses_region): - """Create object to view gas network data.""" - cbgn_no_duplicate=cbgn_no_duplicate.merge(buses_region,left_on='buses_start',right_on='name') - cbgn_no_duplicate=cbgn_no_duplicate.merge(buses_region,left_on='buses_destination',right_on='name') - cbgn_no_duplicate.geometry_x=cbgn_no_duplicate.geometry_x.apply(lambda x: x.centroid) - cbgn_no_duplicate.geometry_y=cbgn_no_duplicate.geometry_y.apply(lambda x: x.centroid) - cbgn_no_duplicate['geometry']=list(zip(cbgn_no_duplicate['geometry_x'], - cbgn_no_duplicate['geometry_y'])) - - final = cbgn_no_duplicate[['buses_start', 'buses_destination', - 'Capacity_GWh_h', 'geometry']] - final['geometry'] = final['geometry'].apply(LineString) - final=gpd.GeoDataFrame(final,crs='EPSG:4326') - - return final - - -def view(cbgn_no_duplicate, buses_region, shapes_path): - """Plot gas network.""" - final = create_view_object(cbgn_no_duplicate,buses_region) - - eu=gpd.read_file(shapes_path) - ax = geoplot.webmap(eu, projection=gcrs.WebMercator(), figsize=(20,20), - alpha=0.5) - geoplot.choropleth(buses_region, hue='name',ax=ax, alpha=0.2, - edgecolor='red', linewidth=2) - geoplot.sankey( final, scale='Capacity_GWh_h', hue='Capacity_GWh_h', - cmap='viridis', ax=ax, legend=True, legend_var='hue') - plt.savefig("../graphics/clustered-gas-network_{}.pdf".format(snakemake.wildcards.clusters), - bbox_inches='tight', pad_inches=0.1) - -#%% if __name__ == "__main__": - # for testing 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') + 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, points = preprocessing(snakemake.input.gas_network) + gas_network, points = load_gas_network(snakemake.input.gas_network) # get clustered bus regions - buses_region = load_region(snakemake.input.regions_onshore, - snakemake.input.regions_offshore) - nodes = pd.Index(buses_region.name.unique()) + 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 = create_points2buses_map(points, buses_region) - # create gas network between pypsa nodes - cross_buses_gas_network = create_cross_regions_network(gas_network, - points2buses_map) + points2buses_map = points2buses(points, bus_regions) - # view(cross_buses_gas_network, buses_region, snakemake.input.country_shapes) + # create gas network between pypsa nodes + gas_connections = build_gas_network_topology(gas_network, points2buses_map) - # check which buses are not connected in gas network - check_missing(nodes, cross_buses_gas_network) + gas_connections = clean_dataset(nodes, gas_connections) - # convert units and save only needed data - gas_pipes = clean_dataset(cross_buses_gas_network) - gas_pipes.to_csv(snakemake.output.clustered_gas_network) + gas_connections.to_csv(snakemake.output.clustered_gas_network)