From 5e8d03cee60a81345bd60b0cd9ee85afcf776ab0 Mon Sep 17 00:00:00 2001 From: toniseibold Date: Wed, 25 Jun 2025 16:33:39 +0200 Subject: [PATCH 01/29] implementing existing ammonia plants --- config/config.default.yaml | 5 +- rules/solve_myopic.smk | 5 ++ rules/solve_perfect.smk | 5 ++ scripts/add_existing_baseyear.py | 79 ++++++++++++++++++++++++++++++++ 4 files changed, 92 insertions(+), 2 deletions(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index e7490e97ce..011cb72215 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -472,6 +472,7 @@ solar_thermal: existing_capacities: grouping_years_power: [1920, 1950, 1955, 1960, 1965, 1970, 1975, 1980, 1985, 1990, 1995, 2000, 2005, 2010, 2015, 2020, 2025] grouping_years_heat: [1980, 1985, 1990, 1995, 2000, 2005, 2010, 2015, 2019] # heat grouping years >= baseyear will be ignored + grouping_years_industry: [1995, 2000, 2005, 2010, 2015, 2020, 2025] threshold_capacity: 10 default_heating_lifetime: 20 conventional_carriers: @@ -812,7 +813,7 @@ sector: var_cf: true sustainability_factor: 0.0025 solid_biomass_import: - enable: false + enable: true price: 54 #EUR/MWh max_amount: 1390 # TWh upstream_emissions_factor: .1 #share of solid biomass CO2 emissions at full combustion @@ -968,7 +969,7 @@ clustering: ramp_limit_down: max temporal: resolution_elec: false - resolution_sector: false + resolution_sector: 8760H # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#adjustments adjustments: diff --git a/rules/solve_myopic.smk b/rules/solve_myopic.smk index df6637039e..96c678a824 100644 --- a/rules/solve_myopic.smk +++ b/rules/solve_myopic.smk @@ -12,6 +12,8 @@ rule add_existing_baseyear: costs=config_provider("costs"), heat_pump_sources=config_provider("sector", "heat_pump_sources"), energy_totals_year=config_provider("energy", "energy_totals_year"), + countries=config_provider("countries"), + MWh_NH3_per_tNH3=config_provider("industry", "MWh_NH3_per_tNH3"), input: network=resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc" @@ -28,6 +30,9 @@ rule add_existing_baseyear: "existing_heating_distribution_base_s_{clusters}_{planning_horizons}.csv" ), heating_efficiencies=resources("heating_efficiencies.csv"), + regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"), + ammonia="data/ammonia_plants.csv", + isi_database="data/1-s2.0-S0196890424010586-mmc2.xlsx", output: resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_brownfield.nc" diff --git a/rules/solve_perfect.smk b/rules/solve_perfect.smk index 8586d44366..8b5af71767 100644 --- a/rules/solve_perfect.smk +++ b/rules/solve_perfect.smk @@ -10,6 +10,8 @@ rule add_existing_baseyear: costs=config_provider("costs"), heat_pump_sources=config_provider("sector", "heat_pump_sources"), energy_totals_year=config_provider("energy", "energy_totals_year"), + countries=config_provider("countries"), + MWh_NH3_per_tNH3=config_provider("industry", "MWh_NH3_per_tNH3"), input: network=resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc" @@ -27,6 +29,9 @@ rule add_existing_baseyear: ), existing_heating="data/existing_infrastructure/existing_heating_raw.csv", heating_efficiencies=resources("heating_efficiencies.csv"), + regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"), + ammonia="data/ammonia_plants.csv", + isi_database="data/1-s2.0-S0196890424010586-mmc2.xlsx", output: resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_brownfield.nc" diff --git a/scripts/add_existing_baseyear.py b/scripts/add_existing_baseyear.py index 811ba0d414..58349f4790 100644 --- a/scripts/add_existing_baseyear.py +++ b/scripts/add_existing_baseyear.py @@ -13,6 +13,7 @@ import country_converter as coco import numpy as np import pandas as pd +import geopandas as gpd import powerplantmatching as pm import pypsa import xarray as xr @@ -716,6 +717,80 @@ def add_heating_capacities_installed_before_baseyear( ) +def prepare_plant_data(): + + # add existing industry + regions = gpd.read_file(snakemake.input.regions_onshore).set_index("name") + + isi_data = pd.read_excel(snakemake.input.isi_database, sheet_name="Database", index_col=1) + # assign bus region to each plant + geometry = gpd.points_from_xy(isi_data["Longitude"], isi_data["Latitude"]) + plant_data = gpd.GeoDataFrame(isi_data, geometry=geometry, crs="EPSG:4326") + plant_data = gpd.sjoin(plant_data, regions, how="inner", predicate="within") + plant_data.rename(columns={"name": "bus"}, inplace=True) + # filter for countries in model scope + plant_data = plant_data[plant_data.Country.isin(snakemake.params.countries)] + # replace UK with GB in Country column + plant_data["Country"] = plant_data["Country"].replace("UK", "GB") + # assign industry grouping year + grouping_years = snakemake.params.existing_capacities["grouping_years_industry"] + plant_data.loc[:, 'Year of last modernisation'] = plant_data['Year of last modernisation'].replace("x", np.nan) + plant_data['grouping_year'] = 0 + valid_mask = plant_data['Year of last modernisation'].notna() + valid_years = plant_data.loc[valid_mask, 'Year of last modernisation'] + indices = np.searchsorted(grouping_years, valid_years, side='right') + plant_data.loc[valid_years.index, 'grouping_year'] = np.array(grouping_years)[indices] + + return plant_data, regions + + +def add_existing_ammonia_plants(n): + + logger.info("Adding existing ammonia plants.") + + plant_data, regions = prepare_plant_data() + + ammonia_plants = plant_data[plant_data.Product=="Ammonia"] + + ammonia_plants = ammonia_plants.groupby(['bus', 'Country', 'grouping_year', 'Product'], as_index=False)['Production in tons (calibrated)'].sum() + + ammonia_plants.index = ammonia_plants['bus'] + " Haber-Bosch-" + ammonia_plants['grouping_year'].astype(str) + + # add dataset + df = pd.read_csv(snakemake.input.ammonia, index_col=0) + + geometry = gpd.points_from_xy(df.Longitude, df.Latitude) + gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326") + + gdf = gpd.sjoin(gdf, regions, how="inner", predicate="within") + + gdf.rename(columns={"name": "bus"}, inplace=True) + gdf["country"] = gdf.bus.str[:2] + # filter for countries that are missing + gdf[~gdf.country.isin(ammonia_plants.Country.unique())] + + n.add( + "Link", + ammonia_plants.index, + bus0=ammonia_plants.bus, + bus1=[bus + " NH3" for bus in ammonia_plants.bus] if snakemake.params.sector["ammonia"] else "EU NH3", + bus2=[bus + " gas" for bus in ammonia_plants.bus] if snakemake.params.sector["gas_network"] else "EU gas", + p_nom=ammonia_plants["Production in tons (calibrated)"].mul(snakemake.params.MWh_NH3_per_tNH3).div(costs.at["Haber-Bosch", "electricity-input"]).div(8760).values, + p_nom_extendable=False, + carrier="Haber-Bosch", + efficiency=1 / costs.at["Haber-Bosch", "electricity-input"], + efficiency2=-costs.at["Haber-Bosch", "hydrogen-input"] + / costs.at["Haber-Bosch", "electricity-input"], + capital_cost=costs.at["Haber-Bosch", "capital_cost"] + / costs.at["Haber-Bosch", "electricity-input"], + marginal_cost=costs.at["Haber-Bosch", "VOM"] + / costs.at["Haber-Bosch", "electricity-input"], + build_year=ammonia_plants["grouping_year"], + lifetime=costs.at["Haber-Bosch", "lifetime"], + ) + + + if __name__ == "__main__": if "snakemake" not in globals(): from scripts._helpers import mock_snakemake @@ -797,6 +872,10 @@ def add_heating_capacities_installed_before_baseyear( if options.get("cluster_heat_buses", False): cluster_heat_buses(n) + # add existing industry plants + if snakemake.sector.ammonia: + add_existing_ammonia_plants(n) + n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) sanitize_custom_columns(n) From 243ef47f14d54b3349dd673809d094410f88ccbb Mon Sep 17 00:00:00 2001 From: toniseibold Date: Fri, 27 Jun 2025 11:45:54 +0200 Subject: [PATCH 02/29] refining existing Haber-Bosch plants + Documentation --- scripts/add_existing_baseyear.py | 86 ++++++++++++++++++++++++-------- 1 file changed, 65 insertions(+), 21 deletions(-) diff --git a/scripts/add_existing_baseyear.py b/scripts/add_existing_baseyear.py index 58349f4790..692e06065f 100644 --- a/scripts/add_existing_baseyear.py +++ b/scripts/add_existing_baseyear.py @@ -717,12 +717,25 @@ def add_heating_capacities_installed_before_baseyear( ) -def prepare_plant_data(): +def prepare_plant_data( + regions_fn: str, + isi_database: str, +) -> tuple[pd.DataFrame, gpd.GeoDataFrame]: + """ + Reads in the Fraunhofer ISI database with high resolution plant data and maps them to the bus regions. + Returns the database as df as well as the regions as gdf. + Parameters + ---------- + regions_fn : str + path to the onshore regions file + isi_database: str + path to the fraunhofer isi database + """ # add existing industry - regions = gpd.read_file(snakemake.input.regions_onshore).set_index("name") + regions = gpd.read_file(regions_fn).set_index("name") - isi_data = pd.read_excel(snakemake.input.isi_database, sheet_name="Database", index_col=1) + isi_data = pd.read_excel(isi_database, sheet_name="Database", index_col=1) # assign bus region to each plant geometry = gpd.points_from_xy(isi_data["Longitude"], isi_data["Latitude"]) plant_data = gpd.GeoDataFrame(isi_data, geometry=geometry, crs="EPSG:4326") @@ -744,19 +757,27 @@ def prepare_plant_data(): return plant_data, regions -def add_existing_ammonia_plants(n): - +def add_existing_ammonia_plants( + n: pypsa.Network, +) -> None: + """ + Adds existing Haber-Bosch plants. + The plants are running on natural gas only since the retrofitting to hydrogen would be associated with costs. Plants are not expected to run at a minimal part load to avoid forcing the use of natural gas in planning horizons with climate targets. + Exhaust heat is not integrated since assuming that heat is integrated to make the current process more efficient. + """ logger.info("Adding existing ammonia plants.") - plant_data, regions = prepare_plant_data() + plant_data, regions = prepare_plant_data( + snakemake.input.regions_onshore, + snakemake.input.isi_database, + ) - ammonia_plants = plant_data[plant_data.Product=="Ammonia"] + fh_ammonia = plant_data[plant_data.Product=="Ammonia"] - ammonia_plants = ammonia_plants.groupby(['bus', 'Country', 'grouping_year', 'Product'], as_index=False)['Production in tons (calibrated)'].sum() + fh_ammonia = fh_ammonia.groupby(['bus', 'Country', 'grouping_year', 'Product'], as_index=False)['Production in tons (calibrated)'].sum() - ammonia_plants.index = ammonia_plants['bus'] + " Haber-Bosch-" + ammonia_plants['grouping_year'].astype(str) - - # add dataset + fh_ammonia.index = fh_ammonia['bus'] + " Haber-Bosch-SMR-" + fh_ammonia['grouping_year'].astype(str) + # add dataset for Non EU27 countries df = pd.read_csv(snakemake.input.ammonia, index_col=0) geometry = gpd.points_from_xy(df.Longitude, df.Latitude) @@ -765,22 +786,45 @@ def add_existing_ammonia_plants(n): gdf = gpd.sjoin(gdf, regions, how="inner", predicate="within") gdf.rename(columns={"name": "bus"}, inplace=True) - gdf["country"] = gdf.bus.str[:2] + gdf["Country"] = gdf.bus.str[:2] # filter for countries that are missing - gdf[~gdf.country.isin(ammonia_plants.Country.unique())] - + gdf = gdf[(~gdf.Country.isin(fh_ammonia.Country.unique())) & (gdf.Country.isin(snakemake.params.countries))] + # following approach from build_industrial_distribution_key.py + for country in gdf.Country: + facilities = gdf.query("Country == @country") + production = facilities["Ammonia [kt/a]"] + # assume 50% of the minimum production for missing values + production = production.fillna(0.5 * facilities["Ammonia [kt/a]"].min()) + + # missing data + gdf.drop(gdf[gdf["Ammonia [kt/a]"].isna()].index, inplace=True) + + # get average plant age: + avg_age = plant_data[plant_data.Product=="Ammonia"]["Year of last modernisation"].mean() + gdf["grouping_year"] = min((y for y in snakemake.params.existing_capacities["grouping_years_industry"] if y > avg_age)) + # match database + gdf.index = gdf["bus"] + " Haber-Bosch-SMR-" + gdf["grouping_year"].values.astype(str) + gdf.rename(columns={"Ammonia [kt/a]": "Production in tons (calibrated)"}, inplace=True) + gdf["Production in tons (calibrated)"] *= 1e3 + + ammonia_plants = pd.concat([fh_ammonia, gdf[["bus", "Country", "grouping_year", "Production in tons (calibrated)"]]]) + + # https://dechema.de/dechema_media/Downloads/Positionspapiere/Technology_study_Low_carbon_energy_and_feedstock_for_the_European_chemical_industry.pdf + # page 56: 1.83 t_CO2/t_NH3 + ch4_per_nh3 = 1.83 / costs.at["gas", "CO2 intensity"] / snakemake.params["MWh_NH3_per_tNH3"] n.add( "Link", ammonia_plants.index, - bus0=ammonia_plants.bus, + bus0=[bus + " gas" for bus in ammonia_plants.bus] if snakemake.params.sector["gas_network"] else "EU gas", bus1=[bus + " NH3" for bus in ammonia_plants.bus] if snakemake.params.sector["ammonia"] else "EU NH3", - bus2=[bus + " gas" for bus in ammonia_plants.bus] if snakemake.params.sector["gas_network"] else "EU gas", - p_nom=ammonia_plants["Production in tons (calibrated)"].mul(snakemake.params.MWh_NH3_per_tNH3).div(costs.at["Haber-Bosch", "electricity-input"]).div(8760).values, + bus2=ammonia_plants.bus, + bus3="co2 atmosphere", + p_nom=ammonia_plants["Production in tons (calibrated)"].mul(snakemake.params.MWh_NH3_per_tNH3).div(ch4_per_nh3).div(8760).values, p_nom_extendable=False, carrier="Haber-Bosch", - efficiency=1 / costs.at["Haber-Bosch", "electricity-input"], - efficiency2=-costs.at["Haber-Bosch", "hydrogen-input"] - / costs.at["Haber-Bosch", "electricity-input"], + efficiency=1 / ch4_per_nh3, + efficiency1=-costs.at["Haber-Bosch", "electricity-input"] / ch4_per_nh3, + efficiency2=costs.at["gas", "CO2 intensity"], capital_cost=costs.at["Haber-Bosch", "capital_cost"] / costs.at["Haber-Bosch", "electricity-input"], marginal_cost=costs.at["Haber-Bosch", "VOM"] @@ -873,7 +917,7 @@ def add_existing_ammonia_plants(n): cluster_heat_buses(n) # add existing industry plants - if snakemake.sector.ammonia: + if snakemake.params.sector["ammonia"]: add_existing_ammonia_plants(n) n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) From 643fcc96b8add573b5e7d5b7f20a20d8e61a868f Mon Sep 17 00:00:00 2001 From: toniseibold Date: Mon, 20 Oct 2025 18:02:32 +0200 Subject: [PATCH 03/29] endogenizing steel and cement production and add existing industry capacities --- config/config.default.yaml | 23 +- config/plotting.default.yaml | 18 ++ rules/build_sector.smk | 3 + rules/solve_myopic.smk | 2 + scripts/add_existing_baseyear.py | 379 +++++++++++++++++++++- scripts/prepare_sector_network.py | 504 +++++++++++++++++++++++++++++- 6 files changed, 914 insertions(+), 15 deletions(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index 011cb72215..7849ee54d6 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -29,19 +29,24 @@ run: use_shadow_directory: false # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#foresight -foresight: overnight +foresight: myopic # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#scenario # Wildcard docs in https://pypsa-eur.readthedocs.io/en/latest/wildcards.html scenario: clusters: - - 50 + - 39 + # - 128 + # - 256 opts: - '' sector_opts: - '' planning_horizons: - - 2050 + # - 2020 + - 2030 + - 2040 + # - 2050 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#countries countries: @@ -728,7 +733,7 @@ sector: - nearshore # within 50 km of sea # - offshore methanol: - regional_methanol_demand: false + regional_methanol_demand: true methanol_reforming: false methanol_reforming_cc: false methanol_to_kerosene: false @@ -739,7 +744,7 @@ sector: allam: false biomass_to_methanol: true biomass_to_methanol_cc: false - ammonia: true + ammonia: regional min_part_load_electrolysis: 0 min_part_load_fischer_tropsch: 0.5 min_part_load_methanolisation: 0.3 @@ -827,6 +832,12 @@ sector: methanol: 121 gas: 122 oil: 125 + endogenous_sectors: + - steel + - cement + - ammonia + - methanol + industry_relocation: false # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#industry industry: @@ -969,7 +980,7 @@ clustering: ramp_limit_down: max temporal: resolution_elec: false - resolution_sector: 8760H + resolution_sector: 365H # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#adjustments adjustments: diff --git a/config/plotting.default.yaml b/config/plotting.default.yaml index 5346f09258..f17e61a064 100644 --- a/config/plotting.default.yaml +++ b/config/plotting.default.yaml @@ -586,3 +586,21 @@ plotting: import NH3: '#e2ed74' import oil: '#93eda2' import methanol: '#87d0e6' + EAF: '#586357' + H2 DRI: '#618ab0' + gas DRI: '#baaf68' + BOF: '#1c1f1b' + steel: '#705238' + steel emission: '#756252' + steel emission CC: '#756252' + cement: '#c98349' + cement emission CC: '#c98349' + hbi: '#87d0e6' + DRI: '#87d0e6' + grey methanol: '#6b7887' + blue methanol: '#496687' + cement emission: '#c98349' + cement kiln: '#c98349' + cement production: '#c98349' + cement finishing: '#c98349' + clinker: '#c98349' \ No newline at end of file diff --git a/rules/build_sector.smk b/rules/build_sector.smk index 0ee8ddea58..8fd2bc60c2 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -1558,6 +1558,9 @@ rule prepare_sector_network: if config_provider("sector", "district_heating", "ates", "enable")(w) else [] ), + industry_sector_ratios=resources( + "industry_sector_ratios_{planning_horizons}.csv" + ), output: resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc" diff --git a/rules/solve_myopic.smk b/rules/solve_myopic.smk index 96c678a824..292aac986d 100644 --- a/rules/solve_myopic.smk +++ b/rules/solve_myopic.smk @@ -14,6 +14,7 @@ rule add_existing_baseyear: energy_totals_year=config_provider("energy", "energy_totals_year"), countries=config_provider("countries"), MWh_NH3_per_tNH3=config_provider("industry", "MWh_NH3_per_tNH3"), + MWh_MeOH_per_tMeOH=config_provider("industry", "MWh_MeOH_per_tMeOH"), input: network=resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}.nc" @@ -33,6 +34,7 @@ rule add_existing_baseyear: regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"), ammonia="data/ammonia_plants.csv", isi_database="data/1-s2.0-S0196890424010586-mmc2.xlsx", + gem_gspt="data/Global-Cement-and-Concrete-Tracker_July-2025.xlsx", output: resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_brownfield.nc" diff --git a/scripts/add_existing_baseyear.py b/scripts/add_existing_baseyear.py index 692e06065f..b607b7b895 100644 --- a/scripts/add_existing_baseyear.py +++ b/scripts/add_existing_baseyear.py @@ -14,6 +14,7 @@ import numpy as np import pandas as pd import geopandas as gpd +import pycountry import powerplantmatching as pm import pypsa import xarray as xr @@ -734,7 +735,104 @@ def prepare_plant_data( """ # add existing industry regions = gpd.read_file(regions_fn).set_index("name") + + isi_data = pd.read_excel(isi_database, sheet_name="Database", index_col=1) + # assign bus region to each plant + geometry = gpd.points_from_xy(isi_data["Longitude"], isi_data["Latitude"]) + plant_data = gpd.GeoDataFrame(isi_data, geometry=geometry, crs="EPSG:4326") + plant_data = gpd.sjoin(plant_data, regions, how="inner", predicate="within") + plant_data.rename(columns={"name": "bus"}, inplace=True) + # filter for countries in model scope + plant_data = plant_data[plant_data.Country.isin(snakemake.params.countries)] + # replace UK with GB in Country column + plant_data["Country"] = plant_data["Country"].replace("UK", "GB") + # assign industry grouping year + grouping_years = snakemake.params.existing_capacities["grouping_years_industry"] + plant_data.loc[:, "Year of last modernisation"] = plant_data[ + "Year of last modernisation" + ].replace("x", np.nan) + plant_data["grouping_year"] = 0 + valid_mask = plant_data["Year of last modernisation"].notna() + valid_years = plant_data.loc[valid_mask, "Year of last modernisation"] + indices = np.searchsorted(grouping_years, valid_years, side="right") + plant_data.loc[valid_years.index, "grouping_year"] = np.array(grouping_years)[ + indices + ] + + return plant_data, regions + + +def country_to_code(name): + try: + return pycountry.countries.lookup(name).alpha_2 + except LookupError: + return None + + +def prepare_gem_database(): + """ + Load GEM database of cement plants and map onto bus regions. + """ + regions = gpd.read_file(snakemake.input.regions_onshore).set_index("name") + + df = pd.read_excel( + snakemake.input.gem_gspt, + sheet_name="Plant Data", + na_values=["N/A", "unknown", ">0"], + ) + + df["country_code"] = df["Country/Area"].apply(country_to_code) + + df = df[df.country_code.isin(snakemake.params.countries)] + + df["Start date"] = pd.to_numeric( + df["Start date"].str.split("-").str[0], errors="coerce" + ) + + latlon = ( + df["Coordinates"] + .str.split(", ", expand=True) + .rename(columns={0: "lat", 1: "lon"}) + ) + geometry = gpd.points_from_xy(latlon["lon"], latlon["lat"]) + gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326") + + gdf = gpd.sjoin(gdf, regions, how="inner", predicate="within") + + gdf.rename(columns={"name": "bus"}, inplace=True) + gdf["country"] = gdf.bus.str[:2] + + grouping_years = snakemake.params.existing_capacities["grouping_years_industry"] + gdf = gdf[gdf["Start date"] < grouping_years[-1]] + avg_age = gdf["Start date"].mean() + gdf.fillna({"Start date": avg_age}, inplace=True) + gdf["grouping_year"] = 0 + indices = np.searchsorted(grouping_years, gdf["Start date"], side="right") + gdf.loc[:, "grouping_year"] = np.array(grouping_years)[indices] + gdf["grouping_year"].replace(2025, 2020, inplace=True) + + return gdf + + +def prepare_plant_data( + regions_fn: str, + isi_database: str, +) -> tuple[pd.DataFrame, gpd.GeoDataFrame]: + """ + Reads in the Fraunhofer ISI database with high resolution plant data and maps them to the bus regions. + Returns the database as df as well as the regions as gdf. + + Parameters + ---------- + regions_fn : str + path to the onshore regions file + isi_database: str + path to the fraunhofer isi database + """ + # add existing industry + regions = gpd.read_file(regions_fn).set_index("name") + isi_data = pd.read_excel(isi_database, sheet_name="Database", index_col=1) # assign bus region to each plant geometry = gpd.points_from_xy(isi_data["Longitude"], isi_data["Latitude"]) @@ -747,16 +845,275 @@ def prepare_plant_data( plant_data["Country"] = plant_data["Country"].replace("UK", "GB") # assign industry grouping year grouping_years = snakemake.params.existing_capacities["grouping_years_industry"] - plant_data.loc[:, 'Year of last modernisation'] = plant_data['Year of last modernisation'].replace("x", np.nan) - plant_data['grouping_year'] = 0 - valid_mask = plant_data['Year of last modernisation'].notna() - valid_years = plant_data.loc[valid_mask, 'Year of last modernisation'] - indices = np.searchsorted(grouping_years, valid_years, side='right') - plant_data.loc[valid_years.index, 'grouping_year'] = np.array(grouping_years)[indices] + plant_data.loc[:, "Year of last modernisation"] = plant_data[ + "Year of last modernisation" + ].replace("x", np.nan) + plant_data["grouping_year"] = 0 + valid_mask = plant_data["Year of last modernisation"].notna() + valid_years = plant_data.loc[valid_mask, "Year of last modernisation"] + indices = np.searchsorted(grouping_years, valid_years, side="right") + plant_data.loc[valid_years.index, "grouping_year"] = np.array(grouping_years)[ + indices + ] return plant_data, regions + +def add_existing_cement_plants(n): + database = prepare_gem_database() + # get plants that are still operating + database = database[database["Operating status"] == "operating"] + + cement = database.groupby( + ["bus", "country", "grouping_year"], as_index=False + )['Cement Capacity (millions metric tonnes per annum)'].sum() + + cement.index = (cement["bus"] + + " cement kiln-" + + cement["grouping_year"].astype(str) + ) + cement["clinker"] = (cement["bus"] + + " cement production-" + + cement["grouping_year"].astype(str) + ) + # get rid of capacities without demand + cement_buses = n.loads[n.loads.carrier=="cement"].bus.to_list() + cement_buses = [x.replace(" cement", "") for x in cement_buses] + cement = cement[cement.bus.isin(cement_buses)] + + logger.info(f"Adding {len(cement)} existing cement links.") + + # clinker production + gas_input = costs.at["cement dry clinker", "gas-input"] + costs.at["cement dry clinker", "heat-input"] + electricity_input = costs.at["cement dry clinker", "electricity-input"] + co2_emission = 0.79/gas_input + costs.at["gas", "CO2 intensity"] + + n.add( + "Link", + cement.index, + bus0=[bus + " gas" for bus in cement.bus] + if snakemake.params.sector["gas_network"] + else "EU gas", + bus1=[bus + " clinker" for bus in cement.bus], + bus2=cement.bus, + bus3=[bus + " cement emission" for bus in cement.bus], + carrier="cement kiln", + p_nom_extendable=False, + p_nom=cement['Cement Capacity (millions metric tonnes per annum)'] + .mul(gas_input) + .mul(costs.at["cement finishing", "clinker-input"]) + .div(8760) + .mul(1e6) + .values, + p_min_pu=0, + capital_cost=costs.at["cement dry clinker", "capital_cost"] / gas_input, + efficiency=1/gas_input, + efficiency2=-electricity_input, + efficiency3=co2_emission, + build_year=cement.grouping_year, + lifetime=costs.at["cement dry clinker", "lifetime"], + ) + + # cement finishing + electricity_input = costs.at["cement finishing", "electricity-input"] / costs.at["cement finishing", "clinker-input"] + clinker_input = costs.at["cement finishing", "clinker-input"] + n.add( + "Link", + cement.clinker.to_list(), + bus0=[bus + " clinker" for bus in cement.bus], + bus1=[bus + " cement" for bus in cement.bus], + bus2=cement.bus.to_list(), + carrier="cement finishing", + p_nom_extendable=False, + p_nom=cement['Cement Capacity (millions metric tonnes per annum)'] + .mul(clinker_input) + .div(8760) + .mul(1e6) + .values, + p_min_pu=0, + capital_cost=costs.at["cement finishing", "capital_cost"] / gas_input, + efficiency=1/clinker_input, + efficiency2=-electricity_input, + build_year=cement.grouping_year.to_list(), + lifetime=costs.at["cement finishing", "lifetime"], + ) + + +def add_existing_steel_plants( + n: pypsa.Network, +) -> None: + """ + Adds existing steel plants. + The plants are running on natural gas only since the retrofitting to hydrogen would be associated with costs. Plants are not expected to run at a minimal part load to avoid forcing the use of natural gas in planning horizons with climate targets. + Exhaust heat is not integrated since assuming that heat is integrated to make the current process more efficient. + """ + logger.info("Adding existing steel plants.") + + plant_data, regions = prepare_plant_data( + snakemake.input.regions_onshore, + snakemake.input.isi_database, + ) + steel_buses = n.loads[n.loads.carrier=="steel"].bus.to_list() + steel_buses = [x.replace(" steel", "") for x in steel_buses] + + fh_drg = plant_data[plant_data["Process status qup"] == "Direct reduction NG"] + drg = fh_drg.groupby( + ["bus", "Country", "grouping_year", "Product"], as_index=False + ).agg({ + "Production in tons (calibrated)": "sum", + "Out": "mean", + "Year of last modernisation": "mean", + }) + + fh_bof = plant_data[plant_data["Process status qup"] == "Blast furnace"] + bof = fh_bof.groupby( + ["bus", "Country", "grouping_year", "Product"], as_index=False + ).agg({ + "Production in tons (calibrated)": "sum", + "Out": "mean", + "Year of last modernisation": "mean", + }) + # fill Year of last modernisation NaNs with mean + mean_year = int(bof["Year of last modernisation"].mean()) + bof.loc[:, "Year of last modernisation"] = bof["Year of last modernisation"].fillna(mean_year) + drg.loc[:, "Year of last modernisation"] = drg["Year of last modernisation"].fillna(2025) + # fill in lifetime + bof.loc[:, "Lifetime"] = bof.loc[:, "Out"] - bof.loc[:, "Year of last modernisation"] + bof.loc[:, "Lifetime"] = bof["Lifetime"].fillna(costs.at["blast furnace-basic oxygen furnace", "lifetime"]) + drg.loc[:, "Lifetime"] = drg.loc[:, "Out"] - drg.loc[:, "Year of last modernisation"] + drg.loc[:, "Lifetime"] = drg["Lifetime"].fillna(costs.at["natural gas direct iron reduction furnace", "lifetime"]) + + drg.index = ( + drg["bus"] + + " gas DRI-" + + drg["grouping_year"].astype(str) + ) + drg = drg[drg.bus.isin(steel_buses)] + bof.index = ( + bof["bus"] + + " BOF-" + + bof["grouping_year"].astype(str) + ) + bof = bof[bof.bus.isin(steel_buses)] + + # add direct reduction with natural gas + gas_input = costs.at["natural gas direct iron reduction furnace", "gas-input"] + marginal_cost = ( + costs.at["iron ore DRI-ready", "commodity"] + * costs.at["natural gas direct iron reduction furnace", "ore-input"] + / gas_input + ) + + logger.info(f"Adding {len(drg)} gas DRI plant.") + + n.add( + "Link", + drg.index, + bus0=[bus + " gas" for bus in drg.bus] + if snakemake.params.sector["gas_network"] + else "EU gas", + bus1=[bus + " hbi" for bus in drg.bus] if not snakemake.params.sector["industry_relocation"] else "EU hbi", + bus2=[bus + " gas DRI emission" for bus in drg.bus], + p_nom=drg["Production in tons (calibrated)"] + .mul(gas_input) + .div(8760) + .values, + p_nom_extendable=False, + carrier="gas DRI", + efficiency=1/gas_input, + efficiency2=costs.at["gas", "CO2 intensity"], + capital_cost=costs.at["natural gas direct iron reduction furnace", "capital_cost"] + / gas_input, + marginal_cost=marginal_cost, + build_year=drg["Year of last modernisation"], + lifetime=drg["Lifetime"], + ) + + # add coal blast furnaces + coal_input = costs.at["blast furnace-basic oxygen furnace", "coal-input"] + marginal_cost = ( + costs.at["iron ore DRI-ready", "commodity"] + * costs.at["blast furnace-basic oxygen furnace", "ore-input"] / coal_input + ) + logger.info(f"Adding {len(bof)} BOF plant.") + + n.add( + "Link", + bof.index, + bus0="EU coal", + bus1=[bus + " steel" for bus in bof.bus], + bus2=[bus + " BOF emission" for bus in bof.bus], + p_nom=bof["Production in tons (calibrated)"] + .mul(coal_input) + .div(8760) + .values, + p_nom_extendable=False, + carrier="BOF", + efficiency=1/coal_input, + efficiency2=costs.at["coal", "CO2 intensity"], + marginal_cost=marginal_cost, + capital_cost=costs.at["blast furnace-basic oxygen furnace", "capital_cost"] / coal_input, + build_year=bof["Year of last modernisation"], + lifetime=bof["Lifetime"], + ) + + + +def add_existing_meoh_plants(n): + + logger.info("Adding existing methanol plants.") + + plant_data, regions = prepare_plant_data( + snakemake.input.regions_onshore, + snakemake.input.isi_database, + ) + + fh_meoh = plant_data[plant_data.Product == "Methanol"] + fh_meoh["grouping_year"].replace(2025, 2020, inplace=True) + fh_meoh = fh_meoh.groupby( + ["bus", "Country", "grouping_year", "Product"], as_index=False + )["Production in tons (calibrated)"].sum() + + fh_meoh.index = ( + fh_meoh["bus"] + + " grey methanol-" + + fh_meoh["grouping_year"].astype(str) + ) + + # grey methanol efficiency hard coded for now + costs.at["grey methanol synthesis", "efficiency"] = 1/1.757469244 + capital_cost = ( + costs.at["SMR", "capital_cost"] + + costs.at["methanolisation", "capital_cost"] + * costs.at["grey methanol synthesis", "efficiency"] + ) + co2_emissions = ( + costs.at["gas", "CO2 intensity"] + - costs.at["grey methanol synthesis", "efficiency"] + * costs.at["methanol", "CO2 intensity"] + ) + n.add( + "Link", + fh_meoh.index, + bus0=[bus + " gas" for bus in fh_meoh.bus] + if snakemake.params.sector["gas_network"] + else "EU gas", + bus1="EU methanol", + bus2="co2 atmosphere", + p_nom_extendable=False, + p_nom=fh_meoh["Production in tons (calibrated)"] + .mul(snakemake.params.MWh_MeOH_per_tMeOH) + .div(costs.at["grey methanol synthesis", "efficiency"]), + carrier="grey methanol", + efficiency=costs.at["grey methanol synthesis", "efficiency"], + efficiency2=co2_emissions, + capital_cost=capital_cost, + build_year=fh_meoh.grouping_year, + lifetime=costs.at["SMR", "lifetime"], + ) + + def add_existing_ammonia_plants( n: pypsa.Network, ) -> None: @@ -765,7 +1122,7 @@ def add_existing_ammonia_plants( The plants are running on natural gas only since the retrofitting to hydrogen would be associated with costs. Plants are not expected to run at a minimal part load to avoid forcing the use of natural gas in planning horizons with climate targets. Exhaust heat is not integrated since assuming that heat is integrated to make the current process more efficient. """ - logger.info("Adding existing ammonia plants.") + logger.info("Adding existing Haber-Bosch plants.") plant_data, regions = prepare_plant_data( snakemake.input.regions_onshore, @@ -917,8 +1274,14 @@ def add_existing_ammonia_plants( cluster_heat_buses(n) # add existing industry plants - if snakemake.params.sector["ammonia"]: + if "steel" in snakemake.params.sector["endogenous_sectors"]: + add_existing_steel_plants(n) + if "cement" in snakemake.params.sector["endogenous_sectors"]: + add_existing_cement_plants(n) + if ("ammonia" in snakemake.params.sector["endogenous_sectors"]) and (snakemake.params.sector["ammonia"]): add_existing_ammonia_plants(n) + if "methanol" in snakemake.params.sector["endogenous_sectors"]: + add_existing_meoh_plants(n) n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index fb0fb3de64..9f10e0aadf 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -220,6 +220,26 @@ def define_spatial(nodes, options): spatial.geothermal_heat.nodes = ["EU enhanced geothermal systems"] spatial.geothermal_heat.locations = ["EU"] + # industry + spatial.hbi = SimpleNamespace() + if options["industry_relocation"]: + spatial.hbi.nodes = ["EU hbi"] + spatial.hbi.location = ["EU"] + else: + spatial.hbi.nodes = nodes + " hbi" + spatial.hbi.location = nodes + + spatial.steel = SimpleNamespace() + spatial.steel.nodes = nodes + " steel" + spatial.steel.location = nodes + + spatial.cement = SimpleNamespace() + spatial.cement.nodes = nodes + " cement" + spatial.cement.location = nodes + spatial.clinker = SimpleNamespace() + spatial.clinker.nodes = nodes + " clinker" + spatial.clinker.location = nodes + return spatial @@ -4424,6 +4444,60 @@ def add_biomass( ) +def adjust_industry_demand(nodes): + # import ratios [MWh/t_Material] + fn = snakemake.input.industry_sector_ratios + sector_ratios = pd.read_csv(fn, header=[0, 1], index_col=0) + + # material demand per node and industry [kt/a] + fn = snakemake.input.industrial_production + nodal_production = pd.read_csv(fn, index_col=0) / 1e3 # kt/a -> Mt/a + + nodal_sector_ratios = pd.concat( + {node: sector_ratios[node[:2]] for node in nodal_production.index}, axis=1 + ) + + nodal_production_stacked = nodal_production.stack() + nodal_production_stacked.index.names = [None, None] + + # final energy consumption per node and industry (TWh/a) + nodal_df = (nodal_sector_ratios.multiply(nodal_production_stacked)).T + # rename the columns to correct unit + nodal_df.columns.name = "TWh/a" + + # add up all industry sectors except the ones endogenously modelled + endogenous_sector = snakemake.params.sector["endogenous_sectors"] + + sector_dict = { + "steel": ["Integrated steelworks", "DRI + Electric arc"], + "cement": "Cement", + } + keys = [sector_dict[k] for k in endogenous_sector if k in sector_dict] + + remaining_sectors = ~nodal_df.index.get_level_values(1).isin( + keys + ) + + remaining_demand = ( + nodal_df.loc[(nodes, remaining_sectors), :].groupby(level=0).sum() + ).mul(1e6) + remaining_demand.columns.name="MWh/a" + remaining_demand.rename(columns={ + "elec": "electricity", + "biomass": "solid biomass", + "heat": "low-temperature heat" + }, inplace=True) + # adjust initial demand + + # get the demand of steel and cement + steel = nodal_production[["Integrated steelworks", "DRI + Electric arc"]].sum(axis=1) + steel.name = "Mt/a" + cement = nodal_production["Cement"] + cement.name = "Mt/a" + + return remaining_demand, steel, cement + + def add_industry( n: pypsa.Network, costs: pd.DataFrame, @@ -4508,7 +4582,11 @@ def add_industry( nyears = nhours / 8760 # 1e6 to convert TWh to MWh - industrial_demand = pd.read_csv(industrial_demand_file, index_col=0) * 1e6 * nyears + industrial_demand_raw = pd.read_csv(industrial_demand_file, index_col=0) * 1e6 * nyears + # adjust the industrial energy demand excluding endogenously modelled carriers + adjusted_demand, steel, cement = adjust_industry_demand(nodes) + adjusted_demand["current electricity"] = industrial_demand_raw["current electricity"] + industrial_demand = adjusted_demand.copy() n.add( "Bus", @@ -5066,6 +5144,430 @@ def add_industry( efficiency2=costs.at["coal", "CO2 intensity"], ) + if "steel" in snakemake.params.sector["endogenous_sectors"]: + # add steel load + logger.info("Adding steel production technologies.") + n.add("Carrier", "steel", unit="t") + n.add("Carrier", "hbi", unit="t") + + steel.index = steel.index + " steel" + + n.add( + "Bus", + spatial.steel.nodes, + location=spatial.steel.location, + carrier="steel", + unit="t", + ) + + n.add( + "Bus", + spatial.hbi.nodes, + location=spatial.hbi.location, + carrier="hbi", + unit="t", + ) + + p_set = steel*1e6/8760 + n.add( + "Load", + spatial.steel.nodes, + bus=spatial.steel.nodes, + carrier="steel", + p_set=p_set, + ) + if investment_year > 2025: + electricity_input = costs.at["hydrogen direct iron reduction furnace", "electricity-input"] + hydrogen_input = costs.at["hydrogen direct iron reduction furnace", "hydrogen-input"] + + marginal_cost = ( + costs.at["iron ore DRI-ready", "commodity"] + * costs.at["hydrogen direct iron reduction furnace", "ore-input"] + / electricity_input + ) + n.add("Carrier", "H2 DRI") + + n.add( + "Link", + spatial.nodes, + suffix=" H2 DRI", + carrier="H2 DRI", + capital_cost=costs.at["hydrogen direct iron reduction furnace", "capital_cost"] + / electricity_input, + marginal_cost=marginal_cost, + p_nom=0, + p_nom_extendable=True, + p_min_pu=0, + bus0=spatial.nodes, + bus1=spatial.hbi.nodes, + bus2=spatial.h2.nodes, + efficiency=1 / electricity_input, + efficiency2=-hydrogen_input / electricity_input, + lifetime=costs.at["hydrogen direct iron reduction furnace", "lifetime"], + ) + + # HBI to steel via electric arc furnace + electricity_input = costs.at["electric arc furnace", "electricity-input"] + n.add("Carrier", "EAF") + + n.add( + "Link", + spatial.nodes, + suffix=" EAF", + carrier="EAF", + capital_cost=costs.at["electric arc furnace", "capital_cost"] / electricity_input, + p_nom=0, + p_nom_extendable=True, + p_min_pu=0, + bus0=spatial.nodes, + bus1=spatial.steel.nodes, + bus2=spatial.hbi.nodes, + efficiency=1 / electricity_input, + efficiency2=-costs.at["electric arc furnace", "hbi-input"] / electricity_input, + lifetime=costs.at["electric arc furnace", "lifetime"] + ) + # fossil options + # bus to decide whether to capture emissions or not + n.add("Carrier", "steel emission") + n.add( + "Bus", + spatial.nodes, + suffix=" gas DRI emission", + location=spatial.nodes, + carrier="steel emission", + unit="t", + ) + # NG DRI + gas_input = costs.at["natural gas direct iron reduction furnace", "gas-input"] + marginal_cost = ( + costs.at["iron ore DRI-ready", "commodity"] + * costs.at["natural gas direct iron reduction furnace", "ore-input"] + / gas_input + ) + + n.add("Carrier", "gas DRI") + n.add( + "Link", + spatial.nodes, + suffix=" gas DRI", + carrier="gas DRI", + capital_cost=costs.at["natural gas direct iron reduction furnace", "capital_cost"] + / gas_input, + marginal_cost=marginal_cost, + p_nom=0, + p_nom_extendable=True, + p_min_pu=0, + bus0=spatial.gas.nodes, + bus1=spatial.hbi.nodes, + bus2=spatial.nodes + " gas DRI emission", + efficiency=1/gas_input, + efficiency2=costs.at["gas", "CO2 intensity"], + lifetime=costs.at["natural gas direct iron reduction furnace", "lifetime"], + ) + + # BOF + coal_input = costs.at["blast furnace-basic oxygen furnace", "coal-input"] + marginal_cost = ( + costs.at["iron ore DRI-ready", "commodity"] + * costs.at["blast furnace-basic oxygen furnace", "ore-input"] / coal_input + ) + n.add("Carrier", "BOF") + n.add( + "Bus", + spatial.nodes, + suffix=" BOF emission", + location=spatial.nodes, + carrier="steel emission", + unit="t", + ) + n.add( + "Link", + spatial.nodes, + suffix=" BOF", + carrier="BOF", + capital_cost=costs.at["blast furnace-basic oxygen furnace", "capital_cost"] / coal_input, + p_nom=0, + p_nom_extendable=True, + p_min_pu=0, + marginal_cost=marginal_cost, + bus0="EU coal", + bus1=spatial.steel.nodes, + bus2=spatial.nodes + " BOF emission", + efficiency=1/coal_input, + efficiency2=costs.at["coal", "CO2 intensity"], + lifetime=costs.at["blast furnace-basic oxygen furnace", "lifetime"], + ) + # either capture or emitt co2 + n.add( + "Link", + spatial.nodes, + suffix=" gas DRI emission", + carrier="steel emission", + bus0=spatial.nodes + " gas DRI emission", + bus1="co2 atmosphere", + efficiency=1, + p_nom=0, + p_nom_extendable=True, + ) + n.add( + "Link", + spatial.nodes, + suffix=" BOF emission", + carrier="steel emission", + bus0=spatial.nodes + " BOF emission", + bus1="co2 atmosphere", + efficiency=1, + p_nom=0, + p_nom_extendable=True, + ) + if investment_year > 2025: + # Assumption: enough waste heat to recover sorbent + capture_rate = costs.at["cement capture", "capture_rate"] + electricity_input=costs.at["cement capture", "compression-electricity-input"] + costs.at["cement capture", "electricity-input"] # MWh/t_CO2 + n.add("Carrier", "steel emission CC") + n.add( + "Link", + spatial.nodes, + suffix=" gas DRI emission CC", + carrier="steel emission CC", + capital_cost=costs.at["cement capture", "capital_cost"], + bus0=spatial.nodes + " gas DRI emission", + bus1=spatial.co2.nodes, + bus2=spatial.nodes, + bus3="co2 atmosphere", + efficiency=capture_rate, + efficiency2=-electricity_input, + efficiency3=1-capture_rate, + p_nom=0, + p_nom_extendable=True, + lifetime=costs.at["cement capture", "lifetime"], + ) + n.add( + "Link", + spatial.nodes, + suffix=" BOF emission CC", + carrier="steel emission CC", + capital_cost=costs.at["cement capture", "capital_cost"], + bus0=spatial.nodes + " BOF emission", + bus1=spatial.co2.nodes, + bus2=spatial.nodes, + bus3="co2 atmosphere", + efficiency=capture_rate, + efficiency2=-electricity_input, + efficiency3=1-capture_rate, + p_nom=0, + p_nom_extendable=True, + lifetime=costs.at["cement capture", "lifetime"], + ) + + if "cement" in snakemake.params.sector["endogenous_sectors"]: + # add cement processes + logger.info("Adding cement production capacities.") + n.add("Carrier", "cement") + n.add("Carrier", "clinker") + + cement.index = cement.index + " cement" + + n.add( + "Bus", + spatial.cement.nodes, + location=spatial.cement.location, + carrier="cement", + unit="t", + ) + n.add( + "Bus", + spatial.clinker.nodes, + location=spatial.clinker.location, + carrier="clinker", + unit="t", + ) + + n.add( + "Bus", + spatial.cement.nodes, + suffix=" emission", + location=spatial.cement.location, + carrier="cement", + unit="t", + ) + + p_set = cement*1e6/8760 + n.add( + "Load", + cement.index, + bus=spatial.cement.nodes, + carrier="cement", + p_set=p_set, + ) + # clinker production + gas_input = costs.at["cement dry clinker", "gas-input"] + costs.at["cement dry clinker", "heat-input"] + electricity_input = costs.at["cement dry clinker", "electricity-input"] + co2_emission = 0.79/gas_input + costs.at["gas", "CO2 intensity"] + n.add("Carrier", "cement kiln") + n.add( + "Link", + spatial.clinker.nodes, + suffix = " kiln", + bus0=spatial.gas.nodes, + bus1=spatial.clinker.nodes, + bus2=spatial.nodes, + bus3=spatial.cement.nodes + " emission", + carrier="cement kiln", + p_nom_extendable=True, + p_min_pu=0, + capital_cost=costs.at["cement dry clinker", "capital_cost"] / gas_input, + efficiency=1/gas_input, + efficiency2=-electricity_input, + efficiency3=co2_emission, + lifetime=costs.at["cement dry clinker", "lifetime"], + ) + + # cement finishing + electricity_input = costs.at["cement finishing", "electricity-input"] / costs.at["cement finishing", "clinker-input"] + clinker_input = costs.at["cement finishing", "clinker-input"] + + n.add("Carrier", "cement finishing") + n.add( + "Link", + spatial.cement.nodes, + suffix = " production", + bus0=spatial.clinker.nodes, + bus1=spatial.cement.nodes, + bus2=spatial.nodes, + carrier="cement finishing", + p_nom_extendable=True, + p_min_pu=0, + capital_cost=costs.at["cement finishing", "capital_cost"] / gas_input, + efficiency=1/clinker_input, + efficiency2=-electricity_input, + lifetime=costs.at["cement finishing", "lifetime"], + ) + + if investment_year > 2025: + capture_rate = costs.at["cement capture", "capture_rate"] + electricity_input=costs.at["cement capture", "compression-electricity-input"] + + # add post combustion retrofit + n.add("Carrier", "cement emission CC") + n.add( + "Link", + spatial.cement.nodes, + suffix = " CC", + bus0=spatial.cement.nodes + " emission", + bus1=spatial.co2.nodes, + bus2=spatial.nodes, + bus3="co2 atmosphere", + carrier="cement emission CC", + p_nom_extendable=True, + p_min_pu=0, + capital_cost=costs.at["cement capture", "capital_cost"], + efficiency=costs.at["cement capture", "capture_rate"], + efficiency2=-electricity_input, + efficiency3=1-costs.at["cement capture", "capture_rate"], + lifetime=25, + ) + n.add("Carrier", "cement emission") + n.add( + "Link", + spatial.cement.nodes, + suffix = " emission", + bus0=spatial.cement.nodes + " emission", + bus1="co2 atmosphere", + carrier="cement emission", + p_nom_extendable=True, + capital_cost=0, + efficiency=1, + ) + + if investment_year > 2025: + # add methanol-to-olefins/aromatics + tech = "methanol-to-olefins/aromatics" + + naphtha_per_t_hvc = 12.622 # MWh per tonne of HVC (taken from sector ratios) + + process_emissions = ( + costs.at[tech, "carbondioxide-output"] / costs.at[tech, "methanol-input"] + ) + n.add("Carrier", "methanol-to-olefins/aromatics") + + n.add( + "Link", + nodes, + suffix=f" {tech}", + bus0=spatial.methanol.nodes, + bus1=spatial.oil.naphtha, + bus2=nodes, + bus3=spatial.co2.process_emissions, + efficiency=1 + / costs.at[tech, "methanol-input"] + * naphtha_per_t_hvc, # because MtO produces t HVC not MWh no MWh HVC + efficiency2=-costs.at[tech, "electricity-input"] + / costs.at[tech, "methanol-input"], + efficiency3=process_emissions, + carrier=tech, + p_nom_extendable=True, + ) + + costs.at["grey methanol synthesis", "efficiency"] = 1/1.757469244 + # grey methanol + capital_cost = ( + costs.at["SMR", "investment"] + + costs.at["methanolisation", "investment"] + * costs.at["grey methanol synthesis", "efficiency"] + ) + co2_emissions = ( + costs.at["gas", "CO2 intensity"] + - costs.at["grey methanol synthesis", "efficiency"] + * costs.at["methanol", "CO2 intensity"] + ) + n.add( + "Link", + spatial.methanol.demand_locations, + suffix=" grey methanol", + bus0=spatial.gas.nodes, + bus1=spatial.methanol.nodes, + bus2="co2 atmosphere", + p_nom_extendable=True, + carrier="grey methanol", + efficiency=costs.at["grey methanol synthesis", "efficiency"], + efficiency2=co2_emissions, + capital_cost=capital_cost, + lifetime=costs.at["SMR", "lifetime"], + ) + + # blue methanol + options["cc_fraction"] = 0.9 + co2_emissions = ( + costs.at["gas", "CO2 intensity"] + - costs.at["grey methanol synthesis", "efficiency"] + * costs.at["methanol", "CO2 intensity"] + ) + capital_cost = ( + costs.at["SMR", "investment"] + + costs.at["methanolisation", "investment"] + * costs.at["grey methanol synthesis", "efficiency"] + + costs.at["cement capture", "investment"] + * co2_emissions + * options["cc_fraction"] + ) + n.add( + "Link", + spatial.methanol.demand_locations, + suffix=" blue methanol", + bus0=spatial.gas.nodes, + bus1=spatial.methanol.nodes, + bus2="co2 atmosphere", + bus3=spatial.co2.nodes, + p_nom_extendable=True, + carrier="blue methanol", + efficiency=costs.at["grey methanol synthesis", "efficiency"], + efficiency2=co2_emissions * (1 - options["cc_fraction"]), + efficiency3=co2_emissions * options["cc_fraction"], + capital_cost=capital_cost, + lifetime=costs.at["SMR", "lifetime"], + ) + def add_aviation( n: pypsa.Network, From ef8e50f5dcbf49cb9eb05157b6170d0ef770add4 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 1 Dec 2025 13:39:50 +0000 Subject: [PATCH 04/29] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- config/config.default.yaml | 8 +- config/plotting.default.yaml | 2 +- scripts/add_existing_baseyear.py | 228 ++++++++++++++++++------------ scripts/prepare_sector_network.py | 180 +++++++++++++---------- 4 files changed, 244 insertions(+), 174 deletions(-) diff --git a/config/config.default.yaml b/config/config.default.yaml index 7849ee54d6..705ec6c471 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -833,10 +833,10 @@ sector: gas: 122 oil: 125 endogenous_sectors: - - steel - - cement - - ammonia - - methanol + - steel + - cement + - ammonia + - methanol industry_relocation: false # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#industry diff --git a/config/plotting.default.yaml b/config/plotting.default.yaml index f17e61a064..80f3d9b4fb 100644 --- a/config/plotting.default.yaml +++ b/config/plotting.default.yaml @@ -603,4 +603,4 @@ plotting: cement kiln: '#c98349' cement production: '#c98349' cement finishing: '#c98349' - clinker: '#c98349' \ No newline at end of file + clinker: '#c98349' diff --git a/scripts/add_existing_baseyear.py b/scripts/add_existing_baseyear.py index b607b7b895..9f3a00b972 100644 --- a/scripts/add_existing_baseyear.py +++ b/scripts/add_existing_baseyear.py @@ -11,11 +11,11 @@ from types import SimpleNamespace import country_converter as coco +import geopandas as gpd import numpy as np import pandas as pd -import geopandas as gpd -import pycountry import powerplantmatching as pm +import pycountry import pypsa import xarray as xr @@ -780,7 +780,7 @@ def prepare_gem_database(): sheet_name="Plant Data", na_values=["N/A", "unknown", ">0"], ) - + df["country_code"] = df["Country/Area"].apply(country_to_code) df = df[df.country_code.isin(snakemake.params.countries)] @@ -806,7 +806,7 @@ def prepare_gem_database(): gdf = gdf[gdf["Start date"] < grouping_years[-1]] avg_age = gdf["Start date"].mean() gdf.fillna({"Start date": avg_age}, inplace=True) - + gdf["grouping_year"] = 0 indices = np.searchsorted(grouping_years, gdf["Start date"], side="right") gdf.loc[:, "grouping_year"] = np.array(grouping_years)[indices] @@ -859,35 +859,33 @@ def prepare_plant_data( return plant_data, regions - def add_existing_cement_plants(n): database = prepare_gem_database() # get plants that are still operating database = database[database["Operating status"] == "operating"] - cement = database.groupby( - ["bus", "country", "grouping_year"], as_index=False - )['Cement Capacity (millions metric tonnes per annum)'].sum() + cement = database.groupby(["bus", "country", "grouping_year"], as_index=False)[ + "Cement Capacity (millions metric tonnes per annum)" + ].sum() - cement.index = (cement["bus"] - + " cement kiln-" - + cement["grouping_year"].astype(str) - ) - cement["clinker"] = (cement["bus"] - + " cement production-" - + cement["grouping_year"].astype(str) + cement.index = cement["bus"] + " cement kiln-" + cement["grouping_year"].astype(str) + cement["clinker"] = ( + cement["bus"] + " cement production-" + cement["grouping_year"].astype(str) ) # get rid of capacities without demand - cement_buses = n.loads[n.loads.carrier=="cement"].bus.to_list() + cement_buses = n.loads[n.loads.carrier == "cement"].bus.to_list() cement_buses = [x.replace(" cement", "") for x in cement_buses] cement = cement[cement.bus.isin(cement_buses)] logger.info(f"Adding {len(cement)} existing cement links.") # clinker production - gas_input = costs.at["cement dry clinker", "gas-input"] + costs.at["cement dry clinker", "heat-input"] + gas_input = ( + costs.at["cement dry clinker", "gas-input"] + + costs.at["cement dry clinker", "heat-input"] + ) electricity_input = costs.at["cement dry clinker", "electricity-input"] - co2_emission = 0.79/gas_input + costs.at["gas", "CO2 intensity"] + co2_emission = 0.79 / gas_input + costs.at["gas", "CO2 intensity"] n.add( "Link", @@ -900,7 +898,7 @@ def add_existing_cement_plants(n): bus3=[bus + " cement emission" for bus in cement.bus], carrier="cement kiln", p_nom_extendable=False, - p_nom=cement['Cement Capacity (millions metric tonnes per annum)'] + p_nom=cement["Cement Capacity (millions metric tonnes per annum)"] .mul(gas_input) .mul(costs.at["cement finishing", "clinker-input"]) .div(8760) @@ -908,7 +906,7 @@ def add_existing_cement_plants(n): .values, p_min_pu=0, capital_cost=costs.at["cement dry clinker", "capital_cost"] / gas_input, - efficiency=1/gas_input, + efficiency=1 / gas_input, efficiency2=-electricity_input, efficiency3=co2_emission, build_year=cement.grouping_year, @@ -916,7 +914,10 @@ def add_existing_cement_plants(n): ) # cement finishing - electricity_input = costs.at["cement finishing", "electricity-input"] / costs.at["cement finishing", "clinker-input"] + electricity_input = ( + costs.at["cement finishing", "electricity-input"] + / costs.at["cement finishing", "clinker-input"] + ) clinker_input = costs.at["cement finishing", "clinker-input"] n.add( "Link", @@ -926,14 +927,14 @@ def add_existing_cement_plants(n): bus2=cement.bus.to_list(), carrier="cement finishing", p_nom_extendable=False, - p_nom=cement['Cement Capacity (millions metric tonnes per annum)'] + p_nom=cement["Cement Capacity (millions metric tonnes per annum)"] .mul(clinker_input) .div(8760) .mul(1e6) .values, p_min_pu=0, capital_cost=costs.at["cement finishing", "capital_cost"] / gas_input, - efficiency=1/clinker_input, + efficiency=1 / clinker_input, efficiency2=-electricity_input, build_year=cement.grouping_year.to_list(), lifetime=costs.at["cement finishing", "lifetime"], @@ -954,56 +955,64 @@ def add_existing_steel_plants( snakemake.input.regions_onshore, snakemake.input.isi_database, ) - steel_buses = n.loads[n.loads.carrier=="steel"].bus.to_list() + steel_buses = n.loads[n.loads.carrier == "steel"].bus.to_list() steel_buses = [x.replace(" steel", "") for x in steel_buses] fh_drg = plant_data[plant_data["Process status qup"] == "Direct reduction NG"] drg = fh_drg.groupby( ["bus", "Country", "grouping_year", "Product"], as_index=False - ).agg({ - "Production in tons (calibrated)": "sum", - "Out": "mean", - "Year of last modernisation": "mean", - }) + ).agg( + { + "Production in tons (calibrated)": "sum", + "Out": "mean", + "Year of last modernisation": "mean", + } + ) fh_bof = plant_data[plant_data["Process status qup"] == "Blast furnace"] bof = fh_bof.groupby( ["bus", "Country", "grouping_year", "Product"], as_index=False - ).agg({ - "Production in tons (calibrated)": "sum", - "Out": "mean", - "Year of last modernisation": "mean", - }) + ).agg( + { + "Production in tons (calibrated)": "sum", + "Out": "mean", + "Year of last modernisation": "mean", + } + ) # fill Year of last modernisation NaNs with mean mean_year = int(bof["Year of last modernisation"].mean()) - bof.loc[:, "Year of last modernisation"] = bof["Year of last modernisation"].fillna(mean_year) - drg.loc[:, "Year of last modernisation"] = drg["Year of last modernisation"].fillna(2025) + bof.loc[:, "Year of last modernisation"] = bof["Year of last modernisation"].fillna( + mean_year + ) + drg.loc[:, "Year of last modernisation"] = drg["Year of last modernisation"].fillna( + 2025 + ) # fill in lifetime - bof.loc[:, "Lifetime"] = bof.loc[:, "Out"] - bof.loc[:, "Year of last modernisation"] - bof.loc[:, "Lifetime"] = bof["Lifetime"].fillna(costs.at["blast furnace-basic oxygen furnace", "lifetime"]) - drg.loc[:, "Lifetime"] = drg.loc[:, "Out"] - drg.loc[:, "Year of last modernisation"] - drg.loc[:, "Lifetime"] = drg["Lifetime"].fillna(costs.at["natural gas direct iron reduction furnace", "lifetime"]) - - drg.index = ( - drg["bus"] - + " gas DRI-" - + drg["grouping_year"].astype(str) + bof.loc[:, "Lifetime"] = ( + bof.loc[:, "Out"] - bof.loc[:, "Year of last modernisation"] ) - drg = drg[drg.bus.isin(steel_buses)] - bof.index = ( - bof["bus"] - + " BOF-" - + bof["grouping_year"].astype(str) + bof.loc[:, "Lifetime"] = bof["Lifetime"].fillna( + costs.at["blast furnace-basic oxygen furnace", "lifetime"] ) + drg.loc[:, "Lifetime"] = ( + drg.loc[:, "Out"] - drg.loc[:, "Year of last modernisation"] + ) + drg.loc[:, "Lifetime"] = drg["Lifetime"].fillna( + costs.at["natural gas direct iron reduction furnace", "lifetime"] + ) + + drg.index = drg["bus"] + " gas DRI-" + drg["grouping_year"].astype(str) + drg = drg[drg.bus.isin(steel_buses)] + bof.index = bof["bus"] + " BOF-" + bof["grouping_year"].astype(str) bof = bof[bof.bus.isin(steel_buses)] # add direct reduction with natural gas gas_input = costs.at["natural gas direct iron reduction furnace", "gas-input"] marginal_cost = ( - costs.at["iron ore DRI-ready", "commodity"] - * costs.at["natural gas direct iron reduction furnace", "ore-input"] - / gas_input - ) + costs.at["iron ore DRI-ready", "commodity"] + * costs.at["natural gas direct iron reduction furnace", "ore-input"] + / gas_input + ) logger.info(f"Adding {len(drg)} gas DRI plant.") @@ -1013,18 +1022,19 @@ def add_existing_steel_plants( bus0=[bus + " gas" for bus in drg.bus] if snakemake.params.sector["gas_network"] else "EU gas", - bus1=[bus + " hbi" for bus in drg.bus] if not snakemake.params.sector["industry_relocation"] else "EU hbi", + bus1=[bus + " hbi" for bus in drg.bus] + if not snakemake.params.sector["industry_relocation"] + else "EU hbi", bus2=[bus + " gas DRI emission" for bus in drg.bus], - p_nom=drg["Production in tons (calibrated)"] - .mul(gas_input) - .div(8760) - .values, + p_nom=drg["Production in tons (calibrated)"].mul(gas_input).div(8760).values, p_nom_extendable=False, carrier="gas DRI", - efficiency=1/gas_input, + efficiency=1 / gas_input, efficiency2=costs.at["gas", "CO2 intensity"], - capital_cost=costs.at["natural gas direct iron reduction furnace", "capital_cost"] - / gas_input, + capital_cost=costs.at[ + "natural gas direct iron reduction furnace", "capital_cost" + ] + / gas_input, marginal_cost=marginal_cost, build_year=drg["Year of last modernisation"], lifetime=drg["Lifetime"], @@ -1034,7 +1044,8 @@ def add_existing_steel_plants( coal_input = costs.at["blast furnace-basic oxygen furnace", "coal-input"] marginal_cost = ( costs.at["iron ore DRI-ready", "commodity"] - * costs.at["blast furnace-basic oxygen furnace", "ore-input"] / coal_input + * costs.at["blast furnace-basic oxygen furnace", "ore-input"] + / coal_input ) logger.info(f"Adding {len(bof)} BOF plant.") @@ -1044,24 +1055,20 @@ def add_existing_steel_plants( bus0="EU coal", bus1=[bus + " steel" for bus in bof.bus], bus2=[bus + " BOF emission" for bus in bof.bus], - p_nom=bof["Production in tons (calibrated)"] - .mul(coal_input) - .div(8760) - .values, + p_nom=bof["Production in tons (calibrated)"].mul(coal_input).div(8760).values, p_nom_extendable=False, carrier="BOF", - efficiency=1/coal_input, + efficiency=1 / coal_input, efficiency2=costs.at["coal", "CO2 intensity"], marginal_cost=marginal_cost, - capital_cost=costs.at["blast furnace-basic oxygen furnace", "capital_cost"] / coal_input, + capital_cost=costs.at["blast furnace-basic oxygen furnace", "capital_cost"] + / coal_input, build_year=bof["Year of last modernisation"], lifetime=bof["Lifetime"], ) - def add_existing_meoh_plants(n): - logger.info("Adding existing methanol plants.") plant_data, regions = prepare_plant_data( @@ -1076,13 +1083,11 @@ def add_existing_meoh_plants(n): )["Production in tons (calibrated)"].sum() fh_meoh.index = ( - fh_meoh["bus"] - + " grey methanol-" - + fh_meoh["grouping_year"].astype(str) + fh_meoh["bus"] + " grey methanol-" + fh_meoh["grouping_year"].astype(str) ) # grey methanol efficiency hard coded for now - costs.at["grey methanol synthesis", "efficiency"] = 1/1.757469244 + costs.at["grey methanol synthesis", "efficiency"] = 1 / 1.757469244 capital_cost = ( costs.at["SMR", "capital_cost"] + costs.at["methanolisation", "capital_cost"] @@ -1115,7 +1120,7 @@ def add_existing_meoh_plants(n): def add_existing_ammonia_plants( - n: pypsa.Network, + n: pypsa.Network, ) -> None: """ Adds existing Haber-Bosch plants. @@ -1129,11 +1134,17 @@ def add_existing_ammonia_plants( snakemake.input.isi_database, ) - fh_ammonia = plant_data[plant_data.Product=="Ammonia"] - - fh_ammonia = fh_ammonia.groupby(['bus', 'Country', 'grouping_year', 'Product'], as_index=False)['Production in tons (calibrated)'].sum() - - fh_ammonia.index = fh_ammonia['bus'] + " Haber-Bosch-SMR-" + fh_ammonia['grouping_year'].astype(str) + fh_ammonia = plant_data[plant_data.Product == "Ammonia"] + + fh_ammonia = fh_ammonia.groupby( + ["bus", "Country", "grouping_year", "Product"], as_index=False + )["Production in tons (calibrated)"].sum() + + fh_ammonia.index = ( + fh_ammonia["bus"] + + " Haber-Bosch-SMR-" + + fh_ammonia["grouping_year"].astype(str) + ) # add dataset for Non EU27 countries df = pd.read_csv(snakemake.input.ammonia, index_col=0) @@ -1143,9 +1154,12 @@ def add_existing_ammonia_plants( gdf = gpd.sjoin(gdf, regions, how="inner", predicate="within") gdf.rename(columns={"name": "bus"}, inplace=True) - gdf["Country"] = gdf.bus.str[:2] + gdf["Country"] = gdf.bus.str[:2] # filter for countries that are missing - gdf = gdf[(~gdf.Country.isin(fh_ammonia.Country.unique())) & (gdf.Country.isin(snakemake.params.countries))] + gdf = gdf[ + (~gdf.Country.isin(fh_ammonia.Country.unique())) + & (gdf.Country.isin(snakemake.params.countries)) + ] # following approach from build_industrial_distribution_key.py for country in gdf.Country: facilities = gdf.query("Country == @country") @@ -1157,26 +1171,51 @@ def add_existing_ammonia_plants( gdf.drop(gdf[gdf["Ammonia [kt/a]"].isna()].index, inplace=True) # get average plant age: - avg_age = plant_data[plant_data.Product=="Ammonia"]["Year of last modernisation"].mean() - gdf["grouping_year"] = min((y for y in snakemake.params.existing_capacities["grouping_years_industry"] if y > avg_age)) + avg_age = plant_data[plant_data.Product == "Ammonia"][ + "Year of last modernisation" + ].mean() + gdf["grouping_year"] = min( + y + for y in snakemake.params.existing_capacities["grouping_years_industry"] + if y > avg_age + ) # match database - gdf.index = gdf["bus"] + " Haber-Bosch-SMR-" + gdf["grouping_year"].values.astype(str) - gdf.rename(columns={"Ammonia [kt/a]": "Production in tons (calibrated)"}, inplace=True) + gdf.index = ( + gdf["bus"] + " Haber-Bosch-SMR-" + gdf["grouping_year"].values.astype(str) + ) + gdf.rename( + columns={"Ammonia [kt/a]": "Production in tons (calibrated)"}, inplace=True + ) gdf["Production in tons (calibrated)"] *= 1e3 - ammonia_plants = pd.concat([fh_ammonia, gdf[["bus", "Country", "grouping_year", "Production in tons (calibrated)"]]]) + ammonia_plants = pd.concat( + [ + fh_ammonia, + gdf[["bus", "Country", "grouping_year", "Production in tons (calibrated)"]], + ] + ) # https://dechema.de/dechema_media/Downloads/Positionspapiere/Technology_study_Low_carbon_energy_and_feedstock_for_the_European_chemical_industry.pdf # page 56: 1.83 t_CO2/t_NH3 - ch4_per_nh3 = 1.83 / costs.at["gas", "CO2 intensity"] / snakemake.params["MWh_NH3_per_tNH3"] + ch4_per_nh3 = ( + 1.83 / costs.at["gas", "CO2 intensity"] / snakemake.params["MWh_NH3_per_tNH3"] + ) n.add( "Link", ammonia_plants.index, - bus0=[bus + " gas" for bus in ammonia_plants.bus] if snakemake.params.sector["gas_network"] else "EU gas", - bus1=[bus + " NH3" for bus in ammonia_plants.bus] if snakemake.params.sector["ammonia"] else "EU NH3", + bus0=[bus + " gas" for bus in ammonia_plants.bus] + if snakemake.params.sector["gas_network"] + else "EU gas", + bus1=[bus + " NH3" for bus in ammonia_plants.bus] + if snakemake.params.sector["ammonia"] + else "EU NH3", bus2=ammonia_plants.bus, bus3="co2 atmosphere", - p_nom=ammonia_plants["Production in tons (calibrated)"].mul(snakemake.params.MWh_NH3_per_tNH3).div(ch4_per_nh3).div(8760).values, + p_nom=ammonia_plants["Production in tons (calibrated)"] + .mul(snakemake.params.MWh_NH3_per_tNH3) + .div(ch4_per_nh3) + .div(8760) + .values, p_nom_extendable=False, carrier="Haber-Bosch", efficiency=1 / ch4_per_nh3, @@ -1191,7 +1230,6 @@ def add_existing_ammonia_plants( ) - if __name__ == "__main__": if "snakemake" not in globals(): from scripts._helpers import mock_snakemake @@ -1278,7 +1316,9 @@ def add_existing_ammonia_plants( add_existing_steel_plants(n) if "cement" in snakemake.params.sector["endogenous_sectors"]: add_existing_cement_plants(n) - if ("ammonia" in snakemake.params.sector["endogenous_sectors"]) and (snakemake.params.sector["ammonia"]): + if ("ammonia" in snakemake.params.sector["endogenous_sectors"]) and ( + snakemake.params.sector["ammonia"] + ): add_existing_ammonia_plants(n) if "methanol" in snakemake.params.sector["endogenous_sectors"]: add_existing_meoh_plants(n) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index 9f10e0aadf..cf7b0abd6e 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -228,7 +228,7 @@ def define_spatial(nodes, options): else: spatial.hbi.nodes = nodes + " hbi" spatial.hbi.location = nodes - + spatial.steel = SimpleNamespace() spatial.steel.nodes = nodes + " steel" spatial.steel.location = nodes @@ -4445,57 +4445,60 @@ def add_biomass( def adjust_industry_demand(nodes): - # import ratios [MWh/t_Material] - fn = snakemake.input.industry_sector_ratios - sector_ratios = pd.read_csv(fn, header=[0, 1], index_col=0) + # import ratios [MWh/t_Material] + fn = snakemake.input.industry_sector_ratios + sector_ratios = pd.read_csv(fn, header=[0, 1], index_col=0) - # material demand per node and industry [kt/a] - fn = snakemake.input.industrial_production - nodal_production = pd.read_csv(fn, index_col=0) / 1e3 # kt/a -> Mt/a + # material demand per node and industry [kt/a] + fn = snakemake.input.industrial_production + nodal_production = pd.read_csv(fn, index_col=0) / 1e3 # kt/a -> Mt/a - nodal_sector_ratios = pd.concat( - {node: sector_ratios[node[:2]] for node in nodal_production.index}, axis=1 - ) + nodal_sector_ratios = pd.concat( + {node: sector_ratios[node[:2]] for node in nodal_production.index}, axis=1 + ) - nodal_production_stacked = nodal_production.stack() - nodal_production_stacked.index.names = [None, None] + nodal_production_stacked = nodal_production.stack() + nodal_production_stacked.index.names = [None, None] - # final energy consumption per node and industry (TWh/a) - nodal_df = (nodal_sector_ratios.multiply(nodal_production_stacked)).T - # rename the columns to correct unit - nodal_df.columns.name = "TWh/a" + # final energy consumption per node and industry (TWh/a) + nodal_df = (nodal_sector_ratios.multiply(nodal_production_stacked)).T + # rename the columns to correct unit + nodal_df.columns.name = "TWh/a" - # add up all industry sectors except the ones endogenously modelled - endogenous_sector = snakemake.params.sector["endogenous_sectors"] - - sector_dict = { - "steel": ["Integrated steelworks", "DRI + Electric arc"], - "cement": "Cement", - } - keys = [sector_dict[k] for k in endogenous_sector if k in sector_dict] + # add up all industry sectors except the ones endogenously modelled + endogenous_sector = snakemake.params.sector["endogenous_sectors"] - remaining_sectors = ~nodal_df.index.get_level_values(1).isin( - keys - ) + sector_dict = { + "steel": ["Integrated steelworks", "DRI + Electric arc"], + "cement": "Cement", + } + keys = [sector_dict[k] for k in endogenous_sector if k in sector_dict] + + remaining_sectors = ~nodal_df.index.get_level_values(1).isin(keys) - remaining_demand = ( - nodal_df.loc[(nodes, remaining_sectors), :].groupby(level=0).sum() - ).mul(1e6) - remaining_demand.columns.name="MWh/a" - remaining_demand.rename(columns={ + remaining_demand = ( + nodal_df.loc[(nodes, remaining_sectors), :].groupby(level=0).sum() + ).mul(1e6) + remaining_demand.columns.name = "MWh/a" + remaining_demand.rename( + columns={ "elec": "electricity", "biomass": "solid biomass", - "heat": "low-temperature heat" - }, inplace=True) - # adjust initial demand + "heat": "low-temperature heat", + }, + inplace=True, + ) + # adjust initial demand - # get the demand of steel and cement - steel = nodal_production[["Integrated steelworks", "DRI + Electric arc"]].sum(axis=1) - steel.name = "Mt/a" - cement = nodal_production["Cement"] - cement.name = "Mt/a" + # get the demand of steel and cement + steel = nodal_production[["Integrated steelworks", "DRI + Electric arc"]].sum( + axis=1 + ) + steel.name = "Mt/a" + cement = nodal_production["Cement"] + cement.name = "Mt/a" - return remaining_demand, steel, cement + return remaining_demand, steel, cement def add_industry( @@ -4582,10 +4585,14 @@ def add_industry( nyears = nhours / 8760 # 1e6 to convert TWh to MWh - industrial_demand_raw = pd.read_csv(industrial_demand_file, index_col=0) * 1e6 * nyears + industrial_demand_raw = ( + pd.read_csv(industrial_demand_file, index_col=0) * 1e6 * nyears + ) # adjust the industrial energy demand excluding endogenously modelled carriers adjusted_demand, steel, cement = adjust_industry_demand(nodes) - adjusted_demand["current electricity"] = industrial_demand_raw["current electricity"] + adjusted_demand["current electricity"] = industrial_demand_raw[ + "current electricity" + ] industrial_demand = adjusted_demand.copy() n.add( @@ -5168,7 +5175,7 @@ def add_industry( unit="t", ) - p_set = steel*1e6/8760 + p_set = steel * 1e6 / 8760 n.add( "Load", spatial.steel.nodes, @@ -5177,8 +5184,12 @@ def add_industry( p_set=p_set, ) if investment_year > 2025: - electricity_input = costs.at["hydrogen direct iron reduction furnace", "electricity-input"] - hydrogen_input = costs.at["hydrogen direct iron reduction furnace", "hydrogen-input"] + electricity_input = costs.at[ + "hydrogen direct iron reduction furnace", "electricity-input" + ] + hydrogen_input = costs.at[ + "hydrogen direct iron reduction furnace", "hydrogen-input" + ] marginal_cost = ( costs.at["iron ore DRI-ready", "commodity"] @@ -5192,7 +5203,9 @@ def add_industry( spatial.nodes, suffix=" H2 DRI", carrier="H2 DRI", - capital_cost=costs.at["hydrogen direct iron reduction furnace", "capital_cost"] + capital_cost=costs.at[ + "hydrogen direct iron reduction furnace", "capital_cost" + ] / electricity_input, marginal_cost=marginal_cost, p_nom=0, @@ -5215,7 +5228,8 @@ def add_industry( spatial.nodes, suffix=" EAF", carrier="EAF", - capital_cost=costs.at["electric arc furnace", "capital_cost"] / electricity_input, + capital_cost=costs.at["electric arc furnace", "capital_cost"] + / electricity_input, p_nom=0, p_nom_extendable=True, p_min_pu=0, @@ -5223,8 +5237,9 @@ def add_industry( bus1=spatial.steel.nodes, bus2=spatial.hbi.nodes, efficiency=1 / electricity_input, - efficiency2=-costs.at["electric arc furnace", "hbi-input"] / electricity_input, - lifetime=costs.at["electric arc furnace", "lifetime"] + efficiency2=-costs.at["electric arc furnace", "hbi-input"] + / electricity_input, + lifetime=costs.at["electric arc furnace", "lifetime"], ) # fossil options # bus to decide whether to capture emissions or not @@ -5251,7 +5266,9 @@ def add_industry( spatial.nodes, suffix=" gas DRI", carrier="gas DRI", - capital_cost=costs.at["natural gas direct iron reduction furnace", "capital_cost"] + capital_cost=costs.at[ + "natural gas direct iron reduction furnace", "capital_cost" + ] / gas_input, marginal_cost=marginal_cost, p_nom=0, @@ -5260,7 +5277,7 @@ def add_industry( bus0=spatial.gas.nodes, bus1=spatial.hbi.nodes, bus2=spatial.nodes + " gas DRI emission", - efficiency=1/gas_input, + efficiency=1 / gas_input, efficiency2=costs.at["gas", "CO2 intensity"], lifetime=costs.at["natural gas direct iron reduction furnace", "lifetime"], ) @@ -5269,7 +5286,8 @@ def add_industry( coal_input = costs.at["blast furnace-basic oxygen furnace", "coal-input"] marginal_cost = ( costs.at["iron ore DRI-ready", "commodity"] - * costs.at["blast furnace-basic oxygen furnace", "ore-input"] / coal_input + * costs.at["blast furnace-basic oxygen furnace", "ore-input"] + / coal_input ) n.add("Carrier", "BOF") n.add( @@ -5285,7 +5303,8 @@ def add_industry( spatial.nodes, suffix=" BOF", carrier="BOF", - capital_cost=costs.at["blast furnace-basic oxygen furnace", "capital_cost"] / coal_input, + capital_cost=costs.at["blast furnace-basic oxygen furnace", "capital_cost"] + / coal_input, p_nom=0, p_nom_extendable=True, p_min_pu=0, @@ -5293,7 +5312,7 @@ def add_industry( bus0="EU coal", bus1=spatial.steel.nodes, bus2=spatial.nodes + " BOF emission", - efficiency=1/coal_input, + efficiency=1 / coal_input, efficiency2=costs.at["coal", "CO2 intensity"], lifetime=costs.at["blast furnace-basic oxygen furnace", "lifetime"], ) @@ -5323,7 +5342,10 @@ def add_industry( if investment_year > 2025: # Assumption: enough waste heat to recover sorbent capture_rate = costs.at["cement capture", "capture_rate"] - electricity_input=costs.at["cement capture", "compression-electricity-input"] + costs.at["cement capture", "electricity-input"] # MWh/t_CO2 + electricity_input = ( + costs.at["cement capture", "compression-electricity-input"] + + costs.at["cement capture", "electricity-input"] + ) # MWh/t_CO2 n.add("Carrier", "steel emission CC") n.add( "Link", @@ -5337,7 +5359,7 @@ def add_industry( bus3="co2 atmosphere", efficiency=capture_rate, efficiency2=-electricity_input, - efficiency3=1-capture_rate, + efficiency3=1 - capture_rate, p_nom=0, p_nom_extendable=True, lifetime=costs.at["cement capture", "lifetime"], @@ -5354,7 +5376,7 @@ def add_industry( bus3="co2 atmosphere", efficiency=capture_rate, efficiency2=-electricity_input, - efficiency3=1-capture_rate, + efficiency3=1 - capture_rate, p_nom=0, p_nom_extendable=True, lifetime=costs.at["cement capture", "lifetime"], @@ -5392,7 +5414,7 @@ def add_industry( unit="t", ) - p_set = cement*1e6/8760 + p_set = cement * 1e6 / 8760 n.add( "Load", cement.index, @@ -5401,14 +5423,17 @@ def add_industry( p_set=p_set, ) # clinker production - gas_input = costs.at["cement dry clinker", "gas-input"] + costs.at["cement dry clinker", "heat-input"] + gas_input = ( + costs.at["cement dry clinker", "gas-input"] + + costs.at["cement dry clinker", "heat-input"] + ) electricity_input = costs.at["cement dry clinker", "electricity-input"] - co2_emission = 0.79/gas_input + costs.at["gas", "CO2 intensity"] + co2_emission = 0.79 / gas_input + costs.at["gas", "CO2 intensity"] n.add("Carrier", "cement kiln") n.add( "Link", spatial.clinker.nodes, - suffix = " kiln", + suffix=" kiln", bus0=spatial.gas.nodes, bus1=spatial.clinker.nodes, bus2=spatial.nodes, @@ -5417,21 +5442,24 @@ def add_industry( p_nom_extendable=True, p_min_pu=0, capital_cost=costs.at["cement dry clinker", "capital_cost"] / gas_input, - efficiency=1/gas_input, + efficiency=1 / gas_input, efficiency2=-electricity_input, efficiency3=co2_emission, lifetime=costs.at["cement dry clinker", "lifetime"], ) # cement finishing - electricity_input = costs.at["cement finishing", "electricity-input"] / costs.at["cement finishing", "clinker-input"] + electricity_input = ( + costs.at["cement finishing", "electricity-input"] + / costs.at["cement finishing", "clinker-input"] + ) clinker_input = costs.at["cement finishing", "clinker-input"] n.add("Carrier", "cement finishing") n.add( "Link", spatial.cement.nodes, - suffix = " production", + suffix=" production", bus0=spatial.clinker.nodes, bus1=spatial.cement.nodes, bus2=spatial.nodes, @@ -5439,21 +5467,23 @@ def add_industry( p_nom_extendable=True, p_min_pu=0, capital_cost=costs.at["cement finishing", "capital_cost"] / gas_input, - efficiency=1/clinker_input, + efficiency=1 / clinker_input, efficiency2=-electricity_input, lifetime=costs.at["cement finishing", "lifetime"], ) if investment_year > 2025: capture_rate = costs.at["cement capture", "capture_rate"] - electricity_input=costs.at["cement capture", "compression-electricity-input"] - + electricity_input = costs.at[ + "cement capture", "compression-electricity-input" + ] + # add post combustion retrofit n.add("Carrier", "cement emission CC") n.add( "Link", spatial.cement.nodes, - suffix = " CC", + suffix=" CC", bus0=spatial.cement.nodes + " emission", bus1=spatial.co2.nodes, bus2=spatial.nodes, @@ -5464,14 +5494,14 @@ def add_industry( capital_cost=costs.at["cement capture", "capital_cost"], efficiency=costs.at["cement capture", "capture_rate"], efficiency2=-electricity_input, - efficiency3=1-costs.at["cement capture", "capture_rate"], + efficiency3=1 - costs.at["cement capture", "capture_rate"], lifetime=25, ) n.add("Carrier", "cement emission") n.add( "Link", spatial.cement.nodes, - suffix = " emission", + suffix=" emission", bus0=spatial.cement.nodes + " emission", bus1="co2 atmosphere", carrier="cement emission", @@ -5500,16 +5530,16 @@ def add_industry( bus2=nodes, bus3=spatial.co2.process_emissions, efficiency=1 - / costs.at[tech, "methanol-input"] - * naphtha_per_t_hvc, # because MtO produces t HVC not MWh no MWh HVC + / costs.at[tech, "methanol-input"] + * naphtha_per_t_hvc, # because MtO produces t HVC not MWh no MWh HVC efficiency2=-costs.at[tech, "electricity-input"] - / costs.at[tech, "methanol-input"], + / costs.at[tech, "methanol-input"], efficiency3=process_emissions, carrier=tech, p_nom_extendable=True, ) - costs.at["grey methanol synthesis", "efficiency"] = 1/1.757469244 + costs.at["grey methanol synthesis", "efficiency"] = 1 / 1.757469244 # grey methanol capital_cost = ( costs.at["SMR", "investment"] From 5af80f01b5643fb57bf3e73dd7afb11ceba1a779 Mon Sep 17 00:00:00 2001 From: toniseibold Date: Mon, 1 Dec 2025 16:29:26 +0100 Subject: [PATCH 05/29] cleaning up plotting colors --- config/plotting.default.yaml | 26 ++++++++++++-------------- 1 file changed, 12 insertions(+), 14 deletions(-) diff --git a/config/plotting.default.yaml b/config/plotting.default.yaml index 80f3d9b4fb..9f038f8ec4 100644 --- a/config/plotting.default.yaml +++ b/config/plotting.default.yaml @@ -586,21 +586,19 @@ plotting: import NH3: '#e2ed74' import oil: '#93eda2' import methanol: '#87d0e6' + steel: '#22333B' + steel emission: '#5E503F' + steel emission CC: '#A9927D' + hbi: '#CCDBDC' + DRI: '#8E7C93' EAF: '#586357' H2 DRI: '#618ab0' gas DRI: '#baaf68' - BOF: '#1c1f1b' - steel: '#705238' - steel emission: '#756252' - steel emission CC: '#756252' + BOF: '#0A0908' + grey methanol: '#7C6C77' cement: '#c98349' - cement emission CC: '#c98349' - hbi: '#87d0e6' - DRI: '#87d0e6' - grey methanol: '#6b7887' - blue methanol: '#496687' - cement emission: '#c98349' - cement kiln: '#c98349' - cement production: '#c98349' - cement finishing: '#c98349' - clinker: '#c98349' + cement emission: '#119DA4' + cement emission CC: '#0C7489' + cement kiln: '#D90368' + cement finishing: '#820263' + clinker: '#F5FDC6' \ No newline at end of file From 018cf57136dec6520c33355850b742d262798826 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 1 Dec 2025 15:30:37 +0000 Subject: [PATCH 06/29] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- config/plotting.default.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/config/plotting.default.yaml b/config/plotting.default.yaml index 9f038f8ec4..703b647506 100644 --- a/config/plotting.default.yaml +++ b/config/plotting.default.yaml @@ -601,4 +601,4 @@ plotting: cement emission CC: '#0C7489' cement kiln: '#D90368' cement finishing: '#820263' - clinker: '#F5FDC6' \ No newline at end of file + clinker: '#F5FDC6' From 870ca70b188b0360d3bf7c6d129bdbf2aee0227e Mon Sep 17 00:00:00 2001 From: toniseibold Date: Thu, 4 Dec 2025 16:00:28 +0100 Subject: [PATCH 07/29] own rule for industry plant data preparation, restructuring of the config for existing cement and steel plants --- config/config.default.yaml | 10 +- doc/configtables/sector.csv | 7 + rules/build_sector.smk | 21 + rules/solve_myopic.smk | 6 +- scripts/add_existing_baseyear.py | 701 +++++++++--------------------- scripts/build_industry_plants.py | 199 +++++++++ scripts/prepare_sector_network.py | 59 +-- 7 files changed, 461 insertions(+), 542 deletions(-) create mode 100644 scripts/build_industry_plants.py diff --git a/config/config.default.yaml b/config/config.default.yaml index 705ec6c471..da98011444 100644 --- a/config/config.default.yaml +++ b/config/config.default.yaml @@ -93,7 +93,7 @@ snapshots: # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#enable enable: - retrieve: auto + retrieve: false retrieve_databundle: true retrieve_cost_data: true build_cutout: false @@ -835,9 +835,11 @@ sector: endogenous_sectors: - steel - cement - - ammonia - - methanol - industry_relocation: false + hbi_relocation: false + steel_bof: + pledge: true + pledge_delay: 0 + default_phase_out: 2035 # docs in https://pypsa-eur.readthedocs.io/en/latest/configuration.html#industry industry: diff --git a/doc/configtables/sector.csv b/doc/configtables/sector.csv index dbc75504a0..fd579b7623 100644 --- a/doc/configtables/sector.csv +++ b/doc/configtables/sector.csv @@ -228,3 +228,10 @@ imports,,, -- limit_sense,--,"{==, <=, >=}",Sense of the limit -- price,,"{H2, NH3, methanol, gas, oil}", -- -- {carrier},currency/MWh,float,Price for importing renewable energy of carrier +endogenous_sectors,,, +-- {steel, cement}, Endogenously model given subsectors +hbi_relocation,--,"{true, false}",Allow the relocation of iron ore reduction step +steel_bof,,, +-- pledge,--,"{true, false}",decommissions coal blast furnaces at their pledged phase out year +-- pledge_delay,--,int,delays the decommission year by number of years +-- default_phase_out,--,{false, int}",decommissions the coal blast furnaces that do not have a pledged phase out year diff --git a/rules/build_sector.smk b/rules/build_sector.smk index 8fd2bc60c2..dcd26044ca 100755 --- a/rules/build_sector.smk +++ b/rules/build_sector.smk @@ -1402,6 +1402,27 @@ def input_heat_source_power(w): } +rule build_industry_plants: + params: + countries=config_provider("countries"), + input: + regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"), + ammonia="data/ammonia_plants.csv", + isi_database="data/1-s2.0-S0196890424010586-mmc2.xlsx", + gem_gcct="data/gem/Global-Cement-and-Concrete-Tracker_July-2025.xlsx", + output: + industry_plants=resources("industry_plants_{clusters}.csv"), + threads: 1 + resources: + mem_mb=2000, + log: + logs("build_industry_plants_{clusters}.log"), + benchmark: + benchmarks("build_industry_plants_{clusters}") + script: + "../scripts/build_industry_plants.py" + + rule prepare_sector_network: params: time_resolution=config_provider("clustering", "temporal", "resolution_sector"), diff --git a/rules/solve_myopic.smk b/rules/solve_myopic.smk index 292aac986d..58f8d4e0e1 100644 --- a/rules/solve_myopic.smk +++ b/rules/solve_myopic.smk @@ -12,7 +12,6 @@ rule add_existing_baseyear: costs=config_provider("costs"), heat_pump_sources=config_provider("sector", "heat_pump_sources"), energy_totals_year=config_provider("energy", "energy_totals_year"), - countries=config_provider("countries"), MWh_NH3_per_tNH3=config_provider("industry", "MWh_NH3_per_tNH3"), MWh_MeOH_per_tMeOH=config_provider("industry", "MWh_MeOH_per_tMeOH"), input: @@ -31,10 +30,7 @@ rule add_existing_baseyear: "existing_heating_distribution_base_s_{clusters}_{planning_horizons}.csv" ), heating_efficiencies=resources("heating_efficiencies.csv"), - regions_onshore=resources("regions_onshore_base_s_{clusters}.geojson"), - ammonia="data/ammonia_plants.csv", - isi_database="data/1-s2.0-S0196890424010586-mmc2.xlsx", - gem_gspt="data/Global-Cement-and-Concrete-Tracker_July-2025.xlsx", + industry_plants=resources("industry_plants_{clusters}.csv"), output: resources( "networks/base_s_{clusters}_{opts}_{sector_opts}_{planning_horizons}_brownfield.nc" diff --git a/scripts/add_existing_baseyear.py b/scripts/add_existing_baseyear.py index 9f3a00b972..74befd62d5 100644 --- a/scripts/add_existing_baseyear.py +++ b/scripts/add_existing_baseyear.py @@ -718,516 +718,251 @@ def add_heating_capacities_installed_before_baseyear( ) -def prepare_plant_data( - regions_fn: str, - isi_database: str, -) -> tuple[pd.DataFrame, gpd.GeoDataFrame]: - """ - Reads in the Fraunhofer ISI database with high resolution plant data and maps them to the bus regions. - Returns the database as df as well as the regions as gdf. - - Parameters - ---------- - regions_fn : str - path to the onshore regions file - isi_database: str - path to the fraunhofer isi database - """ - # add existing industry - regions = gpd.read_file(regions_fn).set_index("name") - - isi_data = pd.read_excel(isi_database, sheet_name="Database", index_col=1) - # assign bus region to each plant - geometry = gpd.points_from_xy(isi_data["Longitude"], isi_data["Latitude"]) - plant_data = gpd.GeoDataFrame(isi_data, geometry=geometry, crs="EPSG:4326") - plant_data = gpd.sjoin(plant_data, regions, how="inner", predicate="within") - plant_data.rename(columns={"name": "bus"}, inplace=True) - # filter for countries in model scope - plant_data = plant_data[plant_data.Country.isin(snakemake.params.countries)] - # replace UK with GB in Country column - plant_data["Country"] = plant_data["Country"].replace("UK", "GB") - # assign industry grouping year - grouping_years = snakemake.params.existing_capacities["grouping_years_industry"] - plant_data.loc[:, "Year of last modernisation"] = plant_data[ - "Year of last modernisation" - ].replace("x", np.nan) - plant_data["grouping_year"] = 0 - valid_mask = plant_data["Year of last modernisation"].notna() - valid_years = plant_data.loc[valid_mask, "Year of last modernisation"] - indices = np.searchsorted(grouping_years, valid_years, side="right") - plant_data.loc[valid_years.index, "grouping_year"] = np.array(grouping_years)[ - indices - ] - - return plant_data, regions - - -def country_to_code(name): - try: - return pycountry.countries.lookup(name).alpha_2 - except LookupError: - return None - - -def prepare_gem_database(): - """ - Load GEM database of cement plants and map onto bus regions. - """ - regions = gpd.read_file(snakemake.input.regions_onshore).set_index("name") - - df = pd.read_excel( - snakemake.input.gem_gspt, - sheet_name="Plant Data", - na_values=["N/A", "unknown", ">0"], +def add_existing_industry(n): + plant_data = pd.read_csv(snakemake.input.industry_plants) + + # fill missing build_year with average + mean_nonzero = (plant_data + .groupby("carrier")["build_year"] + .transform("mean") + .round() ) - - df["country_code"] = df["Country/Area"].apply(country_to_code) - - df = df[df.country_code.isin(snakemake.params.countries)] - - df["Start date"] = pd.to_numeric( - df["Start date"].str.split("-").str[0], errors="coerce" - ) - - latlon = ( - df["Coordinates"] - .str.split(", ", expand=True) - .rename(columns={0: "lat", 1: "lon"}) - ) - geometry = gpd.points_from_xy(latlon["lon"], latlon["lat"]) - gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326") - - gdf = gpd.sjoin(gdf, regions, how="inner", predicate="within") - - gdf.rename(columns={"name": "bus"}, inplace=True) - gdf["country"] = gdf.bus.str[:2] - - grouping_years = snakemake.params.existing_capacities["grouping_years_industry"] - gdf = gdf[gdf["Start date"] < grouping_years[-1]] - avg_age = gdf["Start date"].mean() - gdf.fillna({"Start date": avg_age}, inplace=True) - - gdf["grouping_year"] = 0 - indices = np.searchsorted(grouping_years, gdf["Start date"], side="right") - gdf.loc[:, "grouping_year"] = np.array(grouping_years)[indices] - gdf["grouping_year"].replace(2025, 2020, inplace=True) - - return gdf - - -def prepare_plant_data( - regions_fn: str, - isi_database: str, -) -> tuple[pd.DataFrame, gpd.GeoDataFrame]: - """ - Reads in the Fraunhofer ISI database with high resolution plant data and maps them to the bus regions. - Returns the database as df as well as the regions as gdf. - - Parameters - ---------- - regions_fn : str - path to the onshore regions file - isi_database: str - path to the fraunhofer isi database - """ - # add existing industry - regions = gpd.read_file(regions_fn).set_index("name") - - isi_data = pd.read_excel(isi_database, sheet_name="Database", index_col=1) - # assign bus region to each plant - geometry = gpd.points_from_xy(isi_data["Longitude"], isi_data["Latitude"]) - plant_data = gpd.GeoDataFrame(isi_data, geometry=geometry, crs="EPSG:4326") - plant_data = gpd.sjoin(plant_data, regions, how="inner", predicate="within") - plant_data.rename(columns={"name": "bus"}, inplace=True) - # filter for countries in model scope - plant_data = plant_data[plant_data.Country.isin(snakemake.params.countries)] - # replace UK with GB in Country column - plant_data["Country"] = plant_data["Country"].replace("UK", "GB") + plant_data["build_year"] = plant_data["build_year"].fillna(mean_nonzero) # assign industry grouping year grouping_years = snakemake.params.existing_capacities["grouping_years_industry"] - plant_data.loc[:, "Year of last modernisation"] = plant_data[ - "Year of last modernisation" - ].replace("x", np.nan) plant_data["grouping_year"] = 0 - valid_mask = plant_data["Year of last modernisation"].notna() - valid_years = plant_data.loc[valid_mask, "Year of last modernisation"] + plant_data["Out"] = plant_data["Out"].fillna(0) + valid_mask = plant_data["build_year"].notna() + valid_years = plant_data.loc[valid_mask, "build_year"] indices = np.searchsorted(grouping_years, valid_years, side="right") plant_data.loc[valid_years.index, "grouping_year"] = np.array(grouping_years)[ indices ] - return plant_data, regions - - -def add_existing_cement_plants(n): - database = prepare_gem_database() - # get plants that are still operating - database = database[database["Operating status"] == "operating"] - - cement = database.groupby(["bus", "country", "grouping_year"], as_index=False)[ - "Cement Capacity (millions metric tonnes per annum)" - ].sum() - - cement.index = cement["bus"] + " cement kiln-" + cement["grouping_year"].astype(str) - cement["clinker"] = ( - cement["bus"] + " cement production-" + cement["grouping_year"].astype(str) - ) - # get rid of capacities without demand - cement_buses = n.loads[n.loads.carrier == "cement"].bus.to_list() - cement_buses = [x.replace(" cement", "") for x in cement_buses] - cement = cement[cement.bus.isin(cement_buses)] - - logger.info(f"Adding {len(cement)} existing cement links.") - - # clinker production - gas_input = ( - costs.at["cement dry clinker", "gas-input"] - + costs.at["cement dry clinker", "heat-input"] - ) - electricity_input = costs.at["cement dry clinker", "electricity-input"] - co2_emission = 0.79 / gas_input + costs.at["gas", "CO2 intensity"] - - n.add( - "Link", - cement.index, - bus0=[bus + " gas" for bus in cement.bus] - if snakemake.params.sector["gas_network"] - else "EU gas", - bus1=[bus + " clinker" for bus in cement.bus], - bus2=cement.bus, - bus3=[bus + " cement emission" for bus in cement.bus], - carrier="cement kiln", - p_nom_extendable=False, - p_nom=cement["Cement Capacity (millions metric tonnes per annum)"] - .mul(gas_input) - .mul(costs.at["cement finishing", "clinker-input"]) - .div(8760) - .mul(1e6) - .values, - p_min_pu=0, - capital_cost=costs.at["cement dry clinker", "capital_cost"] / gas_input, - efficiency=1 / gas_input, - efficiency2=-electricity_input, - efficiency3=co2_emission, - build_year=cement.grouping_year, - lifetime=costs.at["cement dry clinker", "lifetime"], - ) - - # cement finishing - electricity_input = ( - costs.at["cement finishing", "electricity-input"] - / costs.at["cement finishing", "clinker-input"] - ) - clinker_input = costs.at["cement finishing", "clinker-input"] - n.add( - "Link", - cement.clinker.to_list(), - bus0=[bus + " clinker" for bus in cement.bus], - bus1=[bus + " cement" for bus in cement.bus], - bus2=cement.bus.to_list(), - carrier="cement finishing", - p_nom_extendable=False, - p_nom=cement["Cement Capacity (millions metric tonnes per annum)"] - .mul(clinker_input) - .div(8760) - .mul(1e6) - .values, - p_min_pu=0, - capital_cost=costs.at["cement finishing", "capital_cost"] / gas_input, - efficiency=1 / clinker_input, - efficiency2=-electricity_input, - build_year=cement.grouping_year.to_list(), - lifetime=costs.at["cement finishing", "lifetime"], - ) - - -def add_existing_steel_plants( - n: pypsa.Network, -) -> None: - """ - Adds existing steel plants. - The plants are running on natural gas only since the retrofitting to hydrogen would be associated with costs. Plants are not expected to run at a minimal part load to avoid forcing the use of natural gas in planning horizons with climate targets. - Exhaust heat is not integrated since assuming that heat is integrated to make the current process more efficient. - """ - logger.info("Adding existing steel plants.") - - plant_data, regions = prepare_plant_data( - snakemake.input.regions_onshore, - snakemake.input.isi_database, - ) - steel_buses = n.loads[n.loads.carrier == "steel"].bus.to_list() - steel_buses = [x.replace(" steel", "") for x in steel_buses] - - fh_drg = plant_data[plant_data["Process status qup"] == "Direct reduction NG"] - drg = fh_drg.groupby( - ["bus", "Country", "grouping_year", "Product"], as_index=False - ).agg( - { - "Production in tons (calibrated)": "sum", - "Out": "mean", - "Year of last modernisation": "mean", - } - ) - - fh_bof = plant_data[plant_data["Process status qup"] == "Blast furnace"] - bof = fh_bof.groupby( - ["bus", "Country", "grouping_year", "Product"], as_index=False - ).agg( - { - "Production in tons (calibrated)": "sum", - "Out": "mean", - "Year of last modernisation": "mean", - } - ) - # fill Year of last modernisation NaNs with mean - mean_year = int(bof["Year of last modernisation"].mean()) - bof.loc[:, "Year of last modernisation"] = bof["Year of last modernisation"].fillna( - mean_year - ) - drg.loc[:, "Year of last modernisation"] = drg["Year of last modernisation"].fillna( - 2025 - ) - # fill in lifetime - bof.loc[:, "Lifetime"] = ( - bof.loc[:, "Out"] - bof.loc[:, "Year of last modernisation"] - ) - bof.loc[:, "Lifetime"] = bof["Lifetime"].fillna( - costs.at["blast furnace-basic oxygen furnace", "lifetime"] - ) - drg.loc[:, "Lifetime"] = ( - drg.loc[:, "Out"] - drg.loc[:, "Year of last modernisation"] - ) - drg.loc[:, "Lifetime"] = drg["Lifetime"].fillna( - costs.at["natural gas direct iron reduction furnace", "lifetime"] - ) - - drg.index = drg["bus"] + " gas DRI-" + drg["grouping_year"].astype(str) - drg = drg[drg.bus.isin(steel_buses)] - bof.index = bof["bus"] + " BOF-" + bof["grouping_year"].astype(str) - bof = bof[bof.bus.isin(steel_buses)] - - # add direct reduction with natural gas - gas_input = costs.at["natural gas direct iron reduction furnace", "gas-input"] - marginal_cost = ( - costs.at["iron ore DRI-ready", "commodity"] - * costs.at["natural gas direct iron reduction furnace", "ore-input"] - / gas_input - ) - - logger.info(f"Adding {len(drg)} gas DRI plant.") - - n.add( - "Link", - drg.index, - bus0=[bus + " gas" for bus in drg.bus] - if snakemake.params.sector["gas_network"] - else "EU gas", - bus1=[bus + " hbi" for bus in drg.bus] - if not snakemake.params.sector["industry_relocation"] - else "EU hbi", - bus2=[bus + " gas DRI emission" for bus in drg.bus], - p_nom=drg["Production in tons (calibrated)"].mul(gas_input).div(8760).values, - p_nom_extendable=False, - carrier="gas DRI", - efficiency=1 / gas_input, - efficiency2=costs.at["gas", "CO2 intensity"], - capital_cost=costs.at[ - "natural gas direct iron reduction furnace", "capital_cost" - ] - / gas_input, - marginal_cost=marginal_cost, - build_year=drg["Year of last modernisation"], - lifetime=drg["Lifetime"], - ) - - # add coal blast furnaces - coal_input = costs.at["blast furnace-basic oxygen furnace", "coal-input"] - marginal_cost = ( - costs.at["iron ore DRI-ready", "commodity"] - * costs.at["blast furnace-basic oxygen furnace", "ore-input"] - / coal_input - ) - logger.info(f"Adding {len(bof)} BOF plant.") - - n.add( - "Link", - bof.index, - bus0="EU coal", - bus1=[bus + " steel" for bus in bof.bus], - bus2=[bus + " BOF emission" for bus in bof.bus], - p_nom=bof["Production in tons (calibrated)"].mul(coal_input).div(8760).values, - p_nom_extendable=False, - carrier="BOF", - efficiency=1 / coal_input, - efficiency2=costs.at["coal", "CO2 intensity"], - marginal_cost=marginal_cost, - capital_cost=costs.at["blast furnace-basic oxygen furnace", "capital_cost"] - / coal_input, - build_year=bof["Year of last modernisation"], - lifetime=bof["Lifetime"], - ) - - -def add_existing_meoh_plants(n): - logger.info("Adding existing methanol plants.") - - plant_data, regions = prepare_plant_data( - snakemake.input.regions_onshore, - snakemake.input.isi_database, - ) + plant_data = plant_data.groupby( + ["bus", "country", "carrier", "grouping_year", "Out"], as_index=False + )["p_set"].sum() + + if "cement" in options["endogenous_sectors"]: + # add cement + cement = plant_data[plant_data.carrier=="cement"] + index = [f"{bus} cement klin-{year}" for bus, year in zip(cement["bus"], cement["grouping_year"])] + cement.index = index + logger.info(f"Adding {len(cement)} existing cement links.") + + # clinker production + gas_input = ( + costs.at["cement dry clinker", "gas-input"] + + costs.at["cement dry clinker", "heat-input"] + ) + electricity_input = costs.at["cement dry clinker", "electricity-input"] + co2_emission = 0.79 / gas_input + costs.at["gas", "CO2 intensity"] + + n.add( + "Link", + index, + bus0=[bus + " gas" for bus in cement.bus] + if snakemake.params.sector["gas_network"] + else "EU gas", + bus1=[bus + " clinker" for bus in cement.bus], + bus2=cement.bus, + bus3=[bus + " cement emission" for bus in cement.bus], + carrier="cement kiln", + p_nom_extendable=False, + p_nom=cement["p_set"] + .mul(gas_input) + .mul(costs.at["cement finishing", "clinker-input"]) + .div(8760) + .values, + p_min_pu=0, + capital_cost=costs.at["cement dry clinker", "capital_cost"] / gas_input, + efficiency=1 / gas_input, + efficiency2=-electricity_input, + efficiency3=co2_emission, + build_year=cement.grouping_year, + lifetime=costs.at["cement dry clinker", "lifetime"], + ) - fh_meoh = plant_data[plant_data.Product == "Methanol"] - fh_meoh["grouping_year"].replace(2025, 2020, inplace=True) - fh_meoh = fh_meoh.groupby( - ["bus", "Country", "grouping_year", "Product"], as_index=False - )["Production in tons (calibrated)"].sum() + # cement finishing + index = [f"{bus} cement production-{year}" for bus, year in zip(cement["bus"], cement["grouping_year"])] + cement.index = index + electricity_input = ( + costs.at["cement finishing", "electricity-input"] + / costs.at["cement finishing", "clinker-input"] + ) + clinker_input = costs.at["cement finishing", "clinker-input"] + n.add( + "Link", + index, + bus0=[bus + " clinker" for bus in cement.bus], + bus1=[bus + " cement" for bus in cement.bus], + bus2=cement.bus, + carrier="cement finishing", + p_nom_extendable=False, + p_nom=cement["p_set"] + .mul(clinker_input) + .div(8760) + .values, + p_min_pu=0, + capital_cost=costs.at["cement finishing", "capital_cost"] / gas_input, + efficiency=1 / clinker_input, + efficiency2=-electricity_input, + build_year=cement.grouping_year, + lifetime=costs.at["cement finishing", "lifetime"], + ) - fh_meoh.index = ( - fh_meoh["bus"] + " grey methanol-" + fh_meoh["grouping_year"].astype(str) - ) + if options["ammonia"]: + ammonia = plant_data[plant_data.carrier=="Haber-Bosch"] + index = [f"{bus} Haber-Bosch-{year}" for bus, year in zip(ammonia["bus"], ammonia["grouping_year"])] + ammonia.index = index + logger.info(f"Adding {len(ammonia)} existing Haber-Bosch links.") + # https://dechema.de/dechema_media/Downloads/Positionspapiere/Technology_study_Low_carbon_energy_and_feedstock_for_the_European_chemical_industry.pdf + # page 56: 1.83 t_CO2/t_NH3 + ch4_per_nh3 = ( + 1.83 / costs.at["gas", "CO2 intensity"] / snakemake.params["MWh_NH3_per_tNH3"] + ) + n.add( + "Link", + index, + bus0=[bus + " gas" for bus in ammonia.bus] + if snakemake.params.sector["gas_network"] + else "EU gas", + bus1=[bus + " NH3" for bus in ammonia.bus] + if snakemake.params.sector["ammonia"] + else "EU NH3", + bus2=ammonia.bus, + bus3="co2 atmosphere", + p_nom=ammonia["p_set"] + .mul(snakemake.params.MWh_NH3_per_tNH3) + .div(ch4_per_nh3) + .div(8760) + .values, + p_nom_extendable=False, + carrier="Haber-Bosch", + efficiency=1 / ch4_per_nh3, + efficiency1=-costs.at["Haber-Bosch", "electricity-input"] / ch4_per_nh3, + efficiency2=costs.at["gas", "CO2 intensity"], + capital_cost=costs.at["Haber-Bosch", "capital_cost"] + / costs.at["Haber-Bosch", "electricity-input"], + marginal_cost=costs.at["Haber-Bosch", "VOM"] + / costs.at["Haber-Bosch", "electricity-input"], + build_year=ammonia["grouping_year"], + lifetime=costs.at["Haber-Bosch", "lifetime"], + ) - # grey methanol efficiency hard coded for now - costs.at["grey methanol synthesis", "efficiency"] = 1 / 1.757469244 + # methanol + meoh = plant_data[plant_data.carrier=="grey methanol"] + index = [f"{bus} grey methanol-{year}" for bus, year in zip(meoh["bus"], meoh["grouping_year"])] + meoh.index = index + logger.info(f"Adding {len(meoh)} existing methanolisation links.") capital_cost = ( - costs.at["SMR", "capital_cost"] - + costs.at["methanolisation", "capital_cost"] - * costs.at["grey methanol synthesis", "efficiency"] + costs.at["grey methanol synthesis", "investment"] + / costs.at["grey methanol synthesis", "efficiency"] ) co2_emissions = ( - costs.at["gas", "CO2 intensity"] - - costs.at["grey methanol synthesis", "efficiency"] - * costs.at["methanol", "CO2 intensity"] + costs.at["grey methanol synthesis", "carbondioxide-output"] ) n.add( "Link", - fh_meoh.index, - bus0=[bus + " gas" for bus in fh_meoh.bus] + index, + bus0=[bus + " gas" for bus in meoh.bus] if snakemake.params.sector["gas_network"] else "EU gas", bus1="EU methanol", bus2="co2 atmosphere", p_nom_extendable=False, - p_nom=fh_meoh["Production in tons (calibrated)"] + p_nom=meoh["p_set"] .mul(snakemake.params.MWh_MeOH_per_tMeOH) + .div(8760) .div(costs.at["grey methanol synthesis", "efficiency"]), carrier="grey methanol", efficiency=costs.at["grey methanol synthesis", "efficiency"], efficiency2=co2_emissions, capital_cost=capital_cost, - build_year=fh_meoh.grouping_year, + build_year=meoh.grouping_year, lifetime=costs.at["SMR", "lifetime"], ) + if "steel" in options["endogenous_sectors"]: + # natural gas DRI + ng_dri = plant_data[plant_data.carrier=="gas DRI"] + logger.info(f"Adding {len(ng_dri)} existing gas DRI links.") + index = [f"{bus} gas DRI-{year}" for bus, year in zip(ng_dri["bus"], ng_dri["grouping_year"])] + ng_dri.index = index + gas_input = costs.at["natural gas direct iron reduction furnace", "gas-input"] + marginal_cost = ( + costs.at["iron ore DRI-ready", "commodity"] + * costs.at["natural gas direct iron reduction furnace", "ore-input"] + / gas_input + ) + n.add( + "Link", + index, + bus0=[bus + " gas" for bus in ng_dri.bus] + if snakemake.params.sector["gas_network"] + else "EU gas", + bus1=[bus + " hbi" for bus in ng_dri.bus] + if not snakemake.params.sector["hbi_relocation"] + else "EU hbi", + bus2=[bus + " gas DRI emission" for bus in ng_dri.bus], + p_nom=ng_dri["p_set"].mul(gas_input).div(8760).values, + p_nom_extendable=False, + carrier="gas DRI", + efficiency=1 / gas_input, + efficiency2=costs.at["gas", "CO2 intensity"], + capital_cost=costs.at[ + "natural gas direct iron reduction furnace", "capital_cost" + ] + / gas_input, + marginal_cost=marginal_cost, + build_year=ng_dri.grouping_year, + lifetime=costs.at["natural gas direct iron reduction furnace", "lifetime"], + ) + # BF-BOF + bof = plant_data[plant_data.carrier=="BOF"] -def add_existing_ammonia_plants( - n: pypsa.Network, -) -> None: - """ - Adds existing Haber-Bosch plants. - The plants are running on natural gas only since the retrofitting to hydrogen would be associated with costs. Plants are not expected to run at a minimal part load to avoid forcing the use of natural gas in planning horizons with climate targets. - Exhaust heat is not integrated since assuming that heat is integrated to make the current process more efficient. - """ - logger.info("Adding existing Haber-Bosch plants.") - - plant_data, regions = prepare_plant_data( - snakemake.input.regions_onshore, - snakemake.input.isi_database, - ) - - fh_ammonia = plant_data[plant_data.Product == "Ammonia"] - - fh_ammonia = fh_ammonia.groupby( - ["bus", "Country", "grouping_year", "Product"], as_index=False - )["Production in tons (calibrated)"].sum() - - fh_ammonia.index = ( - fh_ammonia["bus"] - + " Haber-Bosch-SMR-" - + fh_ammonia["grouping_year"].astype(str) - ) - # add dataset for Non EU27 countries - df = pd.read_csv(snakemake.input.ammonia, index_col=0) - - geometry = gpd.points_from_xy(df.Longitude, df.Latitude) - gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326") - - gdf = gpd.sjoin(gdf, regions, how="inner", predicate="within") - - gdf.rename(columns={"name": "bus"}, inplace=True) - gdf["Country"] = gdf.bus.str[:2] - # filter for countries that are missing - gdf = gdf[ - (~gdf.Country.isin(fh_ammonia.Country.unique())) - & (gdf.Country.isin(snakemake.params.countries)) - ] - # following approach from build_industrial_distribution_key.py - for country in gdf.Country: - facilities = gdf.query("Country == @country") - production = facilities["Ammonia [kt/a]"] - # assume 50% of the minimum production for missing values - production = production.fillna(0.5 * facilities["Ammonia [kt/a]"].min()) - - # missing data - gdf.drop(gdf[gdf["Ammonia [kt/a]"].isna()].index, inplace=True) - - # get average plant age: - avg_age = plant_data[plant_data.Product == "Ammonia"][ - "Year of last modernisation" - ].mean() - gdf["grouping_year"] = min( - y - for y in snakemake.params.existing_capacities["grouping_years_industry"] - if y > avg_age - ) - # match database - gdf.index = ( - gdf["bus"] + " Haber-Bosch-SMR-" + gdf["grouping_year"].values.astype(str) - ) - gdf.rename( - columns={"Ammonia [kt/a]": "Production in tons (calibrated)"}, inplace=True - ) - gdf["Production in tons (calibrated)"] *= 1e3 + if options["steel_bof"]["pledge"]: + lifetime = bof["Out"] - bof["grouping_year"] + options["steel_bof"]["pledge_delay"] + if options["steel_bof"]["default_phase_out"]: + lifetime[lifetime < 0] = options["steel_bof"]["default_phase_out"] + else: + lifetime[lifetime < 0] = costs.at["blast furnace-basic oxygen furnace", "lifetime"] + else: + lifetime = costs.at["blast furnace-basic oxygen furnace", "lifetime"] + # regroup + bof = bof.groupby( + ["bus", "country", "carrier", "grouping_year"], as_index=False + )["p_set"].sum() + + index = [f"{bus} BOF-{year}-{out}" for bus, year, out in zip(bof["bus"], bof["grouping_year"], bof["Out"])] + bof.index = index + logger.info(f"Adding {len(bof)} existing BOF links.") + + coal_input = costs.at["blast furnace-basic oxygen furnace", "coal-input"] + marginal_cost = ( + costs.at["iron ore DRI-ready", "commodity"] + * costs.at["blast furnace-basic oxygen furnace", "ore-input"] + / coal_input + ) - ammonia_plants = pd.concat( - [ - fh_ammonia, - gdf[["bus", "Country", "grouping_year", "Production in tons (calibrated)"]], - ] - ) + n.add( + "Link", + bof.index, + bus0="EU coal", + bus1=[bus + " steel" for bus in bof.bus], + bus2=[bus + " BOF emission" for bus in bof.bus], + p_nom=bof["p_set"].mul(coal_input).div(8760).values, + p_nom_extendable=False, + carrier="BOF", + efficiency=1 / coal_input, + efficiency2=costs.at["coal", "CO2 intensity"], + marginal_cost=marginal_cost, + capital_cost=costs.at["blast furnace-basic oxygen furnace", "capital_cost"] + / coal_input, + build_year=bof["grouping_year"], + lifetime=lifetime.to_list(), + ) - # https://dechema.de/dechema_media/Downloads/Positionspapiere/Technology_study_Low_carbon_energy_and_feedstock_for_the_European_chemical_industry.pdf - # page 56: 1.83 t_CO2/t_NH3 - ch4_per_nh3 = ( - 1.83 / costs.at["gas", "CO2 intensity"] / snakemake.params["MWh_NH3_per_tNH3"] - ) - n.add( - "Link", - ammonia_plants.index, - bus0=[bus + " gas" for bus in ammonia_plants.bus] - if snakemake.params.sector["gas_network"] - else "EU gas", - bus1=[bus + " NH3" for bus in ammonia_plants.bus] - if snakemake.params.sector["ammonia"] - else "EU NH3", - bus2=ammonia_plants.bus, - bus3="co2 atmosphere", - p_nom=ammonia_plants["Production in tons (calibrated)"] - .mul(snakemake.params.MWh_NH3_per_tNH3) - .div(ch4_per_nh3) - .div(8760) - .values, - p_nom_extendable=False, - carrier="Haber-Bosch", - efficiency=1 / ch4_per_nh3, - efficiency1=-costs.at["Haber-Bosch", "electricity-input"] / ch4_per_nh3, - efficiency2=costs.at["gas", "CO2 intensity"], - capital_cost=costs.at["Haber-Bosch", "capital_cost"] - / costs.at["Haber-Bosch", "electricity-input"], - marginal_cost=costs.at["Haber-Bosch", "VOM"] - / costs.at["Haber-Bosch", "electricity-input"], - build_year=ammonia_plants["grouping_year"], - lifetime=costs.at["Haber-Bosch", "lifetime"], - ) if __name__ == "__main__": @@ -1236,8 +971,8 @@ def add_existing_ammonia_plants( snakemake = mock_snakemake( "add_existing_baseyear", - configfiles="config/test/config.myopic.yaml", - clusters="5", + configfiles="config/config.default.yaml", + clusters="39", opts="", sector_opts="", planning_horizons=2030, @@ -1312,16 +1047,8 @@ def add_existing_ammonia_plants( cluster_heat_buses(n) # add existing industry plants - if "steel" in snakemake.params.sector["endogenous_sectors"]: - add_existing_steel_plants(n) - if "cement" in snakemake.params.sector["endogenous_sectors"]: - add_existing_cement_plants(n) - if ("ammonia" in snakemake.params.sector["endogenous_sectors"]) and ( - snakemake.params.sector["ammonia"] - ): - add_existing_ammonia_plants(n) - if "methanol" in snakemake.params.sector["endogenous_sectors"]: - add_existing_meoh_plants(n) + if options["industry"]: + add_existing_industry(n) n.meta = dict(snakemake.config, **dict(wildcards=dict(snakemake.wildcards))) diff --git a/scripts/build_industry_plants.py b/scripts/build_industry_plants.py new file mode 100644 index 0000000000..cba6f8fec6 --- /dev/null +++ b/scripts/build_industry_plants.py @@ -0,0 +1,199 @@ +# SPDX-FileCopyrightText: Contributors to PyPSA-Eur +# +# SPDX-License-Identifier: MIT +""" +Adds existing power and heat generation capacities for initial planning +horizon. +""" + +import logging +import re +from types import SimpleNamespace + +import country_converter as coco +import geopandas as gpd +import numpy as np +import pandas as pd +import pycountry + +from scripts._helpers import ( + configure_logging, + set_scenario_config, +) + +logger = logging.getLogger(__name__) +cc = coco.CountryConverter() +idx = pd.IndexSlice +spatial = SimpleNamespace() + + +def country_to_code(name): + try: + return pycountry.countries.lookup(name).alpha_2 + except LookupError: + return None + + +def prepare_gem_database(regions): + """ + Load GEM database of cement plants and map onto bus regions. + """ + df = pd.read_excel( + snakemake.input.gem_gcct, + sheet_name="Plant Data", + na_values=["N/A", "unknown", ">0"], + ) + + df["country_code"] = df["Country/Area"].apply(country_to_code) + + df = df[df.country_code.isin(snakemake.params.countries)] + + df["Start date"] = pd.to_numeric( + df["Start date"].str.split("-").str[0], errors="coerce" + ) + + latlon = ( + df["Coordinates"] + .str.split(", ", expand=True) + .rename(columns={0: "lat", 1: "lon"}) + ) + geometry = gpd.points_from_xy(latlon["lon"], latlon["lat"]) + gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326") + + gdf = gpd.sjoin(gdf, regions, how="inner", predicate="within") + + gdf.rename(columns={"name": "bus"}, inplace=True) + gdf["country"] = gdf.bus.str[:2] + + # just get bus, country, grouping_year, plant_size, date_in, date_out, + # get plants that are still operating + gdf = gdf[gdf["Operating status"] == "operating"] + + gdf.rename(columns={ + "Cement Capacity (millions metric tonnes per annum)": "p_set", + "Start date": "build_year", + }, inplace=True) + gdf["p_set"] *= 1e6 + gdf["carrier"] = "cement" + gdf["Out"] = 0 + + cement = gdf[['bus', 'country', 'carrier', 'p_set', 'build_year', 'Out']] + + return cement + + +def prepare_plant_data( + regions: gpd, + isi_database: str, +) -> tuple[pd.DataFrame, gpd.GeoDataFrame]: + """ + Reads in the Fraunhofer ISI database with high resolution plant data and maps them to the bus regions. + Returns the database as df. + + Parameters + ---------- + regions : gpd + onshore regions on which the industry sites are mapped to + isi_database: str + path to the fraunhofer isi database + """ + + isi_data = pd.read_excel(isi_database, sheet_name="Database", index_col=1) + # assign bus region to each plant + geometry = gpd.points_from_xy(isi_data["Longitude"], isi_data["Latitude"]) + plant_data = gpd.GeoDataFrame(isi_data, geometry=geometry, crs="EPSG:4326") + plant_data = gpd.sjoin(plant_data, regions, how="inner", predicate="within") + plant_data.rename(columns={"name": "bus"}, inplace=True) + # filter for countries in model scope + plant_data = plant_data[plant_data.Country.isin(snakemake.params.countries)] + # replace UK with GB in Country column + plant_data["country"] = plant_data["Country"].replace("UK", "GB") + # add carrier column + carrier_dict = { + "Blast furnace": "BOF", + "Direct reduction NG": "gas DRI", + "Ammonia SMR": "Haber-Bosch", + "Methanol SMR": "grey methanol", + } + plant_data["carrier"] = plant_data["Process status qup"].replace(carrier_dict) + # get build_year + plant_data["build_year"] = plant_data[ + "Year of last modernisation"].replace("x", np.nan).fillna(plant_data["Last Relining"]) + plant_data.rename(columns={"Production in tons (calibrated)": "p_set"}, inplace=True) + plant_data = plant_data[['bus', 'country', 'carrier', 'p_set', 'build_year', 'Out']] + + return plant_data + + +def prepare_ammonia_data(regions, plant_data): + + df = pd.read_csv(snakemake.input.ammonia, index_col=0) + + geometry = gpd.points_from_xy(df.Longitude, df.Latitude) + gdf = gpd.GeoDataFrame(df, geometry=geometry, crs="EPSG:4326") + + gdf = gpd.sjoin(gdf, regions, how="inner", predicate="within") + + gdf.rename(columns={"name": "bus"}, inplace=True) + gdf["country"] = gdf.bus.str[:2] + # filter for countries that are missing + gdf = gdf[ + (~gdf.country.isin(plant_data.country.unique())) + & (gdf.country.isin(snakemake.params.countries)) + ] + # following approach from build_industrial_distribution_key.py + for country in gdf.Country: + facilities = gdf.query("country == @country") + production = facilities["Ammonia [kt/a]"] + # assume 50% of the minimum production for missing values + production = production.fillna(0.5 * facilities["Ammonia [kt/a]"].min()) + + # missing data + gdf.drop(gdf[gdf["Ammonia [kt/a]"].isna()].index, inplace=True) + + # get average plant age: + avg_age = plant_data[plant_data.carrier == "Haber-Bosch"][ + "build_year"].mean() + # restructure to fit other data sources + gdf["build_year"] = avg_age + gdf["Out"] = 0 + gdf["carrier"] = "Haber-Bosch" + gdf["Ammonia [kt/a]"] *= 1e3 + gdf.rename(columns={"Ammonia [kt/a]": "p_set"}, inplace=True) + gdf = gdf[['bus', 'country', 'carrier', 'p_set', 'build_year', 'Out']] + + return gdf + + + +if __name__ == "__main__": + if "snakemake" not in globals(): + from scripts._helpers import mock_snakemake + + snakemake = mock_snakemake( + "build_industry_plants", + configfiles="config/config.default.yaml", + clusters="39", + ) + + configure_logging(snakemake) + set_scenario_config(snakemake) + + regions = gpd.read_file(snakemake.input.regions_onshore).set_index("name") + + # get cement plant data + cement = prepare_gem_database(regions) + # get data from fraunhofer database + fh_data = prepare_plant_data( + regions, + snakemake.input.isi_database + ) + # ammonia plants non-EU27 + ammonia = prepare_ammonia_data( + regions, + fh_data + ) + + industry_plants = pd.concat([cement, fh_data, ammonia], ignore_index=True) + + industry_plants.to_csv(snakemake.output.industry_plants) diff --git a/scripts/prepare_sector_network.py b/scripts/prepare_sector_network.py index cf7b0abd6e..1f389c0a7a 100755 --- a/scripts/prepare_sector_network.py +++ b/scripts/prepare_sector_network.py @@ -222,7 +222,7 @@ def define_spatial(nodes, options): # industry spatial.hbi = SimpleNamespace() - if options["industry_relocation"]: + if options["hbi_relocation"]: spatial.hbi.nodes = ["EU hbi"] spatial.hbi.location = ["EU"] else: @@ -4588,12 +4588,15 @@ def add_industry( industrial_demand_raw = ( pd.read_csv(industrial_demand_file, index_col=0) * 1e6 * nyears ) - # adjust the industrial energy demand excluding endogenously modelled carriers - adjusted_demand, steel, cement = adjust_industry_demand(nodes) - adjusted_demand["current electricity"] = industrial_demand_raw[ - "current electricity" - ] - industrial_demand = adjusted_demand.copy() + if snakemake.params.sector["endogenous_sectors"]: + # adjust the industrial energy demand excluding endogenously modelled carriers + adjusted_demand, steel, cement = adjust_industry_demand(nodes) + adjusted_demand["current electricity"] = industrial_demand_raw[ + "current electricity" + ] + industrial_demand = adjusted_demand.copy() + else: + industrial_demand = industrial_demand_raw.copy() n.add( "Bus", @@ -5539,17 +5542,13 @@ def add_industry( p_nom_extendable=True, ) - costs.at["grey methanol synthesis", "efficiency"] = 1 / 1.757469244 # grey methanol capital_cost = ( - costs.at["SMR", "investment"] - + costs.at["methanolisation", "investment"] - * costs.at["grey methanol synthesis", "efficiency"] + costs.at["grey methanol synthesis", "investment"] + / costs.at["grey methanol synthesis", "efficiency"] ) co2_emissions = ( - costs.at["gas", "CO2 intensity"] - - costs.at["grey methanol synthesis", "efficiency"] - * costs.at["methanol", "CO2 intensity"] + costs.at["grey methanol synthesis", "carbondioxide-output"] ) n.add( "Link", @@ -5566,38 +5565,6 @@ def add_industry( lifetime=costs.at["SMR", "lifetime"], ) - # blue methanol - options["cc_fraction"] = 0.9 - co2_emissions = ( - costs.at["gas", "CO2 intensity"] - - costs.at["grey methanol synthesis", "efficiency"] - * costs.at["methanol", "CO2 intensity"] - ) - capital_cost = ( - costs.at["SMR", "investment"] - + costs.at["methanolisation", "investment"] - * costs.at["grey methanol synthesis", "efficiency"] - + costs.at["cement capture", "investment"] - * co2_emissions - * options["cc_fraction"] - ) - n.add( - "Link", - spatial.methanol.demand_locations, - suffix=" blue methanol", - bus0=spatial.gas.nodes, - bus1=spatial.methanol.nodes, - bus2="co2 atmosphere", - bus3=spatial.co2.nodes, - p_nom_extendable=True, - carrier="blue methanol", - efficiency=costs.at["grey methanol synthesis", "efficiency"], - efficiency2=co2_emissions * (1 - options["cc_fraction"]), - efficiency3=co2_emissions * options["cc_fraction"], - capital_cost=capital_cost, - lifetime=costs.at["SMR", "lifetime"], - ) - def add_aviation( n: pypsa.Network, From 321cb4b19094c60e987d5aeaf9ad0f12488e859e Mon Sep 17 00:00:00 2001 From: toniseibold Date: Fri, 5 Dec 2025 06:37:16 +0100 Subject: [PATCH 08/29] adding data files --- data/1-s2.0-S0196890424010586-mmc2.xlsx | Bin 0 -> 40060 bytes ...-Cement-and-Concrete-Tracker_July-2025.xlsx | Bin 0 -> 972391 bytes 2 files changed, 0 insertions(+), 0 deletions(-) create mode 100644 data/1-s2.0-S0196890424010586-mmc2.xlsx create mode 100644 data/gem/Global-Cement-and-Concrete-Tracker_July-2025.xlsx diff --git a/data/1-s2.0-S0196890424010586-mmc2.xlsx b/data/1-s2.0-S0196890424010586-mmc2.xlsx new file mode 100644 index 0000000000000000000000000000000000000000..0b613ceea356b3e522792b3ed3783b7af2692e84 GIT binary patch literal 40060 zcmbTcV{~P0(=OWS*ekYe+h)hMZQJVDMkiUZZQC8&>bPSkJI}MvH_jR7{l5KUuOD-) ziK*P~J|F89ZOQVb5&k*x2|m49yi6wy$qC`XjyWyxhRx69bV zb=tiU#HltB?UxhejsIz+V@G^*VY!HpG^0U{g(fX58f!2*b-w9RbH}b`b(Y`5iZoV9 zN8W9p2(^W7V7;lM<~R&EQEpcu6wY-pW|d1B@)@YPlHm{Lugu6r6WRH z_yQq?sPiAgup(yj$*>;PC^r_0KE+&fdI|{uf4AayRd1#r&8Btz>qt^%i~&E*tUk#; z5Pl9tmjVWj7lb(xUZzrLpiHmH!-PzWJdSEP>>s3EBHX2K38b~Cmwu`Viw>UNpu+#u z%6*%vuZFo7HJ+zLWYA%5ii2wYUFcWX0opxgBPzbzD#I(YKT@uM zga!5|)ycLJc29xQtR6Hv98)#*LrlEdd}=wV`WK>=^0nv?UxiLNQC;b_4QDgMWENcp&rP%HVZ z7l8Sya}QJFxBa9WQ!@vM+Uc#!B=*_sQYV_6)|nb7rtCW`Oal${iQym^MGOfmd5c2P zpi#K$JsU6&RpgwDnQL`sQ2quvq-Z5u4ggk7;ZTb>+A7MZzI*ZOy>6MHQN5)CZFtD z?cX{iqpTnlxJiGHTw`D;+axN9LNokQTv^&^tyu=f?pE%UO#NoI2or0Ti$M~Y5pv<@ zI{m>sftim;=82?|DgLtn7e_`HqjC-)NlQCYg8;ScotVd2<^D%Yu+@t=_8BUR|Mfmj z?fD^>=5ORPjMrt5g>MwXH@o51NnN^XaNtY>52)2|s2MxqB) z!54{(3_Tt5Fg*R=Lw?!a+Y6#@c;R{8^R!Y+<;$-C62-`Wk!m-Q;0wi z_xES`JR{Z)Ug?f3e;f(s5XfrY!?!W%4^QOoh<@6ProRZGUK#V}X-6vf&~FVYIIs?V zO(#@=tE#L6gL{`u=GsHOBonZpLZPN*RV!JZI1IHjF9@5IjTMveM#13wJg|WQs2Y}5 zHaSaR*NJ<%N8+M)2~c=Q$g&>^n@BZ(v^u$|ORtoU4QBHgr}HLb=Or#I&Jr73j6?1A z48p&uh2ja9T-VM`zGOLtcGTG|9$~8#RmYX)WTo=jWlwDoD24AqrW{P7|=z_pW z9kkn1I%Ump9EqU}C+_eXce_JdJAd+7mmdGn=KQv7pFn9t3;Khx#PVWWjjCjuk}05r zl=-r3!+AaEw<3@pjV)_1{b&7C1F2g^Uz z<8yGjzvI;NcO8ah>guoTco*f4CCUt)E?&$ZOzTHkQ7w1QtZ?Ie+^R?z-pY5;;L2ic z=>)*}i<2uJ#khW-#31*ee(-H-KIiKWtD4QQ*m(r_ZkFGFuKM-rKj8l<9Jgttfv2y+ zp@9EC70%bqf4aG|i>HmL^FJkXp`#eL#tHE2(MLM0n_~+UBt%oV=X7jOQ(m+yN|N(s z{26G0hr?A2`+jdeLPQ$B`hA3cY*F(fY|_Jw-yfgHEyDfIFAH#nb7VRBXXUcf)fRR> zUW`dM_JRSiK|j1X1xA!a6m!VfJ6!qo{UA+|+6~NK&ybe(B(T`nn{4nUXCJr=LqX{~ zv6L=76?6>-?Kzev*M5X4R>mwHG4HCoEu|oPF&1H$s7b6uV}e6H?A8L%e7Q3xt$XHM zpSvukO;q^kMRgXnDQR}-(OOX*%I2LmX*9S4*jZTjnI!wQbQ2Lz@nonc8~E~xk>=iN zKx0t;9FvP;DcOa&R+;++-@v%Gnjq71dIK8A_J9%&&zxr+x`R{cf-yJG7mQCRyxgN1 zg7K|uHd4h6SxH286ly;;{{34Xv%l)fIdG1h94M(hvlKLD5fjpzJ$(?G;_tUoW$Hlv zTIrUQ@+x#fCeb0+;@>SaE{RKpI$_y1t{HRO7R_*M?461?Z+X;yPz{xwV*-90E#btt zrmlF?mZqr~of=JiW6IZz5dmDG7wu`FXvdHf2V)3{9@syW&H>k^;gp3+9YyrU2^ID3 zrVmPq2{7jAQ|y#Fx9U|soZ8=}+p#QN{`pC}TcfwDUTH|OF9PpJ?Z5fp-}dAV`{eVQ zhavDuD6s-SxzWsc>*AbjSzzqY&6aDna3eYnvGSM}u#$fU*Qyz_8%|C-N3RC&F%_dO zKqS?h10z6t8_np=A3e&=@L(2x4O)u4e(G?kvJ0E=sV*rlCq_t-RNrm}+lrKzCA<(b z^4ha_CkzE2Vv4EnRZ6S&kR>2kp`~9)EOrGMAg1*hoW(z&js;ql0ULC2nGW$iPIV@f zV{g+->su!m_@&lBi>-UV8B||tLqvlZ_hK?aOtKkx=_OqzcgYq0zQYbJUMFQ0O6pZ1 zKMBSe05Lo9W_?UAHizRp1drC6PF7>ib-$1{>)1w94CR}T8DZ2POVLP;E(s$3eY)P-_kiFOL zUEYy+USzv{bGAuY?ybMAoKVI!xQa@`kK8;XCX(e zjFQYnG~T*@S=AbFaNP$51xeZWFdJJ|Gx3-2vs{nSP9VgRz~yrPXbt9blQ_;a6TExq z)GwoX@9AP>YmGXbg;n1{9p7t(bTG6eznrlwcwgZ8&u%JtYZ>u@NoyYI4E0nqXXa7b9Lmm(Mq!*6yvct1X^UQ=ivlhb zh)LnWeI}6x!7^SDc2UJ6#44N+4q45ztnNMNP+1G~fujCA;a}eEaLh%V9(jN)(U#A) zE0H~bz2(Y#Ud|U6?S#2{cPWhJ;F&jI?XC?ja1C`@w6vIUSz39U{}d|dWhrdWLf03U z(!H@-uF$D`Yxq{BLU`pDKAd&Xx4BZC${EFN$d&6wdEZyTuEE=&UwVZ$r$`?7K*cVk zB|#`u2+4iZN0noXH@AoVcIPK|#w7oLxya7YUnorQZ{L)~|4%KIT;^!@S;KD~bnNr3)n8GjiPNMrWhH?dt%HW+T<4fF zPQoOn*c%N3(Zre^Z1n$`B?-nG1*KnZVT1f%JN_>p``7WRgdRCi zCb;28{vj%4k}+M#kn-1oZ<2W9-wr|M2Hn%x(?sQj*IPzZs1?dUL8DAZ9% zt{fA8O5=*gTTws)ZtFWyp0fJZV9n|l;0KDTHR|{wo?zD>yA=QaQO)fwqCqwyMb_GP z154o)v|1sJ`QcZWXyJ6g5dR1<|71|G%Qt>5z76UGuREDVSGv%72%z<$lH}c*lu3kS zol`o+3$2>@BWF5_W4$2tQQ2~nS;>oM5o(s+dg=L*m&-QJw|P0jbbQ=$Wi|vY)045d zw^w~Z@lSNG?~J(Nzaoq9HLd?=ii`Olu=*EJ|HidzUE8^E6bbpO+{-P^kPMah&$@}VE@lOl3M$2HExu)x)a+iNCJzn;YSIf4UV?!IrF zo7XaL*(^|>l=mZ%CF(grBr*3b`!kbYzQK9??U%=Y0vAbLq+OOC?rCo&g`p+t7*DhhD$r70lK=4L)m`BS*vT&bz zdAY^k4yhdTSzEk4#W()TiH!%!zR+UD2kmw9A0!RD&nIaU&lFn^84{;d(fw~#>3==k zwHrUJxw%6WrZ5uM(`SlqvU+U%y)p;puHN9MmYAG|ldVT5Mg9KdHIGsze+~Yg zHG1!CM@v+NvHI}e5OJT!bNo2^m?w!^;Wrnm$qJXjzt+RwmXnuSXav*-`~P%eg4ohe zT>9t__u+^Y6C~L(b3*ssXrM$!!Hav?(@Ao#)rywbPCU44eI&srw`4ziyF{}uq5gjQ zWIt9EA>*o2Sye*%Of1>RgJMZtO!VF7VTij@oC*I~M_+>vyAAfxTVmC}JbTo3T6BJM z+10IauYC8U+Nf{;uQ%p_1V6K7y$0=ZVyDgbp<5$HBv(HdCM3v*G+KE%w$$0z;UpEA zlpr8n&w!rBAP$xID%+Md*O{fgLWef~Y2_1D5dQOVWzur}l~rR3-T$6?5ZgCvaipEzVav*<0qzVoLr+VA8B2@4hdQv6MYi7E?&4xSYx(pR*T zq9C3P+4{4d*CqMTh;V{^oUC$aQOUNF?XJF z#xj|T5>*8T8$2gSs1InMSWCr;sspq2mP5%97JCyU?kQ9#E0q6B={MycR5%!N@Tj04 zefc{|f6=c*J-J+x(}o1uZIjbj$5ByXsKH}`X)r`Wz+H8wLhEB}iPT#6W?ON(0pf1_3-R z$o*C>xoCCAbv$|)>i}OAL{W$pfH(M!Nr)DJGYH}+)LDHO%zkit=1Vca(JZ%U(Ef_US1WZBTA3}%-XrsUtLWBs|Jyj=DqBp5$@7uJI z&Pu?LvPc6$stX?LaiirorFa!mcqs|ycv(`|F>$AO9nu3jE z^a;L{kWvCGDc+bMa{?_X&X|x>LMQ&cGNZ>PZQ3Gmp$kKGp$iM-KLe^hpMo>LL$Ev! z>!YMQ&!g-tmxV$IE+sY>??MVXCR!hFLkc-2cH^}$Gqz*esFMi)^%G4G!>mpe}mEh#dW07r^4CS0F@ONuderYIm4yu5f*E~n&urDZ55o!IeNlxW6E zH6>-M6!8DEG_#u2^l!Vip0UzWdNNy1$baUOO;@UK$Fx}oJRjw4JoevvwvwB#E)^mB(QQQvAOZ0>=Ir z+x3Y3=eW=IZ;3h4r=SJ47D3xfbGD;3M+ zatLMqKjX#%i;@CjBAMJrJmm&B{W&q#uf-*fqI^R%8-Tg%o4P9anzt0E>A+nMrm#P` zLtHXEJmh*Bnc?2NjC7CX*B}Dno%^LLpzK+d6}pCtXS`e(HSE^-T_vQ)*4^*#8a50l zJq%=)BvLG%-HC^*K|NRTq|f&v^JZ}yhrHacya;i-9BP+>CpO<{9;O#xOUw1Ma+LQI zEcS@#ybu$YpKcR|5Pgc@(38V_3RDwkq%uqV(r@!6K0YQG9RI#DO=ewq6cXO^@N7T4 zVYkc68Ks$!KU5EO^->=zyxnIuA=#SdWs4-_@0Ha?N2-AZ=o3Nax^ESr7i+Go4j>Bs z>TYaV+_XgoJFMpOOL5V=O*4%md3__h=(=&e5zRLN$14O(!Ja4(0}|5y^CB;gLGQi2tDCE*{q3U(Yk0Y)IgOxt zdYZxiqQ}tLe%|urH^0dF+`e#d`Z0V2pd!Lu$ za}mQlfxE_=+mm|%zvmlWx>uQtvYtHs?Cp%V92_U?zDhcAmGw^8dC_t%S>?}_#` z{;NwpFHSGF4<#2D{(h_P_ZRv#A4&pGANM2kQCSyf99LW2ZJ#YyTe&&QHOA#D6=xSC zU*YL^p4a~Yh0-AM@%a4El11U~=X!f`+CpFj?%&4QfPfF5 zz`2QX&)?UrQoie!+lF=d?H>AItSFh?y28`t`Zc*6y|W8Q?a#dd296f>JOcjCvo8O| zt3SI*ul9$h)vs~QJ$4LML>y|9mztzonunV6{G7W$Ue;Hk zNl1>N6QHQpLl;tBhEPA;hqqFQbF^S?5=0>ZuIKgI5txukh9sj=Nf@Q06%NBN%iiIc zxB-RotI;9|nh4!0GF0LoZfJ!b=SWJ@(W3d7Wk&E!FKt7A@{xq214c2f=iEz?Fe&tW zAXQ>jNK!}^puFmR^sK_LK zIWael?d4Ob)P>;=Y&gOB(@F@M#L!ll&rua7bVfGOPQvthMz{Pz(lL*Q|m-xdMVkIxRhf4v1p;KNZd+K zBnD~Wk1n(&B~#SuL|8f1p=6i9xasdNmtqC$aKBvQD)!^rT1ETJvy%0171sP_$iO{XTfL>`7cFJ_@8i1)u)(0#tVFDtnzAP z;&%=d+NZT}hnI^YaZW!EWMg04ScNl31ktX9Z)oG08w91p1p}DGs$)}&Liw)E!$8k|G;|wN33OO7R zLauFbo4$_7UUBi59*9s%_E;eY)?wxNlLp2&1)XU5u4~FYFPt!L@4c5@7Lm|P2z}m) zKjBm0wEC`y59>mHD@?#Xz;?%I$;N4YE4N}&qAQl51EXsC23!=5SAtdPDnb{`vSQMr zV^Wz$aP>*B$EVI%5Z)eYu<{H9{j^>9K5ZM{pdNI3r^=lNF{xXH9#p9UyE-WSj}C>` z$sOT$t3%vWjUE)WWM>UfyH2KAFbs&KNEKM^3&nD9#^Wj+nX?~v-=p5vOGDr6i_mWp zqI3x8;ah98{?QO0M?*miHHHP-TvLS3yo6%~(!sl8G0yn}B&SoU#Vd;OlwG;C%GexW zH}IAAHWE^BQ=;98v8n@*c3}r4Yl_4>S=YK2f4mtd3S6Mo*^RndqDYV|3fnc08IL+r z8Yfzv)HR%3?+`K<^WwY>pURVt6Kq1ORqw8)ukzV~>XoKaU=^QEL7J7NQV2=s@=db~ zY2_}or}75cj-zZ`yf9ZJNpM6DNqN@0SAngC;@zeNj(Z3C+Uo7-`)4gvwzb%)Fevrc zZ|Lb_*NQhiK_zw*3B8C($9p9RFN*lCkI|X$ZNx{qKrKP&*qH0)NF(-98%IQvaV5Os z`g-m4oWmZ1Pv3B07xZo=f7bw*xIAHSb~_E-;S-^K6X$QjE*v#bGe+f_G>uR_Fov=k z4yAu#^fg=^wfwUdFr!ld+YGlg8Z!o;QpLfB@LFp;=p!Z6-7R06rM15R+j(YBD)3jX z>1}i?MP|1ML`HKWj*v`%xj*G?v`_7Z1h#`tpz$Q@29)9z>p4HANSgwq%jTj&w@m2yKUnW$*G5HJN1QA27?5#F%OoX@0LVhOu@+N~y@?%nV45 zu;fG?wJP}56SGPhN#iio!!5=oHA_T3_bR^LevNu9guESPoVAAibCAj%NSX^I6O$Rb z$Fi^oK;jySBxL<}Acwj7fc5#~A zBQz3L<`Sc779`>$3&=6_tBC=L6;S0~x?&*~r|qngEd)5HVH!iI1(k5qSCo;nqOj3DzJiv$yO&nBPk?%R(~DMA~KZ*?4so_ zA^OoO;>~?4*^+~p<{;m`=KdDX zStZ%FMlB4{3>}ls<cgGP!^txL!iZ9&>1=E6q|b+%s1HWHxOis{ueI)**gbs4cjF zt=Tv{1Ff=sF$#K+KQ>8oO;Q^N2(6em1slNtUjQ-5DqZ@`T<=$#*$5z(JPRQ#La{Qy zG;-3c-cTDwSf!xoS!oK*`DXpY#r1}Ly!8McP1Y*kZYwmOqxt*!xi&~@rxTV{y`Z$X z2R&Wuv=%+s7#&H_9Z!up!#C^P66m(H;MK?#&IJZqDHmHqoZj|w8UZk_4ZGR`r)4Mt zyVIn@1aR6awztX43vK3}rrd}`mZ~Cfe~q)BgheIz&Uu7Y6tVK=d?4#{C{Hg-%Yjyv zB-)Ath>I%Fk-4!!UV-rT<}JH)&kI@#qEg4^mVI}NLBx{K)y#5?JY4=RV5W%czcTSJDa;D`4Hnsv6w6KLJ>~s}gepyo0 zhg_I6S7?S%*icP*5}06bi$Db@hqa7Djj0wFLZ&ZkM~lC?;rJT6w{}?LE2Lr@%qokA zb^q19$9%1V8nwK>KbD=k+M!3L!Dt{^g(DlKR==zTX4(6?u_G#@R-wWk*GehT-suir zpDfogq++z5ZHv>QF4OnyfGeuf&RIH`WYg3EfHy=5Cl!Oj@yM_jut1_irGvC~Vs?0; zJabJRr)Ay0eF(Z>%oB~H?hz0I{8phv3L_Iyr?M_BV)r^QrCF|}& zW>J78g1F)WYX*n1rec*<&UVMDs9m5P9n(UmHkP;xux7z*s#<|16xtZRPLCg8QPxv7 zX7ZPNlx{jSk+{a5H&!QJK8_*D1zw4r03e#!g*8>=`B_*&F9rww2ybzpiKVPz-#EX( z`9LhRwd57Jr_RPXeg~*E1~6Q%6YkiV51&U^+Ux5mq~*2Hoq9(Z)IcD=4(#RDHE(iVAdq zXm$Gv>VD&LdOY>}wz|`8ecDa#vf%Sx2(_r<%CaJK&d1S?sMoosZ?KBRqx!$+?maHR z#ITG7r&IpS|@qKg+kmJjEEG)%Z?{}^(;cm8$1B; za4tzGto(9w$4PQpx+{J_U)>^ql4-4MFtN1i#w+7RT zMP8h`#1H#z(rb}@UUgYpP%fyxYZejf1Jp|NrRn(8ScRNm4+03kcuBGrXSl4cvDoth zTAgdmVNNutj(2TAP83&%p{!GA`|f2DkIqQQ(2K*m(!`tQ?=jQV;Dhf}SkdiwJMa@z zF^^DVzOX?q?03~Je-5rv`@pSrk`;4l3ijhWqx!_3Pj$VSfk_0u&?3nkqc*ke)O6gC zJ`=3((PLEeOol|oxQF7iU(6z<%-;=mLz%V(ol5Zru5T`wt>TLhAX~I!!(lc^s=#nWR00@dCUvG*~mH;QPI<5lljcp^}%-PD0qgk9FSAu{3 zqhe!$M6$=tyI`v9v=5Oyb``6wrkS{Xfx}X2%)1Rr!@G9LnU)Aq$jct*T-f7Y3H-UDW0pG=`x56d5WbE34il6<Y+ z11~zc=(~%ZPGy33CHv(vRNOKdycmwhRD7%vKt`aGh`$c~V0IP@bhOKWV0KTDO-NMU zJmg=DrxKDb20NKZ{5ID@CN#CqJmki1`XYE$r6$tMend2u_82z4^xX3ApgzEGm&Kew z^==@WUoBh$^&HnK7L`)>!y0W@X9ky?Xr!>pA*55jIsa1f9$PF+<n)<&B8CBH`PAdoVgl@Ph1lu4=2#gWcw9KX#nc;^KlCqmc2CA+##0j(I(c7gXxMF&O8H>5#CYIS> zA0)qGSP&G+V!;pmiUqN$MU=42#(OZ7fR~+;y+#P-6mKKu>D#gybNSFVyG~RkQ+B#M z6zVbkvnoaO6{@r5)#@{|XGS5J;?AVGM5+s*Y}r*o4sy9n#tpetBJ)LmET%G!UkIDY zOHfxb<12J%KugUyCPYDHdocUmc@ZGl1PWA0tC)`hwe(1}*pM0F5Me8Fr(6kRJiW7X z&js&Kvat92f+OAq(gRBbB7<%OQAM%?#nVdaiT*32Z^1i2!jO#|1%iZ}ch&?~5(Sfd53U=$xS z8E^li$V!Y zv-m6u%CIEkmsVl5YqeG@*kBteCVZ%{Z*Tku{!e}BOw)EjxoBu zm4+%Ih1BaXuoK5HLbmDO*YUivg(gZMe`st<@#-s*1YHRZxci&KX!?dZqveKhjMUFW z2{pi3ihs8&u}tUW3hW&PW=b1-eN=A!;PZw98w__KNw|3~C%d80)Fod1<>ZvAX|NYu*wn$yKV>?M;B z)%)zDnz78vsgUToYq zJf_0%kGC5xcq^r`%i=4TfBA>Gg`&jq%~uW7E3{o=>tM)V-b6KCGc$P36b5k&ie9CX z$Ea)hAJXr@q|C1a2wg`ZiP%oelK4w-Ba^W7)m9XE7PwmKH0ha_5)f?C6MvK?c{T0E z>Q*%b6RS};$S19yTbo^Oh#+-26k$#Bk!E#7^QCeLX|b5>P)`AXtOKc@GAtgfMP*oj zC%cG99HJc`F*(sa^58HOZ6>QOYaWD9dNY@)ZO=sR+Q=o_7*xMG=-{$gc`u&dM4l za%*DGLh@E)an73hiL3B+QE-tg^~-2>6d>EYF+0gu)R^w*rVQYGm01kfN7_+qSTZz@ z$COws*e8VJ{!4PwEF;fANCnUQK=EH$(w5zS@ye7{cB<^5__MHqxzbY11zS98B7K08 zF9Xc10TC_x@;+Xm6W*<0x^hNL@z!6~W3Q7QQV8 zAt#yp075>{X`*5lr z(W}SeK_#p$k6ELZ3zeb|FWh1X-lGb51x@^+`;4s!S=l0w`odBY9`$3%2tqI}0 zA4Z@eqpfRORIYFd&b#Myb7B7-`740?@5KK+&-4F0zh7isRwS2|@)-4n*%{-|ym{H& z_17sAnJ;)YnKN%%nh{uAgb@~D1ubv}$K#rpp4v_UHkT$ryYdDW(UOW~hC)ngOPWBU z?|s#nopE|II;H5aSOZL3`ao5a9;}tW8dDRO?n~RW+{$~K~tFv zsok2I&1?I=dX*gxXHwtU;UY8apo+!?^W#%5*chM=G^o)6N~yq*fU~vTAyTY5E227h zIC*2;ETxt?3DIVreyB_b#;(*NecC{MnJaN7=5IJ}3w-w2&ixJM?Dx6?dWUH~enuhY zm!mK4o00p=xZG<|R+n!CXarK1V@Bz^n_cA-Q^~T@o8N6a%zIEU7p1%k5Qx;cW=csl zJ8BMONspD&eg2`MYgJ(g;iJhV+9M$e0g4PP?uDWnQ<%IiofofyYV>QLojAJ=bt*mI zPH?F~nPz0TY7?xCY<%X<;HD)6?ATU~%@JugMWoTBtb^EFmjWl5KW)@`c0Z&ZdAYE=={#?9WP; z^vSJRgW-?*ITXzeS36y;;MI)>C-QxdWc{n}b;O8Bv@Dq2TGU`;UGhP9*xeM(>E(Jw z%w)l70;gzoOABR9&djJrIb+(*L~_vL4ehA7G7wz`-_s48tp;BY+qw>AtC^r6YFJW| z!c8Gs)M3L2$+R~h+bRBFxzKXflNzSE2ZjJzKM zxB#?;VsnTR^<3_+wG}5HuNwydutet1Jc1!^7uK@R z=19ScYrcc!zRq3nL_InVH5R$P#%g=5O)Q=B`T;nKkjF_9iP}N+W_>((H5h1o0tYDp zpa^^(h2e**yQMi+S?ANrG-En2t>f&wUYpFA1dAxcNFRKPED4+`R7AVN_#+A zIr}Y!=%RZ9@11I|J6pZ{77H0-ElSb0I9OQ=k}Y(I;Uv3ciP6_e1WDMfQ=Y4CFB9VI zkWqN@fwf&Qxb?)hI+f$>?tAIX^=^6#eVRUl?8dh~2H2-DLQG(?fft-jNU{N?Qu`NH zGI4@HKG%wsC{}};oT?l8f`k|`n?JG*j@;^CSBjAp3!f~k8OR)HdM@lT-lSZy4OzjXm z0xV&igkZHzpBZ+jf_37RoKD7K1LG56%|`^=zOBW>77aPf=dLqCU?Dc-JMBV2@PBBw zP2yqe{y~tT!+?uP-F2Jz9!50^PDPfMyJHcOti#C!F%JGaOsBWEnoO+^u3OSeok}Ly){BaDR3_++g4srmo(pN(HV-Ky~4{E49 zre4Cuzob&OWdq*J8ck*igPg$4NdjI@jy5iVWr-`tXD8&fp{}5YX8jpWW-Lyq?m%tg zcTQ37TBMjPP4h^-I1K6pMnt{Q*p3}#DsjUm{hl|Y(_wvIVimR^B2QZ(L`gwbyfyU$ z_=lUW{M4J=7M~)t+L69Kp;;vBW@?Z%4LXzOFLb6uQ6{n~N}E2Hx^#6aboJ48D*4(; zrCLAH$KdNZi{dk2B!d9DoC9xs`&Mofj_1L92(3kTa$bYIc+&tn<(DE;H8#q7Y6X-lx1ATN@@a@1hjq{{$7DSnCzF#W5 zMT6y9qaXEPs_C!?8!i4ENz-8$G+K-z*#CsEFEW``mVeF)=35Tyt7?*5-(0Q{4hEB}*uVakK%$`3H*K*6Nc8 zWuhgJxRKX#uBRLbM5Vn_w3^!%{?XYr#P0Gofc%`cA2*vQ9gV^Kp1;>ln8&d|FXCaV zn?#bA5Hb$wE&%Jl5y+A1xQ_eIZb3Ba-B{uhy!Cj@vbA}_vQ>T~IEfWiVCX3GL|9$` zikmgll~C}kgk*R|>)r5wcjs`?W}M^heF*Y6Q=qQ$)yeRrD~lFEqw8q%gao>ciY=}Y zNn3|I%QUz}*jhk#o*g$1RSX8aj( zqK>X~nj%IA*9oCDl8%$sAvICvAN;T(2kH{z~*(6{UCug^SW5Tis$CTaoChn;N@B6VV zt%)Vouy^!A^weIRHj;q9%lBxBTh$_xJTMOtRgJ9NP3Gu2_0F1>@O%0^wI#pu_b)t` z&ch+zuOe&;PG6FfhvQV%;yLi(XbK(^)mj6Yf5p1^+Bdp{;4out?l!Zo;AX)kp--qb zg77C?X|QGTBx81R>;Hvv zFNd61N=s$oh#h`iLIUOKIJ2d4WG#jKVgJmp4^v%*KiKHiE^7~L^x`z=$j7&`Q|3GH zHs={uI95JoWKV_15^Y78sW5ljUJn{6$E^o9Qstw@oO*cPiY{vMq-ANK>@SuOir$yKwl{w1lwK0=BMYy}G(z zLf6((C|lKNAg6oNgu2m14Sam>4jj89ia3Gi04{{#H1!tco7CuQ1^A_fpI2@y=?c`p zg2cG;F`_D$|Hm)is36X^#^S`|f}n6+e}ViZs{z#P!%4Shd730)EWI zM6Gv4B1XPrcz-Ypy-W z#m;~@HVDNGis7Z(j`>5d5f?8)VU{M?SXGI` zG~1OYB-wE+K(RMrRUw+h^S&zWc@adIIjYKZT&{6Lo2>re#;n5 zyN@?=1RYH!Ji`~lA`${zmGE>mU3EW{0zp?)*G(n7+8>LGN|)UAnD*hc{;4?Cm5P}x zzqoE*aWcK7!jL&S#x5zb!hObQ0Tn3hhpt>sGeXW0R7p9Ew%UY6Hn27pxy?JtdjW6L zh=pn03ol5UQHa$#m_|@e8LfSZCul#rE}SQ;U8(z~E^B$lhlqSDzet@wQn!6X&x#c-DRkin1jl2C!!cBn6 z3U`DTl}0xm9$i-g-xVRL)n;Pefha#`HO5WiqXEI3ST>yZiy|x`7&_HE1f=V&5EY@j zEd){Zct6O=`Km$rV}xfbr}oL8NqE@j9`MngmTD}(4=qTTXIX-7oae8{PMYII3#TiR zynjvZO#pck6doZ&`w7fuF}f^E4rLk?jc7@kTzug3*Wk7noR4 zc;0l7D6a4Up#wy56C{!(Y2zz;Z}%q4OJP=yws$&M38tHk)GgS{$tR4x)rcg)D@x#0 z5nR9He7AyJbk+Io$C;{kkbPa7#3Nj-rYzqvyN@_*{fOe-)K;FgZRl=t5M&3UJTRza zV7p0N?7Dtfg96@OgcoG)B7=iC^9P;x0o9h-tds@1N*e*2L%KzGLopb1bFUlZ?18Mh zE9)J7`ioxLEo!>=6Hi`@#z1~4k@fBuY~gEhbOP-M6x@1Sw&T0VY_;Q^V@j}$%WN*k zJ$-B`V9#5ewW~59v}eB#$X%B9orv~4nVQ$^4a>JGYuDfJZz#(kqBkmbh`=7HBt94I zWVV~i;?)t#)|e+KooQML;n|OzSE?heSBQ_mt%mI^WP80wTwP&nkRZz0)fz-P55Ekj z|FDidF-J^%Ouy8MH0skyu#YL)!B-(-ps_X>4Y2-9Es7bv!%->2cJN*O2aNjiFrP0n685T&Y*I5Zmh*IZY|@YpzW@MBt0 z`_qeldh~Cu`lwFlB@ZX@{fGpY)2(U#;u&=jmP2l)MMPg*`sb8lT>v|pjN$VS))6jJ?Vr0)ii% zm(9SN2nP*mk#&i*(9XmccfQPe^g}3G(NLfDN=l9_R(>5cHSV&^A-%pe2VOTqBOPSh z2I?d~CNY37r6L&KNREIQF|N240rMebsR-(~$_Vc|9Mh|7QvVlc-xQ=-w`H4|m9}l$ zm9}kG+P0mQwr!)*wr$(C?VI(Vj@#XFqtA`#hyC!q?1(w$TJxK0?J>se)%q^v`I!tP zNgNnT{838G8)|by-|<;>>z*0RH)Utpy0Cz?y}%K(#W#tDCDyY6h3f|Bi+|!2a zw5jB0wMxPhtMT)DNBhjV6tq!PYdX2U|Ha+zbwe_;iJj~PCd<}_MYMe-n%WgI0*=wo z=X-_Ks(rWfYW?$&m7E*Y;xC-{{|@I2}9B{r;0PNv*r`nrC-;%XKLn4f!(O35Ia2ubgNWNFMK$BAHVG z@^bfGGyc8}KdX|9dmhLJEVWKAFb8I?g}bZ|7|4=OO_WZYI*HkTD_F9VeRtg37C(oi zLX)9aDoK1e+0O%BP(F=7isu1Y~lcKrN9z32?7cAr?5J$mWPn4!)6hm?z z_eYHpxKnOX69AAUpp>dl8*>Q?PNGI)e$h<-OM;h!Fq$#z?hlM}e$a606slJxW4>T5 zs9^dQ1Sj8=PE&Xa?tQsNPgYlh4)>4pm`TBL^OL`aMWEHe6Sk^ZnBwRq8%mCmo!+^# z#XxIgUV8Z9Oc;0PJ}}Vx^PqJnc5ct_Y6{1P-d@@Y3oQpliiana^B8(tZm z>um3|F!!g&z00|Sg#o&3Ix_C^GSv6IBUU&c*M~E+#BtfCwTpo-&xeE8_irBd)!fP? zYxH`0G0DQ$`|Gfa%^7cV`^U}L#mmCMM49Z`T;1JRwhgY)%iY7viu?kJ=k46%BmYcx zW#!bDlarf=d(FD`qC!{L<3iQi+1ef4E1mgAXU1ncH&4gsOFPfA{##8F3+%~5*;(7x zz~xkvs*bbIXsPqXR+aVH)!~_sY~qHt>=({yMt90gE1P=6!R&w{(#D!zH&c`QD*4F34sm+l|Q}!v?)7Q&&c*Pet%FW)x<%x-E zm)7UQ+~e)g>g+31Cr=d&cc%N@+mcRYrA?LWl#Rs#oeb~$-WJ`|=i|w{^V>43$x<<1 zR}=tVBO2OR*;MvfS5=3Xn^VhD%~e;1>?mf}0VU zo<)_;=BKT<*WH7BG~So1v-9)x4A zeK-hNfs38Ne1#35nmQSn#>>W*?ikI(25^=#nmQlG%eHLmaXr(gv(dB&7)SZGonve5 zboIP!45GuhHcwS~JtCD(;_y5(9(^%to3!TV(3CLba<%Su)c^~d)WpCccpQ2?M-r+rBws{-6W1pq~UQ%jV(yrlBh5e%dquVi+k=wK?>%Ms);4Vevx+P8CtC$hRUZc`^~9;+Fxab0aKGDq(6zQ|My$F0f~gCTozBFIhMaQ6@? zCRT;Dqqgyq2pvoim57X7+eD=+sd005g)Pxo7-{=MdORndcJ^cPwr`GOAj{*CJs?wF zJ^k)M3w@<)F(48*%D8~4J#(%Z0I{Z3mtSWVO{zi{`CJFmxUoPO<%#zse5L-fD4EF0 zuYL4tWv?AIvl?>U-0|P+$CxUA1g4sMMvXIV3(KG+T6>sYN_dG)cQ|gTuW3^dRf{;Lczqs-S@S48N5zQT22Ixi+zff9zhWxHaJ3 z1H%*!O4IYN>n^8BrnO{+jJk5F?jRqn9h z`;u=1v}~D|sFv7_YxARsES&CTM|@K)TtI*nP`G5Cg;ahz{S!y2CSJ52yKh13687E= z_(2&YawxH9uyFahrY!%8XDrE7Ssihl?(SiLh8fH`;g!zmvyrvM>$-tidCR!BsLct1 zG^e-lRz44t;y~lH(To+N96N+PrDH||v#g8z+DMGCOMJzeA5cN!_P9d~Vkx#5++r0c za}CSo#w67MO zmFJQ{EkDd7Ge7%Kd8uo(d<`7-W%mqd>QgWU`o$d}+2SGFhn<{_6sHy}Sbl6K(3IUf zfL!G>%@IG==d>=V^@%izM_BVbgop3glBFKH;FPfvCHJx-?vO)Wg}4^*u7alcglXDb zW{bNs@G&DgyDvS$E4_bj!c3G4CyoTSL4P=wd^gVcQbQQ{a+BRH8J%=&=_r(ye*(g6 z?U%uNmq9vMz{4_FLo?f;zWxM)xm%7H`5kaIKnH%^9h{PhH7j z=pmSDd{h{S+Wy@qiCx+}52 zi0n8sBI40gcOWa+xC$2+8ak&f26XLTNLpF*GiMAn&6$v49RT!&V!-8me}=*zTOvpH zi6iD)3A&6Vi0#|4i@-wIL(OK_M zMfek7CoLA!?H;D8v>b_*IW~E1M9zIsFVOZYJg2ADNpRujnRw2;73^COkkoK%+Lf8A z7TOB=42Yvc=OGX3y^){flg#=SeO*w)lK!9unjtgEZNjT0GhKk6HANHv`45It=!r0wnK0JMtw}$Q-agEi&SA;WDg#BLXQ4K z#@~hqL$sPZrC?JX)~$gM$>IuCjQyr(y9h~!gc&cY7pHQ$UMlNk96ZPq>Y<0pmj^{_ z%bc>;n!rM}JkCD~&m+gVMJ2v=M7Wi2qEa)%SXeuE(KLC(tvhhRby0Uv&+LI87DbP3 zV%%3{`&c4c#w4YUA2Xx?sTCTZ+5swfostf+&h-URY>L5)I14vZ6;U)PPamx4&s5mu zstT@k!<#u#u&f5{o5mk#T$Ls1|duVc2{~U{fdyWl*b2l9i#-l`7lw zXU%-e@X0T`iIC!ck@obu8-8K~H#K+a4E3sOI@SHkZ8r0_##(Y9aTY1JSHui)ZJ?7o zOPRkTwQUq%@C!uvpJ}Xe1*&^=u&2{96htcOQ9H=uwQo*`oYdqiM1P2|N1w_U2VT9V zV#g`JlSKK5kskhGl}bVemkJ&EM5?J4IU1N~=2h94dSXOk$0r$!K_pe&_Zrm1YE&L)Due*CF3ZQn}5tjXF0 zY3YD&hn@82mt9D7@C8Do@wNRfC{ht&Xn%p{S-XF>Z42ZBR;4C3yLYZ6 zYtXx6!yK=9fUBsJxc#W2uz%f1Hl?wqWq~YPeoP;mLTh%3waTJcb($?(AEAsidm!7P zDx!LJKubzyz#w;yObef=YDgmnbnd$XplV;bK zZBK$ww5R(PXd`Ce!}=S9Uvx(IEg3nBkM7m=O?CCa2r+8z+mM4E z+@Q6J^JdRN=!ZO&R6A<0w-YiM&ac8Xh_8L>TT9a>(>30PxQ*XKJpeW5dVzKAK>Ya% z*`#klVW@J<_rEqo2d7%rJL|n|!GC#Z8o>4eK=5_vBzV3}ebA1|d|2^{<`9MM$Ve7I z94+dYKCk#+b^80(eB}bx%%BH9OUih0mc^abut2S9X{(O21-y+()hWDesqft5QeZb6 zJ|w@#&Oicf$vn>Y>F`g!e0x+icSh&SlQKx}u&m2GXw!p}A{6vbTdIaGLn#~&`O7e$ z@^i`z+WB$mTOL8cAi2_H4ffR$OO?iH)MF80v%k&IA5XK)Sv~?Ya7Waq-hnDzG0dK? zCxo%Y>PJG7$%68=%L++dG7GHW!|2{ZgN;%IB~Ee+3Ev<0>_tdNs!^N%z`1}mvaD}b zJJaMgpjUlORqLA7imGse`-1x0*;5Jgf2BYH002|}JFW45SP~i9>$&`G>i@%#=s$k` zr!!Av^u%9FqJS$N;pMiKMgI`y3*i9mLLNS`CbQ}ng3vr0o9=KBI9wO!L#(2cEN4Bq zun||XoI@0P`_bq@Gr*;tc5c2U>oh1DMl{JZX|ZxcO+@-5?5Gg&2}7e@o*1g^5m71! z_WVKo1!?nME#(UaZ>AhA26bM@kPo#VSlYk81;S>GVKUkqG%fG>K#?>XAmsT))V>b& z%$P4T7PuLp)N8gsE4_aFjT|MBTA}p&#Detgob&&6&H2BP|JRr^71yjYzHLjVRMos~ zpmR5jXKLl=mGdOz>);GCq^2x3fJg$if+%k0PN5+InL~TKsexONWFae)wk_K@8gpJVZ+i(pSQ<-Ouxe>BOrLNkIqc4glEAjBd9&jn zl&{UT^Ppz6lg5p0I%TDxh}{-`$8?lWGhb>X?nOP7x)Uz@It^prTZ-2_d&l z{HBBZ%28)-Kk_l(3+lO_Z}~I}-2;VSHjvU8Hf-lpEcG(I`^R-`XhNp!ulW|$zg_5< z|2n+?%WCYun)CkSxz!OrVcz#QUN^{m50W)AL5d^?ir<=R$nGe&tB5t;F)UVhkIGn> zobXa`q+~Tj7vsSye-?F$rg-p$E~v>sCd5HiT1D2P;+r-OFJsHkgdc!9x*drrFAS30 zE%${8UdaaYlyrp(CjF&3E1*PlDFuhxV?1U6sqQ=ndYAp@(Xv(=cZe<1{S4dkgTGwc z8W#|=GP%&JT}k}(P6G|e1;B+q@igj^~aUu`_(^R2P3Zsvhl8SXz9D0WZg-kXTS1|%+`He7p-Y3lD zYr?nGeEumH!4}o4pWhhu|Jyb5A3m~HMuuj3bY@n1CPoajw$>&Ut{`$5()9=UzxAg1 z{rl+IfL*{@hJlxKdjzg?2-K##>49DR+S)X(bsZ>n?gPuOWDz-#~l_EjDaE%jMejtQdm_~9`gE1*EvyX` zEk#1dMiR9eC%PSj=4kB^UIN&HX5g3ANiu9H|K3txs({#1ezCKT|XFRY!!d^rZmtyP9gR*^q zE)WvpfZ#)HWZ2I&QIeJqqNTGaq~{HaKC#qhqelG3&4=+!dNtia05{#hH&(#i%& z4g#nmK+GV2(C|0_C)>~0!^?jHJr?@nC=W;JHi1xv83#DwJ*B82Ai9P9_ibDJi+?)9 ztg=zZKs@5dOsT7^_6I%sd6UTx(g%i`;Wc@2($NglJPmkHX)(>T5SX%f@=26~BIKyZ zY5gc9R1rbNnUJA}_jK}QeTBdrgc;GX8{NmK?AL^`&376_;)rm|+k3?WTZEhPA~%i1 zl*LOX1t2$OGg9l5eyb7+zh&dPH;Cbr-h^PH%m4|qkGeI*mC9ZNq}U}?4M4rz2)5B1PCh-4jVS?l!h5q;hx&OZU&%+4 zItNz#m@SjJ{1skarts%w#8}cM=iyd4RNyD+PZoL|3S;~$H%DD3C4x6AA6xVNv$irt z(~`7PzKREOWlx(}OfN{uCu5k<3FB=LMUcCySfSfIB$b+-_Z7!Cm1&-MahGe%U2(IJ z&iE8XkXCjk6MGaU9Mvk8koK}dimbWmJZ8pW$=wydXToqh<*ap=9QYo0R;4{zM{hmX zFbm!tz!bxn;%GN zVfH6L+|I53^wl?|lO~awd!sUfQe98>^~;HLa69k}D-ly57N@SWRzE2}F7`YdIU&Hy z1JN-tYt4f6MKaJ&PEkEuOWX6E&0{yD8Ofv-TCpZ>+n6&h!7j~mAfREByDYIEA@cM= zudA4KtDX=f`liRl#Ki9F+jBo!(cENuZA8RiQkg3(QX9>%`h6H{lwj%IoSfe=;$0D7 ziuCmPwZ3H!duw=vMlX9>FG#(+Ve<3b4}3uF3|B87yjIzWOq6w;Y!) zeYUOHWVf$RYs40_?peSvBPeCRYvnyZ;&2j1dEq6)5gr(rg8lp~RXIO^lghzk?%o=Jx@0`fZ?73=Haf?59#MHy ztC&FAivkZ@tw3XKqK*Cpj4~5T03R+3*>TxuyF9ga>3Rk19KkE7rB11!*8*xQ?PWbn zm4@E41t3f>;@9Z_d3*li>$C7$`1~A@*hse#a3ICofOUU@*?AVHlhuJ5~WR2pG{K_a_n}P)@81c*?WI~2k&;QoJHqkyNV=k1@P+l;-aL6 zcwl6lt;E_$5IXTYOgR1!%l?+x5d;2t0v%Y|&FU~e%Ysv=PW))6GSSn!nJCkZn8Uj# zeg=0j1ENPzX*Ob2PQX&XaPeSAV#t#?L_iuobJE96zwj|Ik*Ck5xdg1;HGefb#$9@< zvcXi--@K6COZih$GS6eg$QY4ORSkz!#N)Mqk1$Gr$H&RsMZ&Vk&@Fs*{E@oHo@y$|do1Nfg>K)FtXYwiBc2Q; z!@ns1MqJ6`yWU8eS|r@LCwUNQl07~mH@nQ`>@A{t!7$IQK7Sh|KH6t`-3~|(>Bu6N zaN)@`my*_pJpTnP7LE^xL=LwVId-z*SEy5l1;Qgtb%|u<1E@DQqW$T;-R?^Dl#d0f@#X*qKoCgaFixEuZXz=-< z;BlZ0R_a>T$oNh^oUlSRa;1vJ4s4Z|G=!o_K*xk*Xws!-!HcFfA1cDDY;_GI_GL#r zR#7|KnfWi_tbtG0_2pEAdNw8L$~}UuVC1mr`^gJWu~2#qr==uKY(TNn07l2;-o0EY zDX9(!+!C=s%*vZn!U{|4k5Zkb_7@S{Z?p9d5NBx2Tyr0vZsjkP$)(%rkQw>M}zNiIN)TqT(J+CAl})sg&xjm9oTRA(y@Nz)=o_814Pa`YZ@WfDOsa=lJ+)dKG%5&24I!Ki}xa(@m zi@#+!f_)S(NF0iYsi_|?iz~dx2mn^tm@tS;AU{#XM6u6mbx$bt{=}}dH4#OVQoc9{yw)@ah?(cj&U=d!mr48s%FJT|yH>;%5448V1WKBBJzowYXa| zmPLZ3&P4(Mm-As^XNgCoFCwB1uZJ*=k3xG~7m;F$|4Z$*#mY*FV5S}sYfQhyX*p!` z#XbFDa8jlf95dt(g8)ze&+D(=Rt!C|h<T)WdG>T)Jv6xH&20=g&6{3+DWRUosEx z?cf=vk$i=i?vQXYt&xVtLy^DKRm5SE&xbXYTKU?c_kgZ(kl%P-oQV14ObhPm z%Lv1<`#Ts9NQOE68a(XG_||0$oDIe^jF)@piRQo2hQaP0sM(dm$9t}2^ERUPK^CWp zfvW^zwP(3vi33*}ktINXhExE>2pD_sb^}5J73jJkb-G#zEeZK)9=F`4FqVUk6$~m6 z>FIdKb|#rTlPhKeOdb2Gfcy`ygu#LAWPrKCi08%R|;LullNg`quGm$h?% zvJsHUq|m)#EyRW7Vw4q5F^kKnEFbM_KZfLXotBBz?xPV}ZwfLIjdm_PFBxFpzgv<8 z=#H}j+iA1smNn*e!?@3ZYXhO}f2|1cZA)WYK#>uI)s%uQI0i8fy7s{xab_2r;@c3d z2n5031NRh{W#Q4Ui*QA1O@xV%A%<4^SdEKQWM25F+dbuMAFv)JF!Qut-w4-ClFF#n z1ccP~*}tvj?X&t-{5VmyZpm~Mk3_>=x!-bHoVpBiP z$S$M3RF0t_KxCRXG$dK77m{5xjKoGH)+Y5a12YX=oAcJ)~pDFp9Kvl*Py`t|RmY z`)GhkfwqTz=u1V=PXqqVE(I*CoIQzZi)fi#Xii%heY|Xb+7=HL+0~IZQWuORqAeOh|{#C4jSzB9+iW*_)cQl53NQRgt z3Stk<4tpaK(E@yeubVIWXr)zq2}Z3H$35Gu1EsvTU+>afCc5c(bIRgGRhl~_#v9^i zpOm(#5bX(>|K#JI5)VlaMP0@rWobibui&ppcaeC!oW?Ejq#qBy}>8&xoaet|i z(S2oi8EYpcg3d^v$$Ug&fV&;p$-{pApnG3N$Ak95Fva^qzGoVzP6Swwcr}G*po?@a zPTVi56%7|Ysa4-aF4Zx8ZQ(0tY+M>{hQs29BYij|=l3RUyE{NNK>nv-e}6?JIJqV~ z5_W)&J4)3w+pwAaKD;N%6Nmdvr89Emc|cQjEYmm^f-QBk0f8iX%+t+r4oVpD{Fvaz zbU?i2N+Wvpv26m3kOHHZoPeTgQ!$Qd(bQWB-@m@GG3g#UXQ0Ry<4_$%NogBG3!xF9 z`Sw~a5aPO&E~2{G8si!iUOqKXyjMjWw=d_brybF{=2*kI5W`)V+)$)TbG>0j8a7m| z4Hk+S2BfIw?(VM24#r+LYDGGSB^>|gQWzydh`WIO*?jA~MoP}`!5={G)!DlpbiNKi zb~%=Z9xUIvxvOQE%%Wwx)YGR~&NSK_h*X1?Wmd~~e`u^FrU4%p2jgj?oqaQ;LWEx`mZ(`l+=8J}e?p03}SV5S_!b18h)N6N)sYdJa(m5Wn9UD@c=pPGu^=j2&yI zT-eHFCU{v{OlYrFJR~?JAmas>b(1>CWmmVaQ>oHQ=lNaaXCFc8I93_sQ!m=B&Tx9- z0exNBADl%~0airjti253YsXa?4H=-GCn7$AMtY!X)nsFw)|jQ3O?y{hfjkcS+(mnkMLZSS)pYAGPuv1O3jj6cnX z{IgNRAKSPoO!2oTI)BYz0t~Hf)%?3-KpjUgt?i+@aK4IE8a)m$++LSVUUh)xOSLH` z!!CRRsC)Ap%|_sD(9sQKr@+|68lACCr+2m;vGZ1?$UiH zLa4Ncujn8pL2K#v6QM0ziw|Iz{V$0lhdF{65NW%RM+xgB@=JDLhy}|! zKIf40b1QW9aTZP90n;F-P*Ge1I5buurfEKkJS)B# z5Z`G~ID(QiS?qD}daSzx-TBJ~N@|dAsg4`eOrnHW1kFMwh3huXfoX|IG0VV0xT*#X zDqvKQe9ooT4-v5ldFF<=JCL)NZ^Aa$q;0 zu=orx2F+9H0m={`^IZ5rr>nDZB7g9`AdiBeJo$&x*0eLUt%4*K4L431fJec5nkTZR zF9a;v``AOabhiAT^~`K+<8B;Gkzo5e0;@!aFkRH0ctl<(W+rbA2_%*rWC(xYEV?(CjxOpdR?WOo*HD26-}%n>`LMi`3xr9 zd>FWV&)dg0OBz_5f>2n&^6g?4azR>R|6H8>-7w;vP~z`Heuk@(7fX@?wu9ZJex?}q zqtuPRi_p4Hk-sVUYzPl?yZkajhd9AXw_%fJ`7v4y^sIhD22izQ#UFiM2MxcN%8B~Y zu|FXy2Xg&yM&{>VxprD`?Up+=fi(#Hx-JXo*(2W~xne^4gAQy2=hk?~k=YQ;EBSlQ;b-Nk3N8-e8KWhmM2wFftsf#6X-M2u|Uz${%11DaX zm%%jdql%z2#|(cBT)*n3lMEN^wHqAuur`jm>KVkz6iEb5uHDH@Q;PfE+kMwFz99X$Q~TlOKC^2+ zCSw4Qy0CsO1<#Y*UND<_RR=vCDa;OXGIoy3DS2zf_Np<;RA(^YjN zHyXBv?qo7$qF``|dT)VG^?RZJQg&j~&&3H4vIG~b2Y?tQb^+{GSaH^mG!)caRR`um=!_i*BPQa-lc!Spf*_38YCA71b zDv+cueG~Vh=oh2?*mD|jcS$HPS~I40-gKzyswFUZm;7`HMQWpo3YK~D2a}wEH9j)a z9^*P#UKx&dMi&htGCDe$29HElv;Jm^Ju^HT2E;DR$X4uG>&7X3zlyxX)k~>el~5v z#@%)#a9xm5c>LCK`iO4zkp0O|*-vvgv32p|&mZcW)Ztj~_V4rXANr5DK|q?3Z4Ykv zkRbV6ceT?i#E}v0CjRFqM}%IX;1x&-ak@=~veZ(=PAlaOjX!(aP(k~iDJx&3nL}0qLu{Mdv-454PWH% zPfH~NYBFOL3x<1misjL?uQB95#q{)YQHCB~5i09n>p9c4rYOtP|RrAwQ3#X+I%VHo^s>B~tSQ2I`<(C6Ucg8LW_U z=<6BUki@0>vEI}LhM|&PpOvY`rDnc>WPza5wZZ8+h8bzL3;f*iBqid}i!%`_vtU$h>$R$0*$C!kGp&8)4eN?EI2S&$0ERvo;hW zt7-QNE+MjaPnUe*v)po*BAcWrYxF+8CH{ejjVUS_b9E!v0RZXpnv8K0kzN1BHWkec zzqx}uIEyASMMDgONwoi`g{JN(*~06H13AOYG}3yVvAHwf-4+L6(1F?2hyIXgj9fOc z^$OFzeq#*SK&*%{$&y!{ljI;Toe@a~19`zGZj~|K@XcGs8f`Ugcf~9Gb8AaW8BWS* zhUZN;3Sy^{ylP0+=qmZgfPFexe3m91r7Mv(u%m+pe+O%?`WN>SiY#Gdk{IeJ&I6-@ z;wt^%^(Dtx4bdfF8~WZt$^wEj1(05?NUo3eb^m#C!&NzZ^SEjf<^n~ZcM(HSS%U5{XhuBeSJq>F5P1EZk>J|!OVfs04!f@$EXyd#bI58RJ2L}zl- zvNqiQixj-tB}3HP{v8)N+t1G{Y8^kQ?*riP-Jl_cv3x>D&W)B58e*YW>&$8ZJbwz$ zV4mn`4QookleC2;HYw9P6s)?Uxo(5v4ZSxcC%uW|3dL6(BV!rSIcGBUkNY|&6D79P zJw`&5keMu}jg|!ps;2YA#Hw#WR1IK(R4pL|2t2mt(Zxw_ZZ1gV?}LE634Zd|6!!RF z6Tz7-6G;=+n4wNvb_yDrfo@V?73&YVot9&X73n}H-p@I()9;Irn2>fY@{-n{Tf$@O z+brB@XMuvMF_t?nO&#a+fb0xvd+H1IO~LphTMMVzSKbE_iMmx*?I*Ht=Ja!N>LqZo zGjr*=51Mhx_gBr_=3S{zUEG#g_Jp9DU}<(^c?!V%MJ#)2T+_EdbR~%nR>~&RxhGqT zCY035WBRBMFWz8*7TK!uw4PCcDGB&1DmUYl&mWNnw9lUisvCUUN#ViBiUD|iy4cDm zuO7WH;`Ux&-G42Lt_wxT%#*|rJwV=c&7kG&w7U{-L#v%Gq-Ta|Dv8a7xA4c8l;)ij zSi_zUsjhe06{j@!9PIDQL8HOLm(W<4;zvS5Cdky*_8WuLJwHu{C-A>`(cx5?=nHai z5Eq;EZ(UF~wALlwaW%uQspqq>MTYpRM;8UIK786sN=S?@T>eHw&V!wq$`^|;wz1JE zczqGQuqeDjo`oxBK-(gB+W36)0ql?GUN(c7EOSJ_peKsq5PFzI#6-iu5a*-ve!6Z7 z*WUf;8OTtFZfj4lR`RQbJoVl&J_St+nzEaF*C8={$u1xHbnSu1U310j&;M2PsiXN$ zL5N-Q>9eT!RW+t3$L9Ec9ek#gY)f2`GN(D{Tl@&h^*xkusT_z32Ah!*M3txsv$Gn} zI~jA_xa6r=fv{e?LN;uV(53A|)}31N;(7REvGJTHIG~`1dV2^MC(6ml7p=;``hV+t8(l>9IiBky3*ux1iCFY(Z>>-tO=F0iVFY!R-| zhszUDkok9}GPi_jLHYYmt#ypCQ%7i)e!?m)z5LDl38m9(tn07i#$jmp(}qHAf2T^| zCl&nlL&#A;2jfk#s(UytXfE8{+dT1uMsC8Ba+_$KT66N-dI<=c3>fQCG$+9^AAV-% zS?+6a6mg3PMBfvNGBMVA_{lQ35vBnY6TF$?uZA$-n+mrLlvt?XnP6TI6Slt?;^1bu;nfOePnYa2);CHN*Q(gj zeQ6~!(+m#b+z2zQgDfok)KH9EsNK*U@mha_{z#{LE&O={8J&S@Hjd!Y$VO))rOmn@ z<0XWT^r6ySzvL8Y3}QAZ-9O9l#`(~F6xZzinC$r$T&l<>Dx|eDPY{0kZI@uG8A(;V zk!``=lK?~(&xXZvv20day?rB|Ec9*4=^lk`5Cw`1)Ck$683lQkVGLFa%W6A$aSsV2 zQGPnCQ{CS)YywjiZX+E9Swqh27_aksK*hI1Pyzaz7KW=O;m4O;=KXoh)g;^5o9ea? zq*pZvLo+)p(uJBN!RnV&fv`DLJ*-4~_Vf)Yu3PyfifbEqJp%tcAx&t{L4VmfOH(`B zDP_8k<|Qk8pM?+yN($;1?bJN$?3T1$>NI+182%BXXR-%T069MU7IJ|of+@i6k1L;Z zwmORx!YBiiR?1C$RVmipI>Ia|UZOViE%)B=pWqWv<*hj{G*%&|Kfzf^coKE3h+|W# zp;*5e4IiH7vGUfyMQKzvYM&v5GDF5-c2?GVF33hLx1K-g7eLJ03|zdAHd&F|`FjiU zFX`M7eW9IxzbCuP&f_1qOG$2`qEy3G0a!KYSB9um*;ZbIN`hu)ssspAhOZGrL56Wf zd~u-ZLPBne(_IMXYt`Qh1IgiO^5dwGnBcx|eUrw(z}>6|lNJa709d~na{n`l?3-T) z@GoM`zsNN6+_K+f8qn>IF2dW*X7OD7gaRmm(;&H6wMefQIsI%S9fEg{a!C2GuoHNv zdatKV1322neB3xz+fG5nT57Fo;0U3G)FQ*Fwcr;T4b{3`2vGjq{-(?a@neQgLz>u* zz4bBF`fG|`^g9J2#UO!&7R@NNXfW5FHEMQGO#aHEsUsqZ<5o(8aTmdX4z)FvoDfgs zE|}1n=r?7A#868BiI&*UoUNZIHL`-_-`e<>m!z(%k)JVge!$7cy4G69&hx){ZIMTZ zsxO{U7^dw~xmOcvho1WzQZ-+D^e8XsHV<2a5W9vX?U}JCecNU_T zg1qcE=g|8}mD`%r0Ov8qgVG;rou)4$4Z4%94ZiPGIs*uc(^jSaS zNDnDVG=xwT&7V;vj6qchXj&>AN2ZJ11yT00Eq^$o10($cs?8wCugI;+AYxZPTvz|| z3oVm>k^))IPe`zjU#zPEyp^J=DQh#NH4-*mM3fY2)=i^xn7Sv zgAvnx!0gpHKFfQrn(v=Z!i4HC54gb<(iXpRSXV}#g@aevViMuCZ z5-rCnu=lvCK~G~o_XlSM(t%zRD$u!DUGx*cQk%Y|^a6-6=JT{Ao~=a7st20(9HAU! zzX|Nj4)BbVT6X!1DC)!3v{*AJo-Jhe`p~j}GpDRvk8D+ATtrJgxlt?q_2UN{uBd%z zq`WSmXfxc`y2&HR-@L_^?3uWI=k4s@P9FSbD*b2P^lWYarY$p;+vYoMK~JAy!>1DG z%>+RE=*oXM2pR;#T64p_098Y+20(kdssK?^Jm$Pk+^z(vqXn~ZY@r8N(vaDSFn%4> z*Q}*^{=ACE7!d$2%XLAuuB?a)!GVRpRcb7qZL#xB z>9v#8TEYmlk+}e?kF_793>ykj$QVD_JTz{$bU#_At9rzLIt>7fC(-@JtrI3siAaDgcCl~%2StzCd~ zg$dzED5%mNx6tZEui6G`KUQSC+wU-Fq|(L1fZjN}LqfH`HU7xY51zM(5fK;t(-R8S zKvm=01}_O+rll#3NKC^Jy~xkz6{Fw`ABi`U=bFEo$H{+Ql|wxQ4%=n8ix~OrNvR-;7)!S?|FBNWeWUiWGpfT4&Ky z08`Eg8riH8=qjq87|H*3de>+<+G85~N2&M>x84ErjbQBm7e#}Az*Dobq%(6gvigU_ zTo~6Q(?<^-bj3B0;E~26o!(&1nJb==63`#hY+P9E%cX+kA)iL!Uc8GW<3 zK?tH$xuwi>gh^UhIYEg(RpC6DYYD`6c%U3r%4%Xd?l@XZZ3#NZsI)E*;9PjCPo^Mo z+&hBU}YOO`!`UYP3-72!+-xD$ zJ`TMy68?<=c~%&>1WBG9c|bCo=xc3A(1VAu{W z8#{lMv9^JjL;q3A19-jU9u$L#Xp>UGp^Pa9fjAyNMKjdN?IB5|JVm#N8$LeC?jYRfMC(}dD8)JV?u&D7a*^SS#_8QGGiXuxIGZxRZEqP|_ zZ(Detfk{6ikggCf@EA;-pea(g=6bw1wY0Dee@cr^+SNy|2ng)s$(>2Fm3d|y$8r5h z;qK>o#_l^|QD-nqSZ-3A9@Sm6|DA|R5|o?Uu=A@fZG+k&N3M~iPg5ZbW0>38siK-0 zB{ActfCNY4NYy~>7GB#s&|kdgd1#O*zwt)=w|M_88UG#ce<;jWlQn-S%!4-@#rPY! zm}mwP>MiIIS7j~!~p8VLbysvADD!Y#B4 z!Q_u4DY5m;r8h! zpG5z^k#97DA(*X(gLi}~U1QE#VOe^MDa`hL3_9;y$PM!1fkXx`JtHVY_ z)ORov1zWhRUgnmvE(97Q#FI#kf+R7L+NS6`g7xE*%@P_e#)@b6xYOu!e1ymhx1mFf zac$*(8S{1M*ESw)uXO?}c3^+>HTLBY<5Xxer#|Nc9U$f{uySX2NJo|zy1^X&b)@e|332nL6jnW#xErGe{(MV|KvRQ(2VS{J0-#e=uaMSkj%>_ z8_o7=<|3!*>B5d*e1WBSK1Aq2gd~tiUiA@IDME!*Al&h-$@PrHsL4E+m{T(5OA&Pw z*S9J~n_hhlD#S_$@Jt)<>;DvXCE!rD-G8iEB2t4SdyOSama+{YJK1+C+l=hl_kAlg zWS6Bx_EH#w49QZBeaV{0nw>!i-|)WwzcHWp{mxw1T-Q9;b$-`5=eh6uInVFReeMUl z5BSNRnzaEL_ws){RG(aIQ=gf2VVk$?Rm4C0Voy;wF13CgTt{bsJo>}_JM&79|Lyrd zETc!a+zQr!`MW;mLFV)LM@^2M4ah>U*508+n$ zh{DDQ1c!&K*|U4tjcoh;OnJi&oSfh&cvrmEykUDE^lJzU~ zMWiuk&Nk@JN9+*G$S%{MW=-w|#p}fOff8Q{^Euw3Y7k}Y@7(0fs?+>3c1k4@?OU?o z6RjRIrVT1(oSxk6Uc)nfyw#iO!bKg9%~gMu;^OW-|Ze(726qs zXSasu1Hq8cz7atvmnSHw<(iPx;lwbh)fyecObqb9PK^Iv`;HcaMI(F19L4kd{z!c1 zEYf9Sg?UpB_i5gt&nW0m<59S7&**?jp!Muk=klQg|EBbYK>N2doATA^aaRd+BuzJZ zRrF>3f7{T8b#tqE)lY;Kk3rwbS)dfD9D!?6{Y)41_vc!We*&AXTbE^*86dBr=^brC z*)?LeuLH@wxaA%d71maDIKBkDjuyN(avxZL+9FWZ#hg_6P>zUIk+~x5s=hRJH6cjJ zP)3_t!)&C6hLxgLRfAXMKH8DoFp0JNV9;oj=JBkISy?hqu=wcfr$u4Ki$aY14~zO1 z{I~BLb1yb{85k{naGpV3|JJ+Av*%yoJ^O3F`nC2l*hn+Bhnc$RX9RPSd736~pPphi zCzP)HH|Esj=!T9jPq@qwG2zQq;Zwg|%u!JxJRL8SY|>rY%7xEpv&ycxKOb3hZ~f7n z`3nw*r{7!;&rN@BV7fQc({iJB&CIjb?vW{>ZAuCQ`CuM#-*j5;K(gIQoqkPe!-7QPsd*X`m+~Af>{63>NE00LTSGc%Lgy-o zWrehwvkB$Wx+c2MK@~l(Fbr?nJ(jei56`P{kzoDTJ-x(=W(@(1Jx*^;UG$Ey1 z9t?U=0Gnw3qG^|}(5K+6UDN{uex{?n*u7TiP5FtPWpv8W5JJ0`c6&HG1Q0sTipbK4 z7HJCkDGfR+PRF@H`>4H#h6B$#Z;qRAvVe?;>dL1qfp#rLt&*Ag?4>iGnuK>cF64Ud z=4KE)69D@hyfi4z4-7Cg+IZ5X$g-J{`es<97k-CilgEB(Hg9V8s)$-6Q-sC$(k_Q6 z!tdY1#yQlBVcE85nV(PS4E7&I4UIunwG-eSZd6?=R9DQ#yT+);cK|$#m0M9B$9OW`lRnU(vMVL zEcW)d*?%ZSewWG!A0Q6xtjYz}zm01KN#rDJ$$?wOeR2(7JHsG;MTxHy8bRUfMu`%i z)?=6_;8`izY8i|mH?Aom`V8K>YnI>`*cB9^mUs0SkpRJPPv&LQ^yYTve!_6#4J3I| ztV2D8{*F+oIx{6dv5Llx$Bi^q`T~eNac9nFY{al|1s9k?kUNic*RCd|e2Hd`Z!2U8 z?D3Quy(~Pm;y=3-0+R#zkwVj&(#^wdXny6u;LY0#k$ubVR&nv^uO>w}5t_=wg;JC# z@8A+ns8_33JX>2n1Uj&G(avIEYMOZAk;|}FE6Q~hwdO14p?o<_zC&kHY-ptKx~ghZ zCPD{~RPFOJi1+g~@39TZZ(Ob#(QuwNUEgzqNz18ZA_bIOrKQX*sho7`QST8Bhym6x zvV@u2?rT|dy9_EyYK^^R$dsD%XE~GPM_&Me1AY45+y{9mjn^tQl1U@Bn_ks+)3!Z^CjTC=z$b6LktM!ow z)ZP2otxm7Mj|y`>+j(Y16`*xj+;z^gw>p1|3ka(gs3g$e;Ay%6D?w2V>yIF&fAPQj z?X=M%?d(P|WxT1Us!NWjLl$<3I`%Tq|H^qL>jp8b8!IgsZXl8>Mwy2iInWdSV6)&e zeD`KS`_1KmHV@L)j;&Y01+K$A3P2y85GUADXg4XUi{8r70RpvoUcKq2!N<3*M3MPh zLsW6)rK$9OopjO`_{}SwIix|>KxXjzM`bD(Ji(cLZb};?mH6BTFJ#^tAyB$|^7;&e zz8|MT$eNygU=K+HYv;Owi%@Jc_tx*VPtGlIyx1Y=4Uz`JN*$224}%_}7V3hRqlko84zjaWZHV zRo&%~G3@`k;h1seA#zlQ#BZ=?jDe)CiC?2WXm?@NN*ccMnsC-5SpWdXL%YX8zLsy; zRvTsM-G%B)0M*Vs68T>40V+YVNH!Z^Afq%=SOf-gB1B18I)t95WZen4gv5~rXYwS8XUjnmUyEqfLuo-o%7cG8e{ z_DTT!Ga-bTL~mDB-f}A%;x97s5SH(qZ;j|q0tGu46W7oSm5h872kT@d4P|Bq^duoI ztX~b>T6#mVbKB=JSMOH_lwehysn&Tma9M@?%SM9!n_>u6(Lhr)0|BHT=^Y6w5RLvy z3xa^^*EBt*!Ix8fmVfD~Fi_Uu9)T@|ZnC(SeeV{!0 z`J!oIpToY4Ut0IX{hPN_XL~+j>cG=N&bDb)drY&c1Z8xI_zEEcchYBRredPjrsceB zU)D-$)&1#q)hwGVeEQb#-jK;P$S41XCvA-Px|@t-#nEsla~Y9u?q^AhcB%Q0bUSx` zG94V{tsSKiFar_*An#vyRR2Tvc9fk|B$zr5GhYZ_-r{Cm5*R2eGZijKdlko3(cSB= zJ~uiN{0{xW0n-{nSI(@|diF)3VVzkSxk>L6j&{DnOIfsiWcmr+ZV1aKL}~#! z+u1V$KPVIr(6?5+Uz?HiDBqedIUDtMCworwNqBw4Lgd@4a$z1wG3m)HO>tV=iC}!G zp>|oPhdpOjTFe)=b}gkCRp_EP5y-k~#<6``XGae*a}QlCD`Zz zy_&w<@aajN$(^**wn)l7KjYTGPaACxEwtkzXIZn^-*!!Fv@>;Sg!8K#Myg2|jU(d| z<0^kXxETuitiKeD2ajf(u-kQ5C@E#%KD0Uc`l+Cixo@2YvzKJ8Ad$c&)ei zh)u#PzqGA)`>e-k=_fF!vJC}(%t)+^_c9q07A_umE8IEf`6mo$1?xEh?CIo@l*Ft&!Ir{ayn>6MOeDQuF2zH`^?|y$NS@jc-;% z{ssd+i>&3i&-<2~AnW9+Ah49csGd`65^6{MYFFSA*G}=`^zC^H(TwQu6!b$5vLDq0 z%C0Q!2KPr20cLvaby1cwALO?wku5!;JEBVcR1b_uG6D+fHy9lV?H^DU#8SrYQ*g-r zG{5u;mZ?(sEco*b=y%8H&)b(*#DAsu$5ZOxApesL_gvpFOK2ImE>>mL(SSL;Yj+Bk z=&a#{%1DTj)KK`2R=)KVu8v;a5Z8Fg;nV=VOi9yy1DeMCa9lRBDVn!Lgc5=XTWVSR zLR+gq98~X(InWx+&+7V!Q97mlvZwx%`Fw5CXS6UXb+TfW$xw6rcIC$RL1}?(v^qGD zSBgs2ST>7Dht@B!=2D?F6u-A(up8l>L{X$Y73Ii7 zTwAc4ZX7nvJC-4ZBNP|h7J9e{*v%_WhNoQn2;u)5R{TYSVYjU~;)8Gjv0EaXbgWe8 zBp@~vE(&(5fs-f&dZ(hC8XoMe_J?f=jwMruT7dxz49*dYScKzn>8YZqfb}8tvcy=snepu>%m(p;xu}dH) z-3A9>Y5zqDA3lt50kCQOlK`QCxK#3~vjm%$J~<5Vf^h+kxqnOpcgPjbWG=n$}7vRW@Q&F%VjVH}mia#}t(oi8HKD?_8 Pn2+=s0N`rk;jjMy>V literal 0 HcmV?d00001 diff --git a/data/gem/Global-Cement-and-Concrete-Tracker_July-2025.xlsx b/data/gem/Global-Cement-and-Concrete-Tracker_July-2025.xlsx new file mode 100644 index 0000000000000000000000000000000000000000..c2a64bb46cf1b989bc7116b0a7d7d5366f8b24c7 GIT binary patch literal 972391 zcmeEv2{@GN`#(`=Q_-T*L?mS?qzzM&eTmaKm8EGDij{v3X8z~r`+d8va~)lccb@lNp69+l_vgN!^;Qy+ z@?z4`(qe1fC7r~+`7;Cj?C+_)(`TEX$6mJsUw>GlKt9VmzATiG^?sGQ@Jx}OXi))*vJw`%U=vvHtLN3NXnVe zGB+c3_H@?^GB+=klzfD2`-rj0)1zHXAv(fAGClt!^o)|Roln<-*3S>q7(JS!K=VS`hR)FJn(jYpT8=OdjdboXfEL4Vp zxOg0=fEEH(qPc7+q>obtwWB#`s2fG$G(o{oHYK>P4A0qjgo1;j&=48T?V)hF(}@i3 zGKg83N91y81XiAR2vmX^N}zGOJF_81$ug*$i0k1HIZeqJs1MW7-$@JUVN_rU+z)ja zXpoMEI;L{DvQWyNV7fNcTSp0&!_kw`SXLzsVvr-)P017vGYTEdra=RF&>$yjBWn=- z;C&vEQ;8z5xCOYOW(udRz#BT0jxu`W9nw>YKlYID>vfPBX-*XlGJ&ATD?r5-J3u;|iWY4C)T{kRc8#p+j^C2q`4E2WMme z(T$d&&!LEn`&cFv>^}tgrH5g-ADWczzv|28+-A^lH!x6J<{~VMQt;ePtv4^6O+vF< zV>u`;6)cm4!qMv}I1<$4MuDhxsE}@mmQ5I{_k()CcW=9oJJyp2u`m$BjZuo`lxT+p zQZR%~I94St7&pXp!(Hfsm|X=>U%R(F6dXhvxJSiIWQWoEmf@J(C<- zH`wWoW9;ef38O$vG}Oz27+NQ z>JFoEo854XJ}2D$20W30!`**M$8iU>xnqc9mBuPfRFYAepfpYiT4cbLTqo@+{j>Bz zX&Y%TX%}g%^jhg1(mzQbkhYTEBkd#|ByA?$K02hQM4js!AB9$I-QgdPJC15&lyC`9 zlPkovBeKIVY!=?Fr+JCNiLyz?3dU276^-SLXBkg7Mioygo>44cJiB;G@yz0B#dC@$ z7Y{K?ml$w-$J;e6VNboFUDiEmcF9lvS|r8m| zLT&Jc9pDRR^1iUM57npiSV=`GO=+i6OLmj}iS;MeA73B#GW6w%m&ad*Ee>6LV)60C zVYfnWow#-UR@lDKeT0zx@!mjWPX6)U;PaEhqC=xkL?4eP9ye3?XE6zO&bl(OZ(ig9Ts3)4(zrmb63e1bFR!mKhu12!k&F(_mOz>mqKZH!E2y%H~Y`u2zB>*Fp?}+oGhHWl=twb7i%g6hhe3-JVvixFq=l zxml`3;i)Y#zP+jMUo5`ADk>znx^sL<@?*16&r!yAx|H(~Nhgt#M3S8@O||=}zNaJJx~V-L8iwd^~ibsmf!H#Y*$=4d$oy%p-iwqsBiyY5g=(h+7CV*|dGD;3s>4E82d;U|j}{&KqS>t`K!wV-+3Pl; z)ImAjVSa0risJXRz+Ol-2$`Q^27VD2;|q3?n^+ahmGc;kIj#!nV)De|GtaZVUDL}R_*JT&p3tQ zEFkYV@x^a_jMo0>cSdJqRuGls=br|@li)TfZIfl%hG%K#tJeQ$VM?ECM z_@0X3J&zfZ9trlI3%6F`M}JNnCE$Cq0=b@;svK>>-t!%~&y5wBfmND3Wu8o|vO?^< z(ckxb+-t%XqdM}q&T-+3_j{X5%|M&I}5)i!*;_ z1#&zw*@BrCIr|;C1DIC2K&C|q0j7nnBQFw6VRpb3T9f^*MtVJYt0hOv^s%l4ZbBb2)1i^1Je?z$T>jdB*8>Z+>N0Dk*5eI zaxM@#m%`_l#I%V);TDhaZr8jFz=p&+Ca_R(={rUHOg5v^E+#lQkIP!)4fkYY@79|-;xdIg=r4g z(j0fE{qiYo%gi+Uq_oYBX+OV9b6J_zXZlENj<&?bM2nc;AO)-A8!XRkusZ$9GUb)k zF+IyTJ*$Wu%j6uZFyDp+Gm1#AVl8^=zd;s0SEDVj#_qlv<9IbL&?DN^BX*ZZ%r72s zhp5pnW{sBcCyYN_Gk=|f20wSzn(CnbvxDYRheZb+v?tUou&L3|tXb$)qor6Akh)b& zJ3)L>Q=0jybQ8;T>u2f43)3yHrJL?fxBZm9W@frYQu;c_benhSYgc;9k8d2%5o4a0 zx40c|el_0W!3Oj64HoxanP%SI^G8dS>O-gfSlR^F7W^^EjjLkuZ^3XutR{y+Z6tprAtDuORlwP9P!E z#SNS-8;_|@66n*9#I#hqA5D=KP;v=V7+YmR1$;p(nc!QJ#sVrkUD9aNs7nIs#9