Merge pull request #835 from PyPSA/gas-storage

add locations, capacities and costs of existing gas storage

.pre-commit-config.yaml

@@ -50,7 +50,7 @@ repos:
   - id: blackdoc
 
   # Formatting with "black" coding style
-- repo: https://github.com/psf/black
+- repo: https://github.com/psf/black-pre-commit-mirror
   rev: 23.12.1
   hooks:
       # Format Python files

doc/release_notes.rst

@@ -75,6 +75,9 @@ Upcoming Release
 * Validate downloads from Zenodo using MD5 checksums. This identifies corrupted
   or incomplete downloads.
 
+* Add locations, capacities and costs of existing gas storage using Global
+  Energy Monitor's `Europe Gas Tracker
+  <https://globalenergymonitor.org/projects/europe-gas-tracker>`_.
 
 **Bugs and Compatibility**
 

rules/build_sector.smk

@@ -85,12 +85,12 @@ if config["sector"]["gas_network"] or config["sector"]["H2_retrofit"]:
 
     rule build_gas_input_locations:
         input:
-            lng=HTTP.remote(
+            gem=HTTP.remote(
                 "https://globalenergymonitor.org/wp-content/uploads/2023/07/Europe-Gas-Tracker-2023-03-v3.xlsx",
                 keep_local=True,
             ),
             entry="data/gas_network/scigrid-gas/data/IGGIELGN_BorderPoints.geojson",
-            production="data/gas_network/scigrid-gas/data/IGGIELGN_Productions.geojson",
+            storage="data/gas_network/scigrid-gas/data/IGGIELGN_Storages.geojson",
             regions_onshore=RESOURCES
             + "regions_onshore_elec_s{simpl}_{clusters}.geojson",
             regions_offshore=RESOURCES

rules/retrieve.smk

@@ -169,6 +169,7 @@ if config["enable"]["retrieve"] and (
         "IGGIELGN_LNGs.geojson",
         "IGGIELGN_BorderPoints.geojson",
         "IGGIELGN_Productions.geojson",
+        "IGGIELGN_Storages.geojson",
         "IGGIELGN_PipeSegments.geojson",
     ]
 

scripts/build_gas_input_locations.py

@@ -23,11 +23,10 @@ def read_scigrid_gas(fn):
     return df
 
 
-def build_gem_lng_data(lng_fn):
+def build_gem_lng_data(fn):
-    df = pd.read_excel(lng_fn[0], sheet_name="LNG terminals - data")
+    df = pd.read_excel(fn[0], sheet_name="LNG terminals - data")
     df = df.set_index("ComboID")
 
-    remove_status = ["Cancelled"]
     remove_country = ["Cyprus", "Turkey"]
     remove_terminal = ["Puerto de la Luz LNG Terminal", "Gran Canaria LNG Terminal"]
 
@@ -42,9 +41,50 @@ def build_gem_lng_data(lng_fn):
     return gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326")
 
 
-def build_gas_input_locations(lng_fn, entry_fn, prod_fn, countries):
+def build_gem_prod_data(fn):
+    df = pd.read_excel(fn[0], sheet_name="Gas extraction - main")
+    df = df.set_index("GEM Unit ID")
+
+    remove_country = ["Cyprus", "Türkiye"]
+    remove_fuel_type = ["oil"]
+
+    df = df.query(
+        "Status != 'shut in' \
+              & 'Fuel type' != 'oil' \
+              & Country != @remove_country \
+              & ~Latitude.isna() \
+              & ~Longitude.isna()"
+    ).copy()
+
+    p = pd.read_excel(fn[0], sheet_name="Gas extraction - production")
+    p = p.set_index("GEM Unit ID")
+    p = p[p["Fuel description"] == "gas"]
+
+    capacities = pd.DataFrame(index=df.index)
+    for key in ["production", "production design capacity", "reserves"]:
+        cap = (
+            p.loc[p["Production/reserves"] == key, "Quantity (converted)"]
+            .groupby("GEM Unit ID")
+            .sum()
+            .reindex(df.index)
+        )
+        # assume capacity such that 3% of reserves can be extracted per year (25% quantile)
+        annualization_factor = 0.03 if key == "reserves" else 1.0
+        capacities[key] = cap * annualization_factor
+
+    df["mcm_per_year"] = (
+        capacities["production"]
+        .combine_first(capacities["production design capacity"])
+        .combine_first(capacities["reserves"])
+    )
+
+    geometry = gpd.points_from_xy(df["Longitude"], df["Latitude"])
+    return gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326")
+
+
+def build_gas_input_locations(gem_fn, entry_fn, sto_fn, countries):
     # LNG terminals
-    lng = build_gem_lng_data(lng_fn)
+    lng = build_gem_lng_data(gem_fn)
 
     # Entry points from outside the model scope
     entry = read_scigrid_gas(entry_fn)
@@ -55,25 +95,30 @@ def build_gas_input_locations(lng_fn, entry_fn, prod_fn, countries):
         | (entry.from_country == "NO")  # malformed datapoint  # entries from NO to GB
     ]
 
+    sto = read_scigrid_gas(sto_fn)
+    remove_country = ["RU", "UA", "TR", "BY"]
+    sto = sto.query("country_code != @remove_country")
+
     # production sites inside the model scope
-    prod = read_scigrid_gas(prod_fn)
+    prod = build_gem_prod_data(gem_fn)
-    prod = prod.loc[
-        (prod.geometry.y > 35) & (prod.geometry.x < 30) & (prod.country_code != "DE")
-    ]
 
     mcm_per_day_to_mw = 437.5  # MCM/day to MWh/h
