From fa03c61187a232c452714979515967910474f14b Mon Sep 17 00:00:00 2001
From: Fabian Neumann <fabian.neumann@outlook.de>
Date: Mon, 31 Jul 2023 10:52:37 +0200
Subject: [PATCH] gas_input: switch production data from scigrid to gem

---
 rules/build_sector.smk               |  3 +-
 scripts/build_gas_input_locations.py | 56 +++++++++++++++++++++-------
 2 files changed, 44 insertions(+), 15 deletions(-)

diff --git a/rules/build_sector.smk b/rules/build_sector.smk
index dd49fc6f..1e8c70ba 100644
--- a/rules/build_sector.smk
+++ b/rules/build_sector.smk
@@ -85,12 +85,11 @@ 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",
             regions_onshore=RESOURCES
             + "regions_onshore_elec_s{simpl}_{clusters}.geojson",
             regions_offshore=RESOURCES
diff --git a/scripts/build_gas_input_locations.py b/scripts/build_gas_input_locations.py
index a3b945ab..07707658 100644
--- a/scripts/build_gas_input_locations.py
+++ b/scripts/build_gas_input_locations.py
@@ -23,11 +23,10 @@ def read_scigrid_gas(fn):
     return df
 
 
-def build_gem_lng_data(lng_fn):
-    df = pd.read_excel(lng_fn[0], sheet_name="LNG terminals - data")
+def build_gem_lng_data(fn):
+    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,43 @@ 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. 
+        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, 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)
@@ -56,16 +89,14 @@ def build_gas_input_locations(lng_fn, entry_fn, prod_fn, countries):
     ]
 
     # production sites inside the model scope
-    prod = read_scigrid_gas(prod_fn)
-    prod = prod.loc[
-        (prod.geometry.y > 35) & (prod.geometry.x < 30) & (prod.country_code != "DE")
-    ]
+    prod = build_gem_prod_data(gem_fn)
 
     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
     entry["p_nom"] = entry["max_cap_from_to_M_m3_per_d"] * mcm_per_day_to_mw
-    prod["p_nom"] = prod["max_supply_M_m3_per_d"] * mcm_per_day_to_mw
+    prod["p_nom"] = prod["mcm_per_year"] * mcm_per_year_to_mw
 
     lng["type"] = "lng"
     entry["type"] = "pipeline"
@@ -83,7 +114,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 +135,8 @@ 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,
         countries,
     )