From 73601f970b7146670b13174706807751a1c40926 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 19 Nov 2023 16:58:24 +0300 Subject: [PATCH 01/50] Initial commit --- build_climate_projections.py | 193 +++++++++++++++++++++++++++++++++++ 1 file changed, 193 insertions(+) create mode 100644 build_climate_projections.py diff --git a/build_climate_projections.py b/build_climate_projections.py new file mode 100644 index 000000000..90cecb7c6 --- /dev/null +++ b/build_climate_projections.py @@ -0,0 +1,193 @@ +# -*- coding: utf-8 -*- +# SPDX-FileCopyrightText: PyPSA-Earth and PyPSA-Eur Authors +# +# SPDX-License-Identifier: AGPL-3.0-or-later + +# -*- coding: utf-8 -*- +""" +Creates a dataset corresponding to the projected climate in a requested year. + +Relevant Settings +----------------- + +.. code:: yaml + + scenario: + future_clim + + future_climate: + ssp_scenario: + reference_year: + future_weather_year: + climate_avr_windows: + +Inputs +------ + +- ``cutouts/cutout.nc``: confer :ref:`cutout`, a cutout file PyPSA Network + +Outputs +------- + +- ``cutouts/cutout_future_{future_weather_year}.nc``: A cutout modified to account for future climate conditions + +Description +----------- + +The rule :mod:`build_climate_projections` creates a cutout file which corresponds to a requested year in the future. Temperature projectons are being calculated combining data for cutout.nc which is assumed to be representative of the past climate and an ensemble of CMIP6 globale climate models to account for the future climate +""" + + +import datetime as dt +import os + +import atlite +import geopandas as gpd +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import pypsa +import xarray as xr +from shapely.geometry import LineString as Line +from shapely.geometry import Point + +# TODO export parameters from the config +scenario_name = "ssp245" +month_in_focus = 5 # a temporary parameter + + +cutout_path = "cutouts/sar-2013-era5.nc" +region_name = "sar" +cmip6_ens_dir = "data/cmip6/Global Atlas/tas_" + scenario_name + "/" +cmip6_fl = "t_CMIP6_" + scenario_name + "_mon_201501-210012.nc" + + +def crop_cmip6(cmip6_xr, cutout_xr): + # spartial margin needed to avoid data losses during further interpolation + d_pad = 1 + + cmip6_region = cmip6_xr.sel( + lat=slice( + min(cutout_xr.coords["y"].values) - d_pad, + max(cutout_xr.coords["y"].values) + d_pad, + ), + lon=slice( + min(cutout_xr.coords["x"].values) - d_pad, + max(cutout_xr.coords["x"].values) + d_pad, + ), + ) + return cmip6_region + + +def interpolate_cmip6_to_cutout_grid(cmip6_xr, cutout_xr): + # TODO read from the cutout file instead of hardcoding + dx_new = 0.3 + + newlon = np.arange( + round(min(cutout_xr.coords["x"].values), 1), + round(max(cutout_xr.coords["x"].values) + dx_new, 1), + dx_new, + ) + newlat = np.arange( + round(min(cutout_xr.coords["y"].values), 1), + round(max(cutout_xr.coords["y"].values) + dx_new, 1), + dx_new, + ) + # Interpolate + cmip6_interp = cmip6_xr.interp(time=cmip6_xr.time, lat=newlat, lon=newlon) + + return cmip6_interp + + +# TODO fix years_window +def subset_by_time(cmip6_xr, month, year, years_window=5): + # TODO add warning in case there are no enough years + cmip6_in_period = cmip6_xr.where( + ( + (cmip6_xr["time.month"] == month) + & (cmip6_xr["time.year"] >= year - years_window) + & (cmip6_xr["time.year"] < year + years_window - 1) + ), + drop=True, + ) + return cmip6_in_period + + +def calculate_proj_of_average(cmip6_xr, month, year0, year1, years_window=5): + cmip6_interp_year0 = subset_by_time( + cmip6_xr, + month=month_in_focus, + year=year0, + ) + cmip6_interp_year1 = subset_by_time(cmip6_xr, month=month_in_focus, year=year1) + dt_interp = cmip6_interp_year1["t"].mean("member").mean( + "time" + ) - cmip6_interp_year0["t"].mean("member").mean("time") + return dt_interp + + +def build_projection_for_month(cutout_xr, dt_xr, month): + k_time = cutout_xr.time.dt.month.isin(month) + + for i in range(0, len(cutout_xr.y)): + for j in range(0, len(cutout_xr.x)): + cutout_xr.temperature[k_time, i, j] = np.add( + cutout_xr.temperature[k_time, i, j], + np.array([dt_xr[i, j].item()] * k_time.sum().item()), + ) + + return cutout_xr + + +def build_cutout_future( + cutout_xr, cmip6_xr, months=5, year0=2020, year1=2070, years_window=5 +): + for k_month in [months]: + dt_k = calculate_proj_of_average( + cmip6_xr=cmip6_xr, month=k_month, year0=year0, year1=year1, years_window=5 + ) + + build_projection_for_month(cutout_xr, dt_k, k_month) + + cutout_xr = cutout_xr.where(cutout_xr.time.dt.month.isin(months), drop=True) + + return cutout_xr + + +cmip6_xr = xr.open_dataset(os.path.join(cmip6_ens_dir, cmip6_fl)) + +cutout_xr = xr.open_dataset(cutout_path) + +cmip6_region = crop_cmip6(cmip6_xr, cutout_xr) + +cmip6_region_interp = interpolate_cmip6_to_cutout_grid( + cmip6_xr=cmip6_region, cutout_xr=cutout_xr +) + +# graphical test of interpolation +fig = plt.figure() +cmip6_region_interp["t"].mean("time").mean("member").plot() +fig.savefig("test_cmip6_interp.png", dpi=700) + +# test of projection +dt_interp_2070vs2020 = calculate_proj_of_average( + cmip6_xr=cmip6_region_interp, + month=5, + year0=2020, + year1=2070, + years_window=month_in_focus, +) +print("dt_interp_2070vs2020") +print(dt_interp_2070vs2020) + +cutout_future = build_cutout_future( + cutout_xr=cutout_xr, + cmip6_xr=cmip6_region_interp, + months=5, + year0=2020, + year1=2070, + years_window=5, +) +cutout_future.to_netcdf( + "cutout_2070_" + str(month_in_focus) + "_" + region_name + ".nc" +) From c68abbea29967299b550541e7b60d8925dcd8c11 Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 28 Nov 2023 22:20:32 +0300 Subject: [PATCH 02/50] Move a projection script into script directory --- .../build_climate_projections.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename build_climate_projections.py => scripts/build_climate_projections.py (100%) diff --git a/build_climate_projections.py b/scripts/build_climate_projections.py similarity index 100% rename from build_climate_projections.py rename to scripts/build_climate_projections.py From 2beb87954b2ca0d3f70cd79e578422dbff5212b2 Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 29 Nov 2023 16:34:30 +0300 Subject: [PATCH 03/50] Restructure into Snakemake workflow --- Snakefile | 23 ++++++ scripts/build_climate_projections.py | 105 ++++++++++++++++----------- 2 files changed, 86 insertions(+), 42 deletions(-) diff --git a/Snakefile b/Snakefile index 24899f6b6..40107653b 100644 --- a/Snakefile +++ b/Snakefile @@ -362,6 +362,29 @@ if config["enable"].get("build_cutout", False): "scripts/build_cutout.py" +if config["enable"].get("modify_cutout", False): + + rule build_climate_projections: + params: + snapshots=config["snapshots"], + climate_scenario=config["projection"]["climate_scenario"], + future_year=config["projection"]["future_year"], + years_window=config["projection"]["years_window"], + input: + "cutouts/" + CDIR + "{cutout}.nc", + output: + "cutouts/" + CDIR + "{cutout}_{future_year}.nc", + log: + "logs/" + RDIR + "build_climate_projections/{cutout}_{future_year}.log", + benchmark: + "benchmarks/" + RDIR + "build_climate_projections_{cutout}_{future_year}" + threads: ATLITE_NPROCESSES + resources: + mem_mb=ATLITE_NPROCESSES * 1000, + script: + "scripts/build_climate_projections.py" + + if config["enable"].get("build_natura_raster", False): rule build_natura_raster: diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 90cecb7c6..9e0a5b214 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -51,11 +51,6 @@ from shapely.geometry import LineString as Line from shapely.geometry import Point -# TODO export parameters from the config -scenario_name = "ssp245" -month_in_focus = 5 # a temporary parameter - - cutout_path = "cutouts/sar-2013-era5.nc" region_name = "sar" cmip6_ens_dir = "data/cmip6/Global Atlas/tas_" + scenario_name + "/" @@ -154,40 +149,66 @@ def build_cutout_future( return cutout_xr -cmip6_xr = xr.open_dataset(os.path.join(cmip6_ens_dir, cmip6_fl)) - -cutout_xr = xr.open_dataset(cutout_path) - -cmip6_region = crop_cmip6(cmip6_xr, cutout_xr) - -cmip6_region_interp = interpolate_cmip6_to_cutout_grid( - cmip6_xr=cmip6_region, cutout_xr=cutout_xr -) - -# graphical test of interpolation -fig = plt.figure() -cmip6_region_interp["t"].mean("time").mean("member").plot() -fig.savefig("test_cmip6_interp.png", dpi=700) - -# test of projection -dt_interp_2070vs2020 = calculate_proj_of_average( - cmip6_xr=cmip6_region_interp, - month=5, - year0=2020, - year1=2070, - years_window=month_in_focus, -) -print("dt_interp_2070vs2020") -print(dt_interp_2070vs2020) - -cutout_future = build_cutout_future( - cutout_xr=cutout_xr, - cmip6_xr=cmip6_region_interp, - months=5, - year0=2020, - year1=2070, - years_window=5, -) -cutout_future.to_netcdf( - "cutout_2070_" + str(month_in_focus) + "_" + region_name + ".nc" -) +if __name__ == "__main__": + if "snakemake" not in globals(): + from _helpers import mock_snakemake + + os.chdir(os.path.dirname(os.path.abspath(__file__))) + snakemake = mock_snakemake( + "build_climate_projections", + snapshots=config["snapshots"], + climate_scenario=config["projection"]["climate_scenario"], + future_year=config["projection"]["future_year"], + years_window=config["projection"]["years_window"], + cutout="africa-2013-era5", + ) + configure_logging(snakemake) + + cutout_params = snakemake.params.cutouts[snakemake.wildcards.cutout] + snapshots = pd.date_range(freq="h", **snakemake.params.snapshots) + time = [snapshots[0], snapshots[-1]] + cutout_params["time"] = slice(*cutout_params.get("time", time)) + onshore_shapes = snakemake.input.onshore_shapes + offshore_shapes = snakemake.input.offshore_shapes + + # TODO export parameters from the config + scenario_name = "ssp245" + month_in_focus = 5 # a temporary parameter + + cmip6_xr = xr.open_dataset(os.path.join(cmip6_ens_dir, cmip6_fl)) + + cutout_xr = xr.open_dataset(cutout_path) + + cmip6_region = crop_cmip6(cmip6_xr, cutout_xr) + + cmip6_region_interp = interpolate_cmip6_to_cutout_grid( + cmip6_xr=cmip6_region, cutout_xr=cutout_xr + ) + + # graphical test of interpolation + fig = plt.figure() + cmip6_region_interp["t"].mean("time").mean("member").plot() + fig.savefig("test_cmip6_interp.png", dpi=700) + + # test of projection + dt_interp_2070vs2020 = calculate_proj_of_average( + cmip6_xr=cmip6_region_interp, + month=5, + year0=2020, + year1=2070, + years_window=month_in_focus, + ) + print("dt_interp_2070vs2020") + print(dt_interp_2070vs2020) + + cutout_future = build_cutout_future( + cutout_xr=cutout_xr, + cmip6_xr=cmip6_region_interp, + months=5, + year0=2020, + year1=2070, + years_window=5, + ) + cutout_future.to_netcdf( + "cutout_2070_" + str(month_in_focus) + "_" + region_name + ".nc" + ) From 1a24bad4b35296ec46daec1269af4c6009d7fe8c Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 30 Nov 2023 00:53:14 +0300 Subject: [PATCH 04/50] Organize code --- scripts/build_climate_projections.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 9e0a5b214..46cdf3b0c 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -48,13 +48,11 @@ import pandas as pd import pypsa import xarray as xr +from _helpers import configure_logging, create_logger from shapely.geometry import LineString as Line from shapely.geometry import Point -cutout_path = "cutouts/sar-2013-era5.nc" -region_name = "sar" -cmip6_ens_dir = "data/cmip6/Global Atlas/tas_" + scenario_name + "/" -cmip6_fl = "t_CMIP6_" + scenario_name + "_mon_201501-210012.nc" +logger = create_logger(__name__) def crop_cmip6(cmip6_xr, cutout_xr): From 174b0b15179ebf5847ac426f567dd3a60846b8fb Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 30 Nov 2023 18:21:41 +0300 Subject: [PATCH 05/50] Add inputs --- scripts/build_climate_projections.py | 30 ++++++++++++---------------- 1 file changed, 13 insertions(+), 17 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 46cdf3b0c..f154615b0 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -154,29 +154,25 @@ def build_cutout_future( os.chdir(os.path.dirname(os.path.abspath(__file__))) snakemake = mock_snakemake( "build_climate_projections", - snapshots=config["snapshots"], - climate_scenario=config["projection"]["climate_scenario"], - future_year=config["projection"]["future_year"], - years_window=config["projection"]["years_window"], + climate_scenario="ssp2-2.6", + future_year=2050, + years_window=5, cutout="africa-2013-era5", ) configure_logging(snakemake) - cutout_params = snakemake.params.cutouts[snakemake.wildcards.cutout] - snapshots = pd.date_range(freq="h", **snakemake.params.snapshots) - time = [snapshots[0], snapshots[-1]] - cutout_params["time"] = slice(*cutout_params.get("time", time)) - onshore_shapes = snakemake.input.onshore_shapes - offshore_shapes = snakemake.input.offshore_shapes - - # TODO export parameters from the config - scenario_name = "ssp245" - month_in_focus = 5 # a temporary parameter + climate_scenario = snakemake.params.climate_scenario + present_year = snakemake.params.present_year + future_year = snakemake.params.future_year + years_window = snakemake.params.years_window + cutout = snakemake.input.cutout + cmip6 = snakemake.input.cmip6_avr - cmip6_xr = xr.open_dataset(os.path.join(cmip6_ens_dir, cmip6_fl)) - - cutout_xr = xr.open_dataset(cutout_path) + snapshots = pd.date_range(freq="h", **snakemake.params.snapshots) + season_in_focus = snapshots.month.unique().to_list() + cmip6_xr = xr.open_dataset(cmip6) + cutout_xr = xr.open_dataset(cutout) cmip6_region = crop_cmip6(cmip6_xr, cutout_xr) cmip6_region_interp = interpolate_cmip6_to_cutout_grid( From 4594541f565a9ec68c075f4b62256e1482173945 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 30 Nov 2023 18:22:10 +0300 Subject: [PATCH 06/50] Fix Snakemake inputs --- Snakefile | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/Snakefile b/Snakefile index 40107653b..f4b0c5327 100644 --- a/Snakefile +++ b/Snakefile @@ -368,10 +368,12 @@ if config["enable"].get("modify_cutout", False): params: snapshots=config["snapshots"], climate_scenario=config["projection"]["climate_scenario"], + present_year=config["projection"]["present_year"], future_year=config["projection"]["future_year"], years_window=config["projection"]["years_window"], input: - "cutouts/" + CDIR + "{cutout}.nc", + cutout="cutouts/" + CDIR + "{cutout}.nc", + cmip6_avr="/Users/ekaterina/Documents/_github_/cmip6/Global Atlas/tas_ssp245/t_CMIP6_ssp245_mon_201501-210012.nc", output: "cutouts/" + CDIR + "{cutout}_{future_year}.nc", log: From 83c6b0b76d15cbbfe2ff0bf93d042fc02820bab5 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 30 Nov 2023 18:27:09 +0300 Subject: [PATCH 07/50] Improve documenting --- scripts/build_climate_projections.py | 19 +++++++++---------- 1 file changed, 9 insertions(+), 10 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index f154615b0..e04bf7ec3 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -12,24 +12,23 @@ .. code:: yaml - scenario: - future_clim - - future_climate: - ssp_scenario: - reference_year: - future_weather_year: - climate_avr_windows: + projection: + climate_scenario: + present_year: + future_year: + years_window: + gcm_selection: Inputs ------ -- ``cutouts/cutout.nc``: confer :ref:`cutout`, a cutout file PyPSA Network +- ``cutouts/cutout.nc``: confer :ref:`cutout`, a cutout file produced by altile +- ``data/cmip6_avr.cn``: Outputs ------- -- ``cutouts/cutout_future_{future_weather_year}.nc``: A cutout modified to account for future climate conditions +- ``cutouts/{cutout}_{future_year}.nc"``: A cutout modified to account for future climate conditions Description ----------- From b196e5a75e0804e5d5ec065f28d2362ca502b186 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 30 Nov 2023 18:27:56 +0300 Subject: [PATCH 08/50] Code clean-up --- scripts/build_climate_projections.py | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index e04bf7ec3..46bf40a1e 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -85,14 +85,13 @@ def interpolate_cmip6_to_cutout_grid(cmip6_xr, cutout_xr): round(max(cutout_xr.coords["y"].values) + dx_new, 1), dx_new, ) - # Interpolate cmip6_interp = cmip6_xr.interp(time=cmip6_xr.time, lat=newlat, lon=newlon) return cmip6_interp # TODO fix years_window -def subset_by_time(cmip6_xr, month, year, years_window=5): +def subset_by_time(cmip6_xr, month, year, years_window): # TODO add warning in case there are no enough years cmip6_in_period = cmip6_xr.where( ( @@ -105,13 +104,13 @@ def subset_by_time(cmip6_xr, month, year, years_window=5): return cmip6_in_period -def calculate_proj_of_average(cmip6_xr, month, year0, year1, years_window=5): +def calculate_proj_of_average(cmip6_xr, month, year0, year1, years_window): cmip6_interp_year0 = subset_by_time( cmip6_xr, - month=month_in_focus, + month, year=year0, ) - cmip6_interp_year1 = subset_by_time(cmip6_xr, month=month_in_focus, year=year1) + cmip6_interp_year1 = subset_by_time(cmip6_xr, month=month, year=year1) dt_interp = cmip6_interp_year1["t"].mean("member").mean( "time" ) - cmip6_interp_year0["t"].mean("member").mean("time") @@ -131,9 +130,7 @@ def build_projection_for_month(cutout_xr, dt_xr, month): return cutout_xr -def build_cutout_future( - cutout_xr, cmip6_xr, months=5, year0=2020, year1=2070, years_window=5 -): +def build_cutout_future(cutout_xr, cmip6_xr, months, year0, year1, years_window): for k_month in [months]: dt_k = calculate_proj_of_average( cmip6_xr=cmip6_xr, month=k_month, year0=year0, year1=year1, years_window=5 @@ -154,6 +151,7 @@ def build_cutout_future( snakemake = mock_snakemake( "build_climate_projections", climate_scenario="ssp2-2.6", + present_year=2020, future_year=2050, years_window=5, cutout="africa-2013-era5", From 369bc9bd93a4596eed67e1ff1715782050844eee Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 30 Nov 2023 18:35:23 +0300 Subject: [PATCH 09/50] Remove hardcoding and add debug plotting --- scripts/build_climate_projections.py | 39 ++++++++++++++++------------ 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 46bf40a1e..b381e4089 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -176,30 +176,35 @@ def build_cutout_future(cutout_xr, cmip6_xr, months, year0, year1, years_window) cmip6_xr=cmip6_region, cutout_xr=cutout_xr ) + # ----------------------------------------------------------------- + # to be replaced after debug + # ----------------------------------------------------------------- # graphical test of interpolation fig = plt.figure() cmip6_region_interp["t"].mean("time").mean("member").plot() fig.savefig("test_cmip6_interp.png", dpi=700) - # test of projection - dt_interp_2070vs2020 = calculate_proj_of_average( - cmip6_xr=cmip6_region_interp, - month=5, - year0=2020, - year1=2070, - years_window=month_in_focus, - ) - print("dt_interp_2070vs2020") - print(dt_interp_2070vs2020) + fig = plt.figure() + (cutout_xr["temperature"] - 273.15).mean("time").plot(cmap="OrRd") + fig.savefig("results/test_present.png", dpi=700) + # ----------------------------------------------------------------- + # TODO Add a check for time match cutout_future = build_cutout_future( cutout_xr=cutout_xr, cmip6_xr=cmip6_region_interp, - months=5, - year0=2020, - year1=2070, - years_window=5, - ) - cutout_future.to_netcdf( - "cutout_2070_" + str(month_in_focus) + "_" + region_name + ".nc" + months=season_in_focus, + year0=present_year, + year1=future_year, + years_window=years_window, ) + + cutout_future.to_netcdf(snakemake.output[0]) + + # ----------------------------------------------------------------- + # to be replaced after debug + # ----------------------------------------------------------------- + # graphical test of projection + fig = plt.figure() + (cutout_future["temperature"] - 273.15).mean("time").plot(cmap="OrRd") + fig.savefig("results/test_cmip6_project.png", dpi=700) From d1fd454f185c67dafc6c25bd42db7a30d126f042 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 30 Nov 2023 18:37:10 +0300 Subject: [PATCH 10/50] Add parameters to the tacked configs --- config.default.yaml | 8 ++++++++ config.tutorial.yaml | 9 +++++++++ 2 files changed, 17 insertions(+) diff --git a/config.default.yaml b/config.default.yaml index ba638576a..22fbde92a 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -204,6 +204,14 @@ atlite: # y: [33., 72] # manual set cutout range +projection: + climate_scenario: ssp245 + present_year: 2020 + future_year: 2050 + years_window: 5 # a half of width currently + gcm_selection: false # false to use all the available global climate models + + renewable: onwind: cutout: cutout-2013-era5 diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 410d93683..527c9c19a 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -218,6 +218,15 @@ atlite: # x: [-12., 35.] # set cutout range manual, instead of automatic by boundaries of country # y: [33., 72] # manual set cutout range + +projection: + climate_scenario: ssp245 + present_year: 2020 + future_year: 2050 + years_window: 5 # a half of width currently + gcm_selection: false # false to use all the available global climate models + + renewable: onwind: cutout: cutout-2013-era5-tutorial From 72272d5cd14d6c13bf805fd1de99f48b1aa13af9 Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 26 Dec 2023 18:30:26 +0300 Subject: [PATCH 11/50] Fix wildcards for the projection --- Snakefile | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/Snakefile b/Snakefile index 3b4214abe..5c7d703a0 100644 --- a/Snakefile +++ b/Snakefile @@ -57,6 +57,8 @@ wildcard_constraints: ll="(v|c)([0-9\.]+|opt|all)|all", opts="[-+a-zA-Z0-9\.]*", unc="[-+a-zA-Z0-9\.]*", + # "^proj-" is reserved for projections, negative lookbehind assertion used in the regexp + cutout=".+(? Date: Tue, 26 Dec 2023 23:09:49 +0300 Subject: [PATCH 12/50] Adjust inputs-outputs in the script docstring --- scripts/build_climate_projections.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index b381e4089..8abefa00b 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -23,12 +23,14 @@ ------ - ``cutouts/cutout.nc``: confer :ref:`cutout`, a cutout file produced by altile -- ``data/cmip6_avr.cn``: +- ``data/cmip6_avr.nc``: confer :ref:`cmip6`, a global-scale dataset of the global climate +CMIP6 projections aggregated into IPCC Global Climate Atlas Outputs ------- -- ``cutouts/{cutout}_{future_year}.nc"``: A cutout modified to account for future climate conditions +- ``cutouts/proj-{cutout}.nc"``: confer :ref:`proj`, a cutout modified to account for future +climate conditions Description ----------- From d77af95cef82a6ec170096170a07aefbcefea68e Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 26 Dec 2023 23:11:20 +0300 Subject: [PATCH 13/50] Add methodology --- scripts/build_climate_projections.py | 26 +++++++++++++++++++++++++- 1 file changed, 25 insertions(+), 1 deletion(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 8abefa00b..db2695e52 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -35,7 +35,31 @@ Description ----------- -The rule :mod:`build_climate_projections` creates a cutout file which corresponds to a requested year in the future. Temperature projectons are being calculated combining data for cutout.nc which is assumed to be representative of the past climate and an ensemble of CMIP6 globale climate models to account for the future climate +The rule :mod:`build_climate_projections` creates a cutout file which corresponds to +a requested year in the future. The parameter projections is calculated combining data +for cutout.nc which is assumed to be representative data sample for the past climate and +an ensemble of CMIP6 global climate models to account for the future climate. + +A morphing approach is implemented meaning that the projection value x_future is calculated +data for a representative period present_year in the past x_past using a combination +of shift dx and stretch a: + +x_future = x_past + dx + a(x_past - |month), + +where +x_future and x_past have the highest temporal resolution, +|month is the monthly average for a considered month. + +The procedure is applied for each month to account for seasonality. + +Methodology notes +----------------- +Originally the morphing approach originates from buildings physics, and is still in active +use for detailed simulation of heat transfer in buildings. See, for example: + +J.B. Dias, G.C. da Graça, P.M.M. Soares (2020) Comparison of methodologies for +generation of future building thermal energy simulation. Energy and Buildings, +vol. 206,2020, 109556, https://doi.org/10.1016/j.enbuild.2019.109556. """ From 51c8a099341648e9f6db5c3df0d318e61d7ee063 Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 26 Dec 2023 23:12:40 +0300 Subject: [PATCH 14/50] Add docstrings --- scripts/build_climate_projections.py | 119 +++++++++++++++++++++++++++ 1 file changed, 119 insertions(+) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index db2695e52..fad9bfa06 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -81,6 +81,17 @@ def crop_cmip6(cmip6_xr, cutout_xr): + """ + Crop the global-scale CMIP6 dataset to match the cutout extend. + + Parameters + ---------- + cmip6_xr: xarray + CMIP6 climate projections loaded as xarray + cutout_xr : xarray of a cutout dataset + Cutout dataset loaded as xarray + """ + # spartial margin needed to avoid data losses during further interpolation d_pad = 1 @@ -98,6 +109,17 @@ def crop_cmip6(cmip6_xr, cutout_xr): def interpolate_cmip6_to_cutout_grid(cmip6_xr, cutout_xr): + """ + Interpolate CMIP6 dataset to the cutout grid. + + Parameters + ---------- + cmip6_xr: xarray + CMIP6 climate projections loaded as xarray + cutout_xr : xarray of a cutout dataset + Cutout dataset loaded as xarray + """ + # TODO read from the cutout file instead of hardcoding dx_new = 0.3 @@ -118,6 +140,27 @@ def interpolate_cmip6_to_cutout_grid(cmip6_xr, cutout_xr): # TODO fix years_window def subset_by_time(cmip6_xr, month, year, years_window): + """ + Filter CMIP6 dataset according to the requested year and the processed + month. + + Parameters + ---------- + cmip6_xr: xarray + CMIP6 climate projections loaded as xarray + month: integer + A month value to be considered further for the morphing procedure + year: integer + The first year of the requested period in the future + years_window: integer + A width of the considered time period + + Example + ------- + month=3, year=2050, years_window=30 + filters CMIP6 cmip6_xr dataset to calculate projection for March at 2050-2079 + """ + # TODO add warning in case there are no enough years cmip6_in_period = cmip6_xr.where( ( @@ -131,6 +174,36 @@ def subset_by_time(cmip6_xr, month, year, years_window): def calculate_proj_of_average(cmip6_xr, month, year0, year1, years_window): + """ + Calculate a change in a monthly average value for a climate parameter for + the requested year periods and the processed month (a single value for each + grid cell) + + Parameters + ---------- + cmip6_xr: xarray + CMIP6 climate projections loaded as xarray + month: float + A month value to be considered further for the morphing procedure + year0: integer + The first year of an earlier period in the future + year1: integer + The first year of a later period in the future + years_window: integer + Width of the considered time period + + Outputs + ------- + dt_interp: xarray of the changes in the monthly averages for each grid cell + Contains a single value for each grid cell. + + Example + ------- + month=3, year0=2020, year0=2070, years_window=30 + calculates a multimodel change in the monthly average for + March in 2020-2049 as compared with 2070-2099 + """ + cmip6_interp_year0 = subset_by_time( cmip6_xr, month, @@ -144,6 +217,26 @@ def calculate_proj_of_average(cmip6_xr, month, year0, year1, years_window): def build_projection_for_month(cutout_xr, dt_xr, month): + """ + Calculate the projection for the cutout data applying the morphing approach + using a dataset dt_xr for the changes in monthly aggregates for the + processed month. + + Parameters + ---------- + cutout_xr: xarray + Cutout dataset loaded as xarray + dt_xr: xarray + Dataset of the changes in the monthly averages for each grid cell + month: float + A month value to be considered further for the morphing procedure + + Output + ------- + cutout_xr: xarray of the projected changes for each grid cell for a processed month + Contains a time series for each grid cell. + """ + k_time = cutout_xr.time.dt.month.isin(month) for i in range(0, len(cutout_xr.y)): @@ -157,6 +250,32 @@ def build_projection_for_month(cutout_xr, dt_xr, month): def build_cutout_future(cutout_xr, cmip6_xr, months, year0, year1, years_window): + """ + Calculate the projection for the cutout data applying the morphing approach + for each month in the months. + + Parameters + ---------- + cutout_xr: xarray + Cutout dataset loaded as xarray + cmip6_xr: xarray + CMIP6 climate projections loaded as xarray + months: list + Months values to be used as a definition of a season to be considered + for building projection time series + year0: integer + The first year of an earlier period in the future + year1: integer + The first year of a later period in the future + years_window: integer + Width of the considered time period + + Output + ------- + cutout_xr: xarray + Dataset on the projected changes for each grid cell for a processed month + Contains a time series for each grid cell. + """ for k_month in [months]: dt_k = calculate_proj_of_average( cmip6_xr=cmip6_xr, month=k_month, year0=year0, year1=year1, years_window=5 From cebf6a338c6351b99590c67a7c6a2c5221b83559 Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 26 Dec 2023 23:12:49 +0300 Subject: [PATCH 15/50] Add TODOs --- scripts/build_climate_projections.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index fad9bfa06..bb4e6850d 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -173,6 +173,8 @@ def subset_by_time(cmip6_xr, month, year, years_window): return cmip6_in_period +# TODO fix hardcoding in the parameter name +# TODO add functionality to customise CMIP6 ensemble def calculate_proj_of_average(cmip6_xr, month, year0, year1, years_window): """ Calculate a change in a monthly average value for a climate parameter for @@ -216,6 +218,7 @@ def calculate_proj_of_average(cmip6_xr, month, year0, year1, years_window): return dt_interp +# TODO fix hardcoding in the parameter name def build_projection_for_month(cutout_xr, dt_xr, month): """ Calculate the projection for the cutout data applying the morphing approach From 0d7ec1d12f7cceedaa647d0d9e45b8fc360e31a5 Mon Sep 17 00:00:00 2001 From: ekatef Date: Tue, 26 Dec 2023 23:46:26 +0300 Subject: [PATCH 16/50] Fix format --- scripts/build_climate_projections.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index bb4e6850d..349436f27 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -207,11 +207,11 @@ def calculate_proj_of_average(cmip6_xr, month, year0, year1, years_window): """ cmip6_interp_year0 = subset_by_time( - cmip6_xr, - month, - year=year0, + cmip6_xr, month, year=year0, years_window=years_window + ) + cmip6_interp_year1 = subset_by_time( + cmip6_xr, month=month, year=year1, years_window=years_window ) - cmip6_interp_year1 = subset_by_time(cmip6_xr, month=month, year=year1) dt_interp = cmip6_interp_year1["t"].mean("member").mean( "time" ) - cmip6_interp_year0["t"].mean("member").mean("time") From 0e507c129432a32b158fc48d171a8e35cf04e4cd Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 27 Dec 2023 22:49:23 +0300 Subject: [PATCH 17/50] Fix input --- Snakefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Snakefile b/Snakefile index 5c7d703a0..a1caeccaa 100644 --- a/Snakefile +++ b/Snakefile @@ -376,7 +376,7 @@ if config["enable"].get("modify_cutout", False): years_window=config["projection"]["years_window"], input: cutout="cutouts/" + CDIR + "{cutout}.nc", - cmip6_avr="/Users/ekaterina/Documents/_github_/cmip6/Global Atlas/tas_ssp245/t_CMIP6_ssp245_mon_201501-210012.nc", + cmip6_avr="data/tas_ssp245/t_CMIP6_ssp245_mon_201501-210012.nc", output: "cutouts/" + CDIR + "proj-{cutout}.nc", log: From 47c628aa69d8a065a6c91b01e4206d2c35c946e1 Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 27 Dec 2023 23:36:09 +0300 Subject: [PATCH 18/50] Remove hardcoding for a parameter name in a projection dataset --- scripts/build_climate_projections.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 349436f27..5983c657b 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -175,7 +175,9 @@ def subset_by_time(cmip6_xr, month, year, years_window): # TODO fix hardcoding in the parameter name # TODO add functionality to customise CMIP6 ensemble -def calculate_proj_of_average(cmip6_xr, month, year0, year1, years_window): +def calculate_proj_of_average( + cmip6_xr, month, year0, year1, years_window, cmip6_param_name="t" +): """ Calculate a change in a monthly average value for a climate parameter for the requested year periods and the processed month (a single value for each @@ -212,9 +214,9 @@ def calculate_proj_of_average(cmip6_xr, month, year0, year1, years_window): cmip6_interp_year1 = subset_by_time( cmip6_xr, month=month, year=year1, years_window=years_window ) - dt_interp = cmip6_interp_year1["t"].mean("member").mean( + dt_interp = cmip6_interp_year1[cmip6_param_name].mean("member").mean( "time" - ) - cmip6_interp_year0["t"].mean("member").mean("time") + ) - cmip6_interp_year0[cmip6_param_name].mean("member").mean("time") return dt_interp From dc735ba94210e612dc997d4ad6f60866e6433d45 Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 27 Dec 2023 23:36:58 +0300 Subject: [PATCH 19/50] Remove hardcoding for a parameter name in a cutout dataset --- scripts/build_climate_projections.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 5983c657b..fc23fb71d 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -221,7 +221,9 @@ def calculate_proj_of_average( # TODO fix hardcoding in the parameter name -def build_projection_for_month(cutout_xr, dt_xr, month): +def build_projection_for_month( + cutout_xr, dt_xr, month, cutout_param_name="temperature" +): """ Calculate the projection for the cutout data applying the morphing approach using a dataset dt_xr for the changes in monthly aggregates for the @@ -246,8 +248,8 @@ def build_projection_for_month(cutout_xr, dt_xr, month): for i in range(0, len(cutout_xr.y)): for j in range(0, len(cutout_xr.x)): - cutout_xr.temperature[k_time, i, j] = np.add( - cutout_xr.temperature[k_time, i, j], + cutout_xr[cutout_param_name][k_time, i, j] = np.add( + cutout_xr[cutout_param_name][k_time, i, j], np.array([dt_xr[i, j].item()] * k_time.sum().item()), ) From 4500f451f915671aeae6e4684dacf51b99f097c7 Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 27 Dec 2023 23:49:55 +0300 Subject: [PATCH 20/50] Wrap plots into functions --- scripts/build_climate_projections.py | 30 +++++++++++++++++++--------- 1 file changed, 21 insertions(+), 9 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index fc23fb71d..4144ebe7c 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -295,6 +295,23 @@ def build_cutout_future(cutout_xr, cmip6_xr, months, year0, year1, years_window) return cutout_xr +def plot_cmpi6(cmip6_xr, param_name="t"): + fig = plt.figure() + cmip6_xr[param_name].mean("time").mean("member").plot() + fig.savefig("results/test_cmip6_interp.png", dpi=700) + + +def plot_cutout( + cutout_xr, param_name="temperature", cmap="OrRd", fl_name="results/test_present.png" +): + fig = plt.figure() + if param_name == "temperature": + (cutout_xr[param_name] - 273.15).mean("time").plot(cmap=cmap) + else: + (cutout_xr[param_name]).mean("time").plot(cmap=cmap) + fig.savefig(fl_name, dpi=700) + + if __name__ == "__main__": if "snakemake" not in globals(): from _helpers import mock_snakemake @@ -332,13 +349,10 @@ def build_cutout_future(cutout_xr, cmip6_xr, months, year0, year1, years_window) # to be replaced after debug # ----------------------------------------------------------------- # graphical test of interpolation - fig = plt.figure() - cmip6_region_interp["t"].mean("time").mean("member").plot() - fig.savefig("test_cmip6_interp.png", dpi=700) + plot_cmpi6(cmip6_region_interp) + + plot_cutout(cutout_xr, fl_name="results/test_present.png") - fig = plt.figure() - (cutout_xr["temperature"] - 273.15).mean("time").plot(cmap="OrRd") - fig.savefig("results/test_present.png", dpi=700) # ----------------------------------------------------------------- # TODO Add a check for time match @@ -357,6 +371,4 @@ def build_cutout_future(cutout_xr, cmip6_xr, months, year0, year1, years_window) # to be replaced after debug # ----------------------------------------------------------------- # graphical test of projection - fig = plt.figure() - (cutout_future["temperature"] - 273.15).mean("time").plot(cmap="OrRd") - fig.savefig("results/test_cmip6_project.png", dpi=700) + plot_cutout(cutout_xr, fl_name="results/test_cmip6_project.png") From 7d8e17e9d9225452f5d883360900dfca59c66cb6 Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 27 Dec 2023 23:52:43 +0300 Subject: [PATCH 21/50] Update docstrings --- scripts/build_climate_projections.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 4144ebe7c..208621c8e 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -195,6 +195,8 @@ def calculate_proj_of_average( The first year of a later period in the future years_window: integer Width of the considered time period + cmip6_param_name: string + A name of CMIP6 parameter of interest Outputs ------- @@ -237,6 +239,8 @@ def build_projection_for_month( Dataset of the changes in the monthly averages for each grid cell month: float A month value to be considered further for the morphing procedure + cutout_param_name: string + A name of a cutout parameter for which a projection is being calculated Output ------- From 42da0b8bd952e17b91b303aa846c289df37f83af Mon Sep 17 00:00:00 2001 From: ekatef Date: Wed, 27 Dec 2023 23:56:16 +0300 Subject: [PATCH 22/50] Draft stretch implementation --- Snakefile | 1 + scripts/build_climate_projections.py | 5 +++++ 2 files changed, 6 insertions(+) diff --git a/Snakefile b/Snakefile index a1caeccaa..94a1bef88 100644 --- a/Snakefile +++ b/Snakefile @@ -374,6 +374,7 @@ if config["enable"].get("modify_cutout", False): present_year=config["projection"]["present_year"], future_year=config["projection"]["future_year"], years_window=config["projection"]["years_window"], + stretch_input=config["projection"]["stretch_input"], input: cutout="cutouts/" + CDIR + "{cutout}.nc", cmip6_avr="data/tas_ssp245/t_CMIP6_ssp245_mon_201501-210012.nc", diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 208621c8e..65cfdede7 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -335,6 +335,11 @@ def plot_cutout( present_year = snakemake.params.present_year future_year = snakemake.params.future_year years_window = snakemake.params.years_window + stretch_input = snakemake.params.stretch_input + + # if stretch_input: + # pass + cutout = snakemake.input.cutout cmip6 = snakemake.input.cmip6_avr From e4780a9cb9e0bb10d2f12df86e559d7d20147afa Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 08:59:48 +0300 Subject: [PATCH 23/50] Remove a TODO --- scripts/build_climate_projections.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 65cfdede7..bd88e52f3 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -173,7 +173,6 @@ def subset_by_time(cmip6_xr, month, year, years_window): return cmip6_in_period -# TODO fix hardcoding in the parameter name # TODO add functionality to customise CMIP6 ensemble def calculate_proj_of_average( cmip6_xr, month, year0, year1, years_window, cmip6_param_name="t" @@ -222,7 +221,6 @@ def calculate_proj_of_average( return dt_interp -# TODO fix hardcoding in the parameter name def build_projection_for_month( cutout_xr, dt_xr, month, cutout_param_name="temperature" ): From dbe81ac89aaebd894b7da77584873437108a4191 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 09:01:34 +0300 Subject: [PATCH 24/50] Use numpy functions --- scripts/build_climate_projections.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index bd88e52f3..6cc46c236 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -124,13 +124,13 @@ def interpolate_cmip6_to_cutout_grid(cmip6_xr, cutout_xr): dx_new = 0.3 newlon = np.arange( - round(min(cutout_xr.coords["x"].values), 1), - round(max(cutout_xr.coords["x"].values) + dx_new, 1), + round(np.min(cutout_xr.coords["x"].values), 1), + round(np.max(cutout_xr.coords["x"].values) + dx_new, 1), dx_new, ) newlat = np.arange( - round(min(cutout_xr.coords["y"].values), 1), - round(max(cutout_xr.coords["y"].values) + dx_new, 1), + round(np.min(cutout_xr.coords["y"].values), 1), + round(np.max(cutout_xr.coords["y"].values) + dx_new, 1), dx_new, ) cmip6_interp = cmip6_xr.interp(time=cmip6_xr.time, lat=newlat, lon=newlon) From 5fa9712fdd056c28ba9552612b2ea97fdacc6be3 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 11:10:43 +0300 Subject: [PATCH 25/50] Clarify a docstring --- scripts/build_climate_projections.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 6cc46c236..ee141d7b8 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -50,7 +50,8 @@ x_future and x_past have the highest temporal resolution, |month is the monthly average for a considered month. -The procedure is applied for each month to account for seasonality. +The procedure is applied for each month to account for seasonality. Snapshots settings +are used to define for which months to be used in calculating a climate projection. Methodology notes ----------------- From b4e1893f3bd14a0103abf57d041b5eb8a850c20f Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 12:14:26 +0300 Subject: [PATCH 26/50] Improve inputs --- Snakefile | 3 ++- config.default.yaml | 2 ++ config.tutorial.yaml | 2 ++ scripts/build_climate_projections.py | 6 ++---- 4 files changed, 8 insertions(+), 5 deletions(-) diff --git a/Snakefile b/Snakefile index 94a1bef88..f3397ba5e 100644 --- a/Snakefile +++ b/Snakefile @@ -374,7 +374,8 @@ if config["enable"].get("modify_cutout", False): present_year=config["projection"]["present_year"], future_year=config["projection"]["future_year"], years_window=config["projection"]["years_window"], - stretch_input=config["projection"]["stretch_input"], + param_nn_fl=config["projection"]["param_nn_fl"], + param_xx_fl=config["projection"]["param_xx_fl"], input: cutout="cutouts/" + CDIR + "{cutout}.nc", cmip6_avr="data/tas_ssp245/t_CMIP6_ssp245_mon_201501-210012.nc", diff --git a/config.default.yaml b/config.default.yaml index 8c608534f..e5ffb4bb8 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -209,6 +209,8 @@ projection: future_year: 2050 years_window: 5 # a half of width currently gcm_selection: false # false to use all the available global climate models + param_nn_fl: data/cmip6/tnn_CMIP6_ssp245_mon_201501-210012.nc # false or names of the files with trends for extrema values + param_xx_fl: data/cmip6/txx_CMIP6_ssp245_mon_201501-210012.nc renewable: diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 70e0732a9..9cd0f2530 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -224,6 +224,8 @@ projection: future_year: 2050 years_window: 5 # a half of width currently gcm_selection: false # false to use all the available global climate models + param_nn_fl: data/cmip6/tnn_CMIP6_ssp245_mon_201501-210012.nc # false or names of the files with trends for extrema values + param_xx_fl: data/cmip6/txx_CMIP6_ssp245_mon_201501-210012.nc renewable: diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index ee141d7b8..3802bbf7a 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -334,10 +334,8 @@ def plot_cutout( present_year = snakemake.params.present_year future_year = snakemake.params.future_year years_window = snakemake.params.years_window - stretch_input = snakemake.params.stretch_input - - # if stretch_input: - # pass + param_nn_fl = snakemake.params.param_nn_fl + param_xx_fl = snakemake.params.param_xx_fl cutout = snakemake.input.cutout cmip6 = snakemake.input.cmip6_avr From a292e146019b93ba9b3ecd9a4c3974111d495dfd Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 12:14:44 +0300 Subject: [PATCH 27/50] Fix names of the inputs --- Snakefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Snakefile b/Snakefile index f3397ba5e..f7fd6fe0d 100644 --- a/Snakefile +++ b/Snakefile @@ -378,7 +378,7 @@ if config["enable"].get("modify_cutout", False): param_xx_fl=config["projection"]["param_xx_fl"], input: cutout="cutouts/" + CDIR + "{cutout}.nc", - cmip6_avr="data/tas_ssp245/t_CMIP6_ssp245_mon_201501-210012.nc", + cmip6_avr="data/cmip6/t_CMIP6_ssp245_mon_201501-210012.nc", output: "cutouts/" + CDIR + "proj-{cutout}.nc", log: From 3715f9fc28d4f141ee52afc9098a8a56b5047301 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 14:28:57 +0300 Subject: [PATCH 28/50] Improve naming --- Snakefile | 4 ++-- config.default.yaml | 4 ++-- config.tutorial.yaml | 4 ++-- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/Snakefile b/Snakefile index f7fd6fe0d..a7a83cc6c 100644 --- a/Snakefile +++ b/Snakefile @@ -374,8 +374,8 @@ if config["enable"].get("modify_cutout", False): present_year=config["projection"]["present_year"], future_year=config["projection"]["future_year"], years_window=config["projection"]["years_window"], - param_nn_fl=config["projection"]["param_nn_fl"], - param_xx_fl=config["projection"]["param_xx_fl"], + cmip6_nn_fl=config["projection"]["cmip6_nn_fl"], + cmip6_xx_fl=config["projection"]["cmip6_xx_fl"], input: cutout="cutouts/" + CDIR + "{cutout}.nc", cmip6_avr="data/cmip6/t_CMIP6_ssp245_mon_201501-210012.nc", diff --git a/config.default.yaml b/config.default.yaml index e5ffb4bb8..621bae3ba 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -209,8 +209,8 @@ projection: future_year: 2050 years_window: 5 # a half of width currently gcm_selection: false # false to use all the available global climate models - param_nn_fl: data/cmip6/tnn_CMIP6_ssp245_mon_201501-210012.nc # false or names of the files with trends for extrema values - param_xx_fl: data/cmip6/txx_CMIP6_ssp245_mon_201501-210012.nc + cmip6_nn_fl: data/cmip6/tnn_CMIP6_ssp245_mon_201501-210012.nc # false or names of the files with trends for extrema values + cmip6_xx_fl: data/cmip6/txx_CMIP6_ssp245_mon_201501-210012.nc renewable: diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 9cd0f2530..f5d0c72a1 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -224,8 +224,8 @@ projection: future_year: 2050 years_window: 5 # a half of width currently gcm_selection: false # false to use all the available global climate models - param_nn_fl: data/cmip6/tnn_CMIP6_ssp245_mon_201501-210012.nc # false or names of the files with trends for extrema values - param_xx_fl: data/cmip6/txx_CMIP6_ssp245_mon_201501-210012.nc + cmip6_nn_fl: data/cmip6/tnn_CMIP6_ssp245_mon_201501-210012.nc # false or names of the files with trends for extrema values + cmip6_xx_fl: data/cmip6/txx_CMIP6_ssp245_mon_201501-210012.nc renewable: From a33b9ecc159d2d6cec49b8d308967bc7a5ea4ebf Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 14:29:35 +0300 Subject: [PATCH 29/50] Add a stretch-reference to the docstring --- scripts/build_climate_projections.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 3802bbf7a..e45a84185 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -60,7 +60,9 @@ J.B. Dias, G.C. da Graça, P.M.M. Soares (2020) Comparison of methodologies for generation of future building thermal energy simulation. Energy and Buildings, -vol. 206,2020, 109556, https://doi.org/10.1016/j.enbuild.2019.109556. +vol. 206, 109556, https://doi.org/10.1016/j.enbuild.2019.109556. + +J. Pulkkinen, J.-N. Louis (2021) Near- and medium-term hourly morphed mean and extreme future temperature datasets for Jyväskylä, Finland, for building thermal energy demand simulations. Data in Brief, vol. 37, 107209, https://doi.org/10.1016/j.dib.2021.107209. """ From d061ef316c4bfb6d51928eb9efe414e2b172c052 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 14:31:02 +0300 Subject: [PATCH 30/50] Wrap preparation of cmip6-xarrays into a function --- scripts/build_climate_projections.py | 39 +++++++++++++++++++++++----- 1 file changed, 32 insertions(+), 7 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index e45a84185..e043165e3 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -141,6 +141,25 @@ def interpolate_cmip6_to_cutout_grid(cmip6_xr, cutout_xr): return cmip6_interp +def prepare_cmip6(cmip6_fl, cutout_xr): + """ + Prepare the CMIP6 for calculation projections for the area of interest. + + Parameters + ---------- + cmip6_xr: xarray + CMIP6 climate projections loaded as xarray + cutout_xr : xarray of a cutout dataset + Cutout dataset loaded as xarray + """ + cmip6_xr = xr.open_dataset(cmip6_fl) + cmip6_cropped_xr = crop_cmip6(cmip6_xr, cutout_xr) + cmip6_interp = interpolate_cmip6_to_cutout_grid( + cmip6_xr=cmip6_cropped_xr, cutout_xr=cutout_xr + ) + return cmip6_interp + + # TODO fix years_window def subset_by_time(cmip6_xr, month, year, years_window): """ @@ -336,8 +355,8 @@ def plot_cutout( present_year = snakemake.params.present_year future_year = snakemake.params.future_year years_window = snakemake.params.years_window - param_nn_fl = snakemake.params.param_nn_fl - param_xx_fl = snakemake.params.param_xx_fl + cmip6_nn_fl = snakemake.params.cmip6_nn_fl + cmip6_xx_fl = snakemake.params.cmip6_xx_fl cutout = snakemake.input.cutout cmip6 = snakemake.input.cmip6_avr @@ -345,13 +364,19 @@ def plot_cutout( snapshots = pd.date_range(freq="h", **snakemake.params.snapshots) season_in_focus = snapshots.month.unique().to_list() - cmip6_xr = xr.open_dataset(cmip6) cutout_xr = xr.open_dataset(cutout) - cmip6_region = crop_cmip6(cmip6_xr, cutout_xr) - cmip6_region_interp = interpolate_cmip6_to_cutout_grid( - cmip6_xr=cmip6_region, cutout_xr=cutout_xr - ) + cmip6_region_interp = prepare_cmip6(cmip6_fl=cmip6, cutout_xr=cutout_xr) + + if param_nn_fl: + cmip6_nn_region_interp = prepare_cmip6( + cmip6_fl=param_nn_fl, cutout_xr=cutout_xr + ) + + if param_xx_fl: + cmip6_xx_region_interp = prepare_cmip6( + cmip6_fl=param_xx_fl, cutout_xr=cutout_xr + ) # ----------------------------------------------------------------- # to be replaced after debug From 1c46763a044e66efe4758f01a5cbaf8610b8886f Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 22:41:22 +0300 Subject: [PATCH 31/50] Add stretch --- scripts/build_climate_projections.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index e043165e3..d15b3b787 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -244,7 +244,7 @@ def calculate_proj_of_average( def build_projection_for_month( - cutout_xr, dt_xr, month, cutout_param_name="temperature" + cutout_xr, dt_xr, month, a, t_mean_month, cutout_param_name="temperature" ): """ Calculate the projection for the cutout data applying the morphing approach @@ -270,6 +270,9 @@ def build_projection_for_month( k_time = cutout_xr.time.dt.month.isin(month) + if a & t_mean_month: + dt_xr = dt_xr + a * (cutout_xr[cutout_param_name][k_time, :, :] - t_mean_month) + for i in range(0, len(cutout_xr.y)): for j in range(0, len(cutout_xr.x)): cutout_xr[cutout_param_name][k_time, i, j] = np.add( From f3825f25a95919b0f7733e5c698492b62cbe05ab Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 22:42:46 +0300 Subject: [PATCH 32/50] Implement stretch --- scripts/build_climate_projections.py | 47 ++++++++++++++++++++++++++-- 1 file changed, 45 insertions(+), 2 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index d15b3b787..539c1130f 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -283,7 +283,18 @@ def build_projection_for_month( return cutout_xr -def build_cutout_future(cutout_xr, cmip6_xr, months, year0, year1, years_window): +def build_cutout_future( + cutout_xr, + cmip6_xr, + months, + year0, + year1, + years_window, + cutout_param_name="temperature", + add_stretch=False, + cmip6_nn_xr=False, + cmip6_xx_xr=False, +): """ Calculate the projection for the cutout data applying the morphing approach for each month in the months. @@ -315,7 +326,39 @@ def build_cutout_future(cutout_xr, cmip6_xr, months, year0, year1, years_window) cmip6_xr=cmip6_xr, month=k_month, year0=year0, year1=year1, years_window=5 ) - build_projection_for_month(cutout_xr, dt_k, k_month) + a = False + t_mean_month = False + + if add_stretch: + # TODO Check cmip6_nn_xr and cmip6_xx_xr are not empty + dt_nn_k = calculate_proj_of_average( + cmip6_xr=cmip6_nn_xr, + month=k_month, + year0=year0, + year1=year1, + years_window=5, + ) + dt_xx_k = calculate_proj_of_average( + cmip6_xr=cmip6_xx_xr, + month=k_month, + year0=year0, + year1=year1, + years_window=5, + ) + ddt_nx = dt_xx_k - dt_nn_k + k_time = cutout_xr.time.dt.month.isin(month) + t_mean_month = cutout_xr[cutout_param_name][k_time, :, :].mean() + + a = ddt_nx / t_mean_month + + build_projection_for_month( + cutout_xr, + dt_k, + k_month, + a=a, + t_mean_month=t_mean_month, + cutout_param_name=cutout_param_name, + ) cutout_xr = cutout_xr.where(cutout_xr.time.dt.month.isin(months), drop=True) From 989ae50a2083047ecc0cdb1072ec3a0872edd965 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 22:43:00 +0300 Subject: [PATCH 33/50] Fix names --- scripts/build_climate_projections.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 539c1130f..8a372e4fc 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -414,14 +414,14 @@ def plot_cutout( cmip6_region_interp = prepare_cmip6(cmip6_fl=cmip6, cutout_xr=cutout_xr) - if param_nn_fl: + if cmip6_nn_fl: cmip6_nn_region_interp = prepare_cmip6( - cmip6_fl=param_nn_fl, cutout_xr=cutout_xr + cmip6_fl=cmip6_nn_fl, cutout_xr=cutout_xr ) - if param_xx_fl: + if cmip6_xx_fl: cmip6_xx_region_interp = prepare_cmip6( - cmip6_fl=param_xx_fl, cutout_xr=cutout_xr + cmip6_fl=cmip6_xx_fl, cutout_xr=cutout_xr ) # ----------------------------------------------------------------- From f68bf21cb096cffd407ee26c8175a6e2dbcbd6f7 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 22:43:47 +0300 Subject: [PATCH 34/50] Test stretch --- scripts/build_climate_projections.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 8a372e4fc..34c43f9ff 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -450,4 +450,22 @@ def plot_cutout( # to be replaced after debug # ----------------------------------------------------------------- # graphical test of projection - plot_cutout(cutout_xr, fl_name="results/test_cmip6_project.png") + plot_cutout(cutout_future, fl_name="results/test_cmip6_project_shift.png") + + cutout_future2 = build_cutout_future( + cutout_xr=cutout_xr, + cmip6_xr=cmip6_region_interp, + months=season_in_focus, + year0=present_year, + year1=future_year, + years_window=years_window, + cutout_param_name="temperature", + cmip6_nn_xr=cmip6_nn_fl, + cmip6_xx_xr=cmip6_xx_fl, + ) + + # ----------------------------------------------------------------- + # to be replaced after debug + # ----------------------------------------------------------------- + # graphical test of projection + plot_cutout(cutout_future2, fl_name="results/test_cmip6_project_full.png") From ba3d1a76106d41b82bddbe1a54abbe4a81b66f00 Mon Sep 17 00:00:00 2001 From: ekatef Date: Thu, 28 Dec 2023 22:44:32 +0300 Subject: [PATCH 35/50] Fix plot parameters --- scripts/build_climate_projections.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 34c43f9ff..67421f94a 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -365,10 +365,10 @@ def build_cutout_future( return cutout_xr -def plot_cmpi6(cmip6_xr, param_name="t"): +def plot_cmpi6(cmip6_xr, param_name="t", fl_name="results/test_cmip6.png"): fig = plt.figure() cmip6_xr[param_name].mean("time").mean("member").plot() - fig.savefig("results/test_cmip6_interp.png", dpi=700) + fig.savefig(fl_name, dpi=700) def plot_cutout( @@ -428,7 +428,13 @@ def plot_cutout( # to be replaced after debug # ----------------------------------------------------------------- # graphical test of interpolation - plot_cmpi6(cmip6_region_interp) + plot_cmpi6(cmip6_region_interp, fl_name="results/test_cmip6.png") + plot_cmpi6( + cmip6_nn_region_interp, param_name="tnn", fl_name="results/test_cmip6_nn.png" + ) + plot_cmpi6( + cmip6_xx_region_interp, param_name="txx", fl_name="results/test_cmip6_xx.png" + ) plot_cutout(cutout_xr, fl_name="results/test_present.png") From c1bfe0c79e9643bbc251a4a538e15d5443453120 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 20 Jan 2024 21:39:59 +0300 Subject: [PATCH 36/50] Improve comments to the config parameters --- config.default.yaml | 4 ++-- config.tutorial.yaml | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index 621bae3ba..307f74cf5 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -209,8 +209,8 @@ projection: future_year: 2050 years_window: 5 # a half of width currently gcm_selection: false # false to use all the available global climate models - cmip6_nn_fl: data/cmip6/tnn_CMIP6_ssp245_mon_201501-210012.nc # false or names of the files with trends for extrema values - cmip6_xx_fl: data/cmip6/txx_CMIP6_ssp245_mon_201501-210012.nc + cmip6_nn_fl: data/cmip6/tnn_CMIP6_ssp245_mon_201501-210012.nc # false or names of min-extremal files of climate inputs + cmip6_xx_fl: data/cmip6/txx_CMIP6_ssp245_mon_201501-210012.nc # false or names of max-extremal files of climate inputs renewable: diff --git a/config.tutorial.yaml b/config.tutorial.yaml index f5d0c72a1..837aa6c8c 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -224,8 +224,8 @@ projection: future_year: 2050 years_window: 5 # a half of width currently gcm_selection: false # false to use all the available global climate models - cmip6_nn_fl: data/cmip6/tnn_CMIP6_ssp245_mon_201501-210012.nc # false or names of the files with trends for extrema values - cmip6_xx_fl: data/cmip6/txx_CMIP6_ssp245_mon_201501-210012.nc + cmip6_nn_fl: data/cmip6/tnn_CMIP6_ssp245_mon_201501-210012.nc # false or names of min-extremal files of climate inputs + cmip6_xx_fl: data/cmip6/txx_CMIP6_ssp245_mon_201501-210012.nc # false or names of max-extremal files of climate inputs renewable: From 9d21fdeb4abb588a002a6ce904f10e32e6a928e5 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 20 Jan 2024 22:20:02 +0300 Subject: [PATCH 37/50] Add comments to the projection calculations --- scripts/build_climate_projections.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 67421f94a..143a3652a 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -270,6 +270,7 @@ def build_projection_for_month( k_time = cutout_xr.time.dt.month.isin(month) + # stretch calculations if a & t_mean_month: dt_xr = dt_xr + a * (cutout_xr[cutout_param_name][k_time, :, :] - t_mean_month) @@ -277,6 +278,8 @@ def build_projection_for_month( for j in range(0, len(cutout_xr.x)): cutout_xr[cutout_param_name][k_time, i, j] = np.add( cutout_xr[cutout_param_name][k_time, i, j], + # spreading dt_xr which is of monthly resolution across the whole time dimension of a cutout + # which is normally of hourly resolution np.array([dt_xr[i, j].item()] * k_time.sum().item()), ) From a8f8636e9f467a609121b6f993d2f27218b4c4ab Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 9 Mar 2024 01:09:10 +0300 Subject: [PATCH 38/50] Fix formatting --- scripts/build_climate_projections.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 143a3652a..0d38ef0d4 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -62,7 +62,10 @@ generation of future building thermal energy simulation. Energy and Buildings, vol. 206, 109556, https://doi.org/10.1016/j.enbuild.2019.109556. -J. Pulkkinen, J.-N. Louis (2021) Near- and medium-term hourly morphed mean and extreme future temperature datasets for Jyväskylä, Finland, for building thermal energy demand simulations. Data in Brief, vol. 37, 107209, https://doi.org/10.1016/j.dib.2021.107209. +J. Pulkkinen, J.-N. Louis (2021) Near- and medium-term hourly morphed mean and +extreme future temperature datasets for Jyväskylä, Finland, for building thermal +energy demand simulations. Data in Brief, vol. 37, 107209, +https://doi.org/10.1016/j.dib.2021.107209. """ From 33d90d236a38bc1a8e7902c6e382ba8a1a4f64d0 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 9 Mar 2024 01:21:04 +0300 Subject: [PATCH 39/50] Enhance docstrings --- scripts/build_climate_projections.py | 24 ++++++++++++++++++++++-- 1 file changed, 22 insertions(+), 2 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 0d38ef0d4..f0459ceb9 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -96,6 +96,11 @@ def crop_cmip6(cmip6_xr, cutout_xr): CMIP6 climate projections loaded as xarray cutout_xr : xarray of a cutout dataset Cutout dataset loaded as xarray + + Returns + ------- + cmip6_region: xarray + Clipped according to an area of interest """ # spartial margin needed to avoid data losses during further interpolation @@ -124,6 +129,11 @@ def interpolate_cmip6_to_cutout_grid(cmip6_xr, cutout_xr): CMIP6 climate projections loaded as xarray cutout_xr : xarray of a cutout dataset Cutout dataset loaded as xarray + + Returns + ------- + cmip6_interp: xarray + Interpolated to the grid identical to the cutout one """ # TODO read from the cutout file instead of hardcoding @@ -154,6 +164,11 @@ def prepare_cmip6(cmip6_fl, cutout_xr): CMIP6 climate projections loaded as xarray cutout_xr : xarray of a cutout dataset Cutout dataset loaded as xarray + + Returns + ------- + cmip6_interp: xarray + Now clipped to the region and interpolated to match with the cutour grid """ cmip6_xr = xr.open_dataset(cmip6_fl) cmip6_cropped_xr = crop_cmip6(cmip6_xr, cutout_xr) @@ -180,6 +195,11 @@ def subset_by_time(cmip6_xr, month, year, years_window): years_window: integer A width of the considered time period + Returns + ------- + cmip6_in_period: xarray + Sub-setted for a month and a years range of interest + Example ------- month=3, year=2050, years_window=30 @@ -225,7 +245,7 @@ def calculate_proj_of_average( Outputs ------- dt_interp: xarray of the changes in the monthly averages for each grid cell - Contains a single value for each grid cell. + Contains a single value for each grid cell. Example ------- @@ -268,7 +288,7 @@ def build_projection_for_month( Output ------- cutout_xr: xarray of the projected changes for each grid cell for a processed month - Contains a time series for each grid cell. + Contains a time series for each grid cell. """ k_time = cutout_xr.time.dt.month.isin(month) From f5cb70e0dd8579b470aaef7aa987c430107df5f7 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 9 Mar 2024 01:23:02 +0300 Subject: [PATCH 40/50] Improve naming --- scripts/build_climate_projections.py | 50 +++++++++++++++------------- 1 file changed, 27 insertions(+), 23 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index f0459ceb9..9de678d11 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -220,7 +220,7 @@ def subset_by_time(cmip6_xr, month, year, years_window): # TODO add functionality to customise CMIP6 ensemble def calculate_proj_of_average( - cmip6_xr, month, year0, year1, years_window, cmip6_param_name="t" + cmip6_xr, month, base_year, predict_year, years_window, cmip6_param_name="t" ): """ Calculate a change in a monthly average value for a climate parameter for @@ -233,9 +233,9 @@ def calculate_proj_of_average( CMIP6 climate projections loaded as xarray month: float A month value to be considered further for the morphing procedure - year0: integer + base_year: integer The first year of an earlier period in the future - year1: integer + predict_year: integer The first year of a later period in the future years_window: integer Width of the considered time period @@ -249,20 +249,20 @@ def calculate_proj_of_average( Example ------- - month=3, year0=2020, year0=2070, years_window=30 + month=3, base_year=2020, base_year=2070, years_window=30 calculates a multimodel change in the monthly average for March in 2020-2049 as compared with 2070-2099 """ - cmip6_interp_year0 = subset_by_time( - cmip6_xr, month, year=year0, years_window=years_window + cmip6_interp_base_year = subset_by_time( + cmip6_xr, month, year=base_year, years_window=years_window ) - cmip6_interp_year1 = subset_by_time( - cmip6_xr, month=month, year=year1, years_window=years_window + cmip6_interp_predict_year = subset_by_time( + cmip6_xr, month=month, year=predict_year, years_window=years_window ) - dt_interp = cmip6_interp_year1[cmip6_param_name].mean("member").mean( + dt_interp = cmip6_interp_predict_year[cmip6_param_name].mean("member").mean( "time" - ) - cmip6_interp_year0[cmip6_param_name].mean("member").mean("time") + ) - cmip6_interp_base_year[cmip6_param_name].mean("member").mean("time") return dt_interp @@ -313,8 +313,8 @@ def build_cutout_future( cutout_xr, cmip6_xr, months, - year0, - year1, + base_year, + predict_year, years_window, cutout_param_name="temperature", add_stretch=False, @@ -334,9 +334,9 @@ def build_cutout_future( months: list Months values to be used as a definition of a season to be considered for building projection time series - year0: integer + base_year: integer The first year of an earlier period in the future - year1: integer + predict_year: integer The first year of a later period in the future years_window: integer Width of the considered time period @@ -349,7 +349,11 @@ def build_cutout_future( """ for k_month in [months]: dt_k = calculate_proj_of_average( - cmip6_xr=cmip6_xr, month=k_month, year0=year0, year1=year1, years_window=5 + cmip6_xr=cmip6_xr, + month=k_month, + base_year=base_year, + predict_year=predict_year, + years_window=5, ) a = False @@ -360,15 +364,15 @@ def build_cutout_future( dt_nn_k = calculate_proj_of_average( cmip6_xr=cmip6_nn_xr, month=k_month, - year0=year0, - year1=year1, + base_year=base_year, + predict_year=predict_year, years_window=5, ) dt_xx_k = calculate_proj_of_average( cmip6_xr=cmip6_xx_xr, month=k_month, - year0=year0, - year1=year1, + base_year=base_year, + predict_year=predict_year, years_window=5, ) ddt_nx = dt_xx_k - dt_nn_k @@ -471,8 +475,8 @@ def plot_cutout( cutout_xr=cutout_xr, cmip6_xr=cmip6_region_interp, months=season_in_focus, - year0=present_year, - year1=future_year, + base_year=present_year, + predict_year=future_year, years_window=years_window, ) @@ -488,8 +492,8 @@ def plot_cutout( cutout_xr=cutout_xr, cmip6_xr=cmip6_region_interp, months=season_in_focus, - year0=present_year, - year1=future_year, + base_year=present_year, + predict_year=future_year, years_window=years_window, cutout_param_name="temperature", cmip6_nn_xr=cmip6_nn_fl, From 2177bb19d9ebd92b70470c2b16fc53e7caf859af Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 9 Mar 2024 01:27:33 +0300 Subject: [PATCH 41/50] Fix parameters in a function call --- scripts/build_climate_projections.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 9de678d11..35f0d1026 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -353,7 +353,7 @@ def build_cutout_future( month=k_month, base_year=base_year, predict_year=predict_year, - years_window=5, + years_window=years_window, ) a = False @@ -366,14 +366,14 @@ def build_cutout_future( month=k_month, base_year=base_year, predict_year=predict_year, - years_window=5, + years_window=years_window, ) dt_xx_k = calculate_proj_of_average( cmip6_xr=cmip6_xx_xr, month=k_month, base_year=base_year, predict_year=predict_year, - years_window=5, + years_window=years_window, ) ddt_nx = dt_xx_k - dt_nn_k k_time = cutout_xr.time.dt.month.isin(month) From 19b65fc5a4bb899b8e7817743b9d5e7a6b52cf04 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 9 Mar 2024 01:40:25 +0300 Subject: [PATCH 42/50] Improve explanation for stretch calculations --- scripts/build_climate_projections.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 35f0d1026..35d9c8cea 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -50,6 +50,9 @@ x_future and x_past have the highest temporal resolution, |month is the monthly average for a considered month. +A value of the stretch is calculated using projections for minimal and maximal values +of the climate parameter. + The procedure is applied for each month to account for seasonality. Snapshots settings are used to define for which months to be used in calculating a climate projection. @@ -340,6 +343,14 @@ def build_cutout_future( The first year of a later period in the future years_window: integer Width of the considered time period + cutout_param_name: string + A name of the variable to be read from the cutout + add_stretch: boolean + Add a stretch transformation to a shift when calculating a projection + cmip6_nn_xr: xarray + CMIP6 projection for the minimum of the climate parameter + cmip6_xx_xr: xarray + CMIP6 projection for the maximum of the climate parameter Output ------- From 7ca6216d4323de61d09da15ba69b728ab4915eba Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 9 Mar 2024 18:42:56 +0300 Subject: [PATCH 43/50] Improve naming and add enable parameters to the configs --- Snakefile | 2 +- config.default.yaml | 1 + config.tutorial.yaml | 1 + 3 files changed, 3 insertions(+), 1 deletion(-) diff --git a/Snakefile b/Snakefile index 1bc09ab6b..9937ff1bb 100644 --- a/Snakefile +++ b/Snakefile @@ -365,7 +365,7 @@ if config["enable"].get("build_cutout", False): "scripts/build_cutout.py" -if config["enable"].get("modify_cutout", False): +if config["enable"].get("climate_projection_cutout", False): rule build_climate_projections: params: diff --git a/config.default.yaml b/config.default.yaml index 8e876b872..493bf2e71 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -21,6 +21,7 @@ enable: # If "build_cutout" : true, then environmental data is extracted according to `snapshots` date range and `countries` # requires cds API key https://cds.climate.copernicus.eu/api-how-to # More information https://atlite.readthedocs.io/en/latest/introduction.html#datasets + climate_projection_cutout: false custom_rules: [] # Default empty [] or link to custom rule file e.g. ["my_folder/my_rules.smk"] that add rules to Snakefile diff --git a/config.tutorial.yaml b/config.tutorial.yaml index c40b1645b..475e5080d 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -30,6 +30,7 @@ enable: # If "build_cutout" : true # requires cds API key https://cds.climate.copernicus.eu/api-how-to # More information https://atlite.readthedocs.io/en/latest/introduction.html#datasets build_cutout: false + climate_projection_cutout: false build_natura_raster: true # If True, then build_natura_raster can be run custom_rules: [] # Default empty [] or link to custom rule file e.g. ["my_folder/my_rules.smk"] that add rules to Snakefile From 39fd2d1d193ced7497883941c09b538d7b6f6d1f Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 9 Mar 2024 18:43:17 +0300 Subject: [PATCH 44/50] Fix typo --- scripts/build_climate_projections.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 35d9c8cea..06bea9cd2 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -387,7 +387,7 @@ def build_cutout_future( years_window=years_window, ) ddt_nx = dt_xx_k - dt_nn_k - k_time = cutout_xr.time.dt.month.isin(month) + k_time = cutout_xr.time.dt.month.isin(k_month) t_mean_month = cutout_xr[cutout_param_name][k_time, :, :].mean() a = ddt_nx / t_mean_month From 517e4228d1dcf25d068bdd956a7460554a5dfa46 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 9 Mar 2024 18:43:43 +0300 Subject: [PATCH 45/50] Replace hard-coding --- scripts/build_climate_projections.py | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 06bea9cd2..69052a18f 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -138,19 +138,15 @@ def interpolate_cmip6_to_cutout_grid(cmip6_xr, cutout_xr): cmip6_interp: xarray Interpolated to the grid identical to the cutout one """ - - # TODO read from the cutout file instead of hardcoding - dx_new = 0.3 - newlon = np.arange( round(np.min(cutout_xr.coords["x"].values), 1), - round(np.max(cutout_xr.coords["x"].values) + dx_new, 1), - dx_new, + round(np.max(cutout_xr.coords["x"].values) + cutout_xr.attrs["dx"], 1), + cutout_xr.attrs["dx"], ) newlat = np.arange( round(np.min(cutout_xr.coords["y"].values), 1), - round(np.max(cutout_xr.coords["y"].values) + dx_new, 1), - dx_new, + round(np.max(cutout_xr.coords["y"].values) + cutout_xr.attrs["dy"], 1), + cutout_xr.attrs["dy"], ) cmip6_interp = cmip6_xr.interp(time=cmip6_xr.time, lat=newlat, lon=newlon) From 1f3520d25ac7b96e860d7c01939eee7ad40350c3 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sat, 9 Mar 2024 19:09:40 +0300 Subject: [PATCH 46/50] Improve naming --- config.default.yaml | 2 +- config.tutorial.yaml | 2 +- scripts/build_climate_projections.py | 12 ++++++------ 3 files changed, 8 insertions(+), 8 deletions(-) diff --git a/config.default.yaml b/config.default.yaml index 493bf2e71..abcc5e717 100644 --- a/config.default.yaml +++ b/config.default.yaml @@ -206,7 +206,7 @@ atlite: projection: climate_scenario: ssp245 - present_year: 2020 + base_year: 2020 future_year: 2050 years_window: 5 # a half of width currently gcm_selection: false # false to use all the available global climate models diff --git a/config.tutorial.yaml b/config.tutorial.yaml index 475e5080d..8603a2780 100644 --- a/config.tutorial.yaml +++ b/config.tutorial.yaml @@ -221,7 +221,7 @@ atlite: projection: climate_scenario: ssp245 - present_year: 2020 + base_year: 2020 future_year: 2050 years_window: 5 # a half of width currently gcm_selection: false # false to use all the available global climate models diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 69052a18f..604e49af0 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -14,7 +14,7 @@ projection: climate_scenario: - present_year: + base_year: future_year: years_window: gcm_selection: @@ -41,7 +41,7 @@ an ensemble of CMIP6 global climate models to account for the future climate. A morphing approach is implemented meaning that the projection value x_future is calculated -data for a representative period present_year in the past x_past using a combination +data for a representative period base_year in the past x_past using a combination of shift dx and stretch a: x_future = x_past + dx + a(x_past - |month), @@ -427,7 +427,7 @@ def plot_cutout( snakemake = mock_snakemake( "build_climate_projections", climate_scenario="ssp2-2.6", - present_year=2020, + base_year=2020, future_year=2050, years_window=5, cutout="africa-2013-era5", @@ -435,7 +435,7 @@ def plot_cutout( configure_logging(snakemake) climate_scenario = snakemake.params.climate_scenario - present_year = snakemake.params.present_year + base_year = snakemake.params.base_year future_year = snakemake.params.future_year years_window = snakemake.params.years_window cmip6_nn_fl = snakemake.params.cmip6_nn_fl @@ -482,7 +482,7 @@ def plot_cutout( cutout_xr=cutout_xr, cmip6_xr=cmip6_region_interp, months=season_in_focus, - base_year=present_year, + base_year=base_year, predict_year=future_year, years_window=years_window, ) @@ -499,7 +499,7 @@ def plot_cutout( cutout_xr=cutout_xr, cmip6_xr=cmip6_region_interp, months=season_in_focus, - base_year=present_year, + base_year=base_year, predict_year=future_year, years_window=years_window, cutout_param_name="temperature", From dd4b0600abbdf7ba3f268945e93a2fccd5c68bb6 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 10 Mar 2024 00:12:23 +0300 Subject: [PATCH 47/50] Add a check for time match --- scripts/build_climate_projections.py | 52 +++++++++++++++++++++++++++- 1 file changed, 51 insertions(+), 1 deletion(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 604e49af0..2fb53bbf2 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -170,6 +170,49 @@ def prepare_cmip6(cmip6_fl, cutout_xr): Now clipped to the region and interpolated to match with the cutour grid """ cmip6_xr = xr.open_dataset(cmip6_fl) + + # All the requested data from the projection period should present in the CMIP6 datasets + year_month_df = pd.DataFrame( + { + "year": cmip6_xr.time.dt.year.to_series(), + "month": cmip6_xr.time.dt.month.to_series(), + } + ) + cmip6_in_base = year_month_df.apply( + lambda row: (row["year"] >= (base_year)) + & (row["year"] < (base_year + years_window)) + & (row["month"] in months), + axis=1, + ) + cmip6_in_future = year_month_df.apply( + lambda row: (row["year"] <= (predict_year)) + & (row["year"] > (predict_year - years_window)) + & (row["month"] in months), + axis=1, + ) + + if sum(cmip6_in_base) == 0: + raise Exception( + f"No data to build a projection are available in the provided CMIP6 data " + f"for the base year {base_year}-{base_year + years_window - 1}.\n\r" + ) + if sum(cmip6_in_future) == 0: + raise Exception( + f"No data to build a projection are available in the provided CMIP6 data " + f"for the base year {predict_year}-{predict_year - years_window + 1}.\n\r" + ) + + if sum(cmip6_in_base) < (len(months) * years_window): + logger.warning( + f"The requested projection period is not fully covered by the provided CMIP6 data. " + f"The requested base period is {base_year}-{base_year + years_window - 1}\n\r" + ) + if sum(cmip6_in_future) < (len(months) * years_window): + logger.warning( + f"The requested projection period is not fully covered by the provided CMIP6 data. " + f"The requested base period is {predict_year}-{predict_year - years_window + 1}\n\r" + ) + cmip6_cropped_xr = crop_cmip6(cmip6_xr, cutout_xr) cmip6_interp = interpolate_cmip6_to_cutout_grid( cmip6_xr=cmip6_cropped_xr, cutout_xr=cutout_xr @@ -449,7 +492,14 @@ def plot_cutout( cutout_xr = xr.open_dataset(cutout) - cmip6_region_interp = prepare_cmip6(cmip6_fl=cmip6, cutout_xr=cutout_xr) + cmip6_region_interp = prepare_cmip6( + cmip6_fl=cmip6, + cutout_xr=cutout_xr, + months=season_in_focus, + base_year=base_year, + predict_year=future_year, + years_window=years_window, + ) if cmip6_nn_fl: cmip6_nn_region_interp = prepare_cmip6( From f71390f994d0c2ef526dab5506c29f666bbc7ec6 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 10 Mar 2024 00:12:35 +0300 Subject: [PATCH 48/50] Fix typo --- Snakefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Snakefile b/Snakefile index 9937ff1bb..975029b75 100644 --- a/Snakefile +++ b/Snakefile @@ -371,7 +371,7 @@ if config["enable"].get("climate_projection_cutout", False): params: snapshots=config["snapshots"], climate_scenario=config["projection"]["climate_scenario"], - present_year=config["projection"]["present_year"], + base_year=config["projection"]["base_year"], future_year=config["projection"]["future_year"], years_window=config["projection"]["years_window"], cmip6_nn_fl=config["projection"]["cmip6_nn_fl"], From 14707c8e02c040f1c89a82451596847da0718b28 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 10 Mar 2024 00:13:10 +0300 Subject: [PATCH 49/50] Update the docstring --- scripts/build_climate_projections.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index 2fb53bbf2..c8735e03a 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -153,7 +153,14 @@ def interpolate_cmip6_to_cutout_grid(cmip6_xr, cutout_xr): return cmip6_interp -def prepare_cmip6(cmip6_fl, cutout_xr): +def prepare_cmip6( + cmip6_fl, + cutout_xr, + months, + base_year, + predict_year, + years_window, +): """ Prepare the CMIP6 for calculation projections for the area of interest. @@ -163,6 +170,14 @@ def prepare_cmip6(cmip6_fl, cutout_xr): CMIP6 climate projections loaded as xarray cutout_xr : xarray of a cutout dataset Cutout dataset loaded as xarray + months: list + A list of month numbers to be considered further projection + base_year: integer + The first year of an earlier period in the future + predict_year: integer + The first year of a later period in the future + years_window: integer + Width of the considered time period Returns ------- From 47e7946e6b3729369d9af7d57961d02db8b56921 Mon Sep 17 00:00:00 2001 From: ekatef Date: Sun, 10 Mar 2024 00:13:21 +0300 Subject: [PATCH 50/50] Remove TODOs --- scripts/build_climate_projections.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/scripts/build_climate_projections.py b/scripts/build_climate_projections.py index c8735e03a..5f58430c8 100644 --- a/scripts/build_climate_projections.py +++ b/scripts/build_climate_projections.py @@ -235,7 +235,6 @@ def prepare_cmip6( return cmip6_interp -# TODO fix years_window def subset_by_time(cmip6_xr, month, year, years_window): """ Filter CMIP6 dataset according to the requested year and the processed @@ -263,7 +262,6 @@ def subset_by_time(cmip6_xr, month, year, years_window): filters CMIP6 cmip6_xr dataset to calculate projection for March at 2050-2079 """ - # TODO add warning in case there are no enough years cmip6_in_period = cmip6_xr.where( ( (cmip6_xr["time.month"] == month) @@ -542,7 +540,6 @@ def plot_cutout( # ----------------------------------------------------------------- - # TODO Add a check for time match cutout_future = build_cutout_future( cutout_xr=cutout_xr, cmip6_xr=cmip6_region_interp,