+    mcm_per_year_to_mw = 1.199  #  MCM/year to MWh/h
     mtpa_to_mw = 1649.224  # mtpa to MWh/h
-    lng["p_nom"] = lng["CapacityInMtpa"] * mtpa_to_mw
+    mcm_to_gwh = 11.36  # MCM to GWh
-    entry["p_nom"] = entry["max_cap_from_to_M_m3_per_d"] * mcm_per_day_to_mw
+    lng["capacity"] = lng["CapacityInMtpa"] * mtpa_to_mw
-    prod["p_nom"] = prod["max_supply_M_m3_per_d"] * mcm_per_day_to_mw
+    entry["capacity"] = entry["max_cap_from_to_M_m3_per_d"] * mcm_per_day_to_mw
+    prod["capacity"] = prod["mcm_per_year"] * mcm_per_year_to_mw
+    sto["capacity"] = sto["max_cushionGas_M_m3"] * mcm_to_gwh
 
     lng["type"] = "lng"
     entry["type"] = "pipeline"
     prod["type"] = "production"
+    sto["type"] = "storage"
 
-    sel = ["geometry", "p_nom", "type"]
+    sel = ["geometry", "capacity", "type"]
 
-    return pd.concat([prod[sel], entry[sel], lng[sel]], ignore_index=True)
+    return pd.concat([prod[sel], entry[sel], lng[sel], sto[sel]], ignore_index=True)
 
 
 if __name__ == "__main__":
@@ -83,7 +128,7 @@ if __name__ == "__main__":
         snakemake = mock_snakemake(
             "build_gas_input_locations",
             simpl="",
-            clusters="37",
+            clusters="128",
         )
 
     logging.basicConfig(level=snakemake.config["logging"]["level"])
@@ -104,9 +149,9 @@ if __name__ == "__main__":
     countries = regions.index.str[:2].unique().str.replace("GB", "UK")
 
     gas_input_locations = build_gas_input_locations(
-        snakemake.input.lng,
+        snakemake.input.gem,
         snakemake.input.entry,
-        snakemake.input.production,
+        snakemake.input.storage,
         countries,
     )
 
@@ -117,8 +162,8 @@ if __name__ == "__main__":
     gas_input_nodes.to_file(snakemake.output.gas_input_nodes, driver="GeoJSON")
 
     gas_input_nodes_s = (
-        gas_input_nodes.groupby(["bus", "type"])["p_nom"].sum().unstack()
+        gas_input_nodes.groupby(["bus", "type"])["capacity"].sum().unstack()
     )
-    gas_input_nodes_s.columns.name = "p_nom"
+    gas_input_nodes_s.columns.name = "capacity"
 
     gas_input_nodes_s.to_csv(snakemake.output.gas_input_nodes_simplified)

scripts/prepare_sector_network.py

@@ -454,10 +454,11 @@ def add_carrier_buses(n, carrier, nodes=None):
     n.add("Carrier", carrier)
 
     unit = "MWh_LHV" if carrier == "gas" else "MWh_th"
+    # preliminary value for non-gas carriers to avoid zeros
+    capital_cost = costs.at["gas storage", "fixed"] if carrier == "gas" else 0.02
 
     n.madd("Bus", nodes, location=location, carrier=carrier, unit=unit)
 
-    # capital cost could be corrected to e.g. 0.2 EUR/kWh * annuity and O&M
     n.madd(
         "Store",
         nodes + " Store",
@@ -465,8 +466,7 @@ def add_carrier_buses(n, carrier, nodes=None):
         e_nom_extendable=True,
         e_cyclic=True,
         carrier=carrier,
-        capital_cost=0.2
+        capital_cost=capital_cost,
-        * costs.at[carrier, "discount rate"],  # preliminary value to avoid zeros
     )
 
     n.madd(
@@ -1162,7 +1162,7 @@ def add_storage_and_grids(n, costs):
 
     if options["gas_network"]:
         logger.info(
-            "Add natural gas infrastructure, incl. LNG terminals, production and entry-points."
+            "Add natural gas infrastructure, incl. LNG terminals, production, storage and entry-points."
         )
 
         if options["H2_retrofit"]:
@@ -1207,10 +1207,25 @@ def add_storage_and_grids(n, costs):
         remove_i = n.generators[gas_i & internal_i].index
         n.generators.drop(remove_i, inplace=True)
 
-        p_nom = gas_input_nodes.sum(axis=1).rename(lambda x: x + " gas")
+        input_types = ["lng", "pipeline", "production"]
+        p_nom = gas_input_nodes[input_types].sum(axis=1).rename(lambda x: x + " gas")
         n.generators.loc[gas_i, "p_nom_extendable"] = False
         n.generators.loc[gas_i, "p_nom"] = p_nom
 
+        # add existing gas storage capacity
+        gas_i = n.stores.carrier == "gas"
+        e_nom = (
+            gas_input_nodes["storage"]
+            .rename(lambda x: x + " gas Store")
+            .reindex(n.stores.index)
+            .fillna(0.0)
+            * 1e3
+        )  # MWh_LHV
+        e_nom.clip(
+            upper=e_nom.quantile(0.98), inplace=True
+        )  # limit extremely large storage
+        n.stores.loc[gas_i, "e_nom_min"] = e_nom
+
         # add candidates for new gas pipelines to achieve full connectivity
 
         G = nx.Graph()