CMIP6 Data generation & merging#

CMIP6 data will be generated from 1995-2014, loaded and merged into one file

Note: Depending on if you want to run historical, or an ssp, changes you need to make are:

  1. Cell 5: Select scenario

  2. Cell 6: Define period

Start-up#

# general python
import warnings
warnings.filterwarnings("ignore", category=UserWarning)

import numpy as np
from pathlib import Path
import pandas as pd
import matplotlib.pyplot as plt
import xarray as xr

#niceties
from rich import print

import yaml
# General eWaterCycle
import ewatercycle
import ewatercycle.models
import ewatercycle.forcing
---------------------------------------------------------------------------
OSError                                   Traceback (most recent call last)
Cell In[1], line 2
      1 # General eWaterCycle
----> 2 import ewatercycle
      3 import ewatercycle.models
      4 import ewatercycle.forcing

File /opt/conda/envs/ewatercycle2/lib/python3.12/site-packages/ewatercycle/__init__.py:3
      1 """ewatercycle package."""
----> 3 from .config import CFG as CFG
      4 from .version import __version__ as __version__

File /opt/conda/envs/ewatercycle2/lib/python3.12/site-packages/ewatercycle/config/__init__.py:321
    317 _SOURCES = (USER_HOME_CONFIG, SYSTEM_CONFIG)
    319 USER_CONFIG = _find_user_config(_SOURCES)
--> 321 CFG = Configuration._load_user_config(USER_CONFIG) if USER_CONFIG else Configuration()
    322 """eWaterCycle configuration object.
    323 
    324 The configuration is loaded from:
   (...)
    341     # apptainer pull docker://ewatercycle/wflow-grpc4bmi:2020.1.1
    342 """

File /opt/conda/envs/ewatercycle2/lib/python3.12/site-packages/ewatercycle/config/__init__.py:197, in Configuration._load_user_config(cls, filename)
    195 mapping = _read_config_file(filename)
    196 try:
--> 197     return Configuration(ewatercycle_config=filename, **mapping)
    198 except ValidationError as e:
    199     # Append filename to error locs
    200     for error in e.errors():

    [... skipping hidden 1 frame]

File /opt/conda/envs/ewatercycle2/lib/python3.12/site-packages/pydantic/types.py:1282, in PathType.validate_directory(path, _)
   1280 @staticmethod
   1281 def validate_directory(path: Path, _: core_schema.ValidationInfo) -> Path:
-> 1282     if path.is_dir():
   1283         return path
   1284     else:

File /opt/conda/envs/ewatercycle2/lib/python3.12/pathlib.py:875, in Path.is_dir(self)
    871 """
    872 Whether this path is a directory.
    873 """
    874 try:
--> 875     return S_ISDIR(self.stat().st_mode)
    876 except OSError as e:
    877     if not _ignore_error(e):

File /opt/conda/envs/ewatercycle2/lib/python3.12/pathlib.py:840, in Path.stat(self, follow_symlinks)
    835 def stat(self, *, follow_symlinks=True):
    836     """
    837     Return the result of the stat() system call on this path, like
    838     os.stat() does.
    839     """
--> 840     return os.stat(self, follow_symlinks=follow_symlinks)

OSError: [Errno 5] Input/output error: '/data/shared/observation/grdc/dailies'
# Choose what you want to run:
Scenario = "ssp370"    # Options: historical, ssp119, ssp126, ssp245, ssp,370, ssp585

Define range & forcing paths#

# Defining period range

# For CMIP trial:
# periods = [
#     [1989, 1995, 2000, 2005, 2010],
#     [1994, 1999, 2004, 2009, 2014]
# ]

# periods = [
#     [1970, 1975, 1980, 1985, 1990, 1995],
#     [1974, 1979, 1984, 1989, 1994, 1999]
# ]

# Future projections 2025-2050:
periods = [
    [2025, 2031, 2036, 2041, 2046],
    [2030, 2035, 2040, 2045, 2050]
]

# Future projections 2050-2075:
# periods = [
#     [2050, 2056, 2061, 2066, 2071],
#     [2055, 2060, 2065, 2070, 2075]
# ]

# Future projections 2075-2099:
# periods = [
#     [2075, 2081, 2086, 2091, 2096],
#     [2080, 2085, 2090, 2095, 2099]
# ]


folder_names = []
forcing_paths = {}

# Generate file names & paths
for i in range(len(periods[0])):
    start_year = periods[0][i]
    end_year = periods[1][i]
    
    folder_name = f'CMIP6-{start_year}-{end_year}'                        
    # forcing_path = Path.home() / "BEP-maxime" / "Workyard" / "forcings" / "CMIP6" / Scenario / folder_name
    # forcing_path = Path.home() / "BEP-maxime" / "Workyard" / "forcings" / "CMIP6" / "SSP" / folder_name
    # forcing_path = Path.home() / "BEP-maxime" / "Workyard" / "forcings" / "CMIP6" / "SSP" / folder_name
    # forcing_path = Path.home() / "BEP-maxime" / "Workyard" / "forcings" / "CMIP6" / "SSP" / folder_name
    forcing_path = Path.home() / "BEP-maxime" / "book" / "thesis_projects" / "BSc" / "2026_Q4_MaximedeBekker_CEG" / "Workyard" / "forcings" / "CMIP6" / Scenario / folder_name

    folder_names.append(folder_name)
    forcing_paths[folder_name] = forcing_path

    forcing_path.mkdir(exist_ok=True, parents=True)

# Call path test
print(forcing_paths[folder_names[0]])

start_year = periods[0][0]
end_year = periods[1][-1]

print(start_year)
print(end_year)
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
25-2030
2025
2050

Generate CMIP6 data#

# Shapefile path
# shape_file = Path.home() / "BEP-maxime" / "Workyard" / "Shapefiles" / "07DA001_basin.shp"
shape_file = Path.home() / "BEP-maxime" / "book" / "thesis_projects" / "BSc" / "2026_Q4_MaximedeBekker_CEG" / "Workyard" / "Shapefiles" / "07DA001_basin.shp"

# Choose CMIP6 dataset
# CMIP_dataset = {'dataset': 'MPI-ESM1-2-HR', 'project': 'CMIP6', 'grid' : 'gn', 'exp': Scenario, 'ensemble': 'r1i1p1f1'}
CMIP_dataset =  {'project': 'CMIP6', 'activity': 'ScenarioMIP', 'exp': Scenario, 'mip': 'day', 'dataset': 'MPI-ESM1-2-HR', 'ensemble': 'r1i1p1f1', 'institute': 'DKRZ', 'grid': 'gn'}

# Generate forcing data
# for i in range(1):
for i in range(len(folder_names)):
    # Change start year to start time so ERA5 can interpret it
    start_time_CMIP = f'{periods[0][i]}-01-01T00:00:00Z'
    end_time_CMIP = f'{periods[1][i]}-12-31T00:00:00Z'

    # Skip if folder is not empty to prevent errors
    if any(forcing_paths[folder_names[i]].iterdir()):
        print(f"Skipping {folder_names[i]} because the folder is not empty.")
        continue

    # Load forcings
    CMIP_forcing = ewatercycle.forcing.sources['LumpedMakkinkForcing'].generate(
        dataset=CMIP_dataset,
        start_time=start_time_CMIP,
        end_time=end_time_CMIP,
        shape=shape_file,
        directory=forcing_paths[folder_names[i]]
    )

    # I need info because my patience is limited
    print(f"Finished {folder_names[i]}")
Skipping CMIP6-2025-2030 because the folder is not empty.
Finished CMIP6-2031-2035
Finished CMIP6-2036-2040
Finished CMIP6-2041-2045
Finished CMIP6-2046-2050

Load CMIP6 data#

CMIP_forcings = {}

for i in range(len(folder_names)):
    # Create new path to find relevant files in the folder
    directory_new = forcing_paths[folder_names[i]] / "work" / "diagnostic" / "script"
    
    # Load forcings
    CMIP_forcings[folder_names[i]] = ewatercycle.forcing.sources["LumpedMakkinkForcing"].load(directory=directory_new)

# Test to see if it works
print(folder_names)
CMIP_forcings[folder_names[0]]
['CMIP6-2025-2030', 'CMIP6-2031-2035', 'CMIP6-2036-2040', 'CMIP6-2041-2045', 'CMIP6-2046-2050']
LumpedMakkinkForcing(start_time='2025-01-01T00:00:00Z', end_time='2030-12-31T00:00:00Z', directory=PosixPath('/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-2025-2030/work/diagnostic/script'), shape=PosixPath('/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-2025-2030/work/diagnostic/script/07DA001_basin.shp'), filenames={'pr': 'CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_pr_gn_2025-2030.nc', 'tas': 'CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_tas_gn_2025-2030.nc', 'rsds': 'CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_rsds_gn_2025-2030.nc', 'evspsblpot': 'Derived_Makkink_evspsblpot.nc'})

Merge CMIP6 data#

# Function to seperate files based on file names and combine same file names from different folders
def combine_variable(var_name):
    datasets = []

    for i in range(len(folder_names)):

        if var_name == "evspsblpot":
            file_name = "Derived_Makkink_evspsblpot.nc"
        else:
            file_name = f"CMIP6_MPI-ESM1-2-HR_day_{Scenario}_r1i1p1f1_{var_name}_gn_{periods[0][i]}-{periods[1][i]}.nc"

        directory = forcing_paths[folder_names[i]] / "work" / "diagnostic" / "script" / file_name

        print(f"Loading {directory}")

        datasets.append(xr.open_dataset(directory))

    combined = xr.concat(datasets, dim="time")
    return combined
# Running function for various file name types
combined_variables = {}
var_names = ["pr", "tas", "rsds", "evspsblpot"]

# combined_path = Path.home() / "BEP-maxime" / "Workyard" / "forcings" / "CMIP6" / Scenario / f'CMIP6-{start_year}-{end_year}'
combined_path = Path.home() / "BEP-maxime" / "book" / "thesis_projects" / "BSc" / "2026_Q4_MaximedeBekker_CEG" / "Workyard" / "forcings" / "CMIP6" / Scenario / f'CMIP6-{start_year}-{end_year}'
combined_path.mkdir(parents=True, exist_ok=True)

for i in range(len(var_names)):
    combined_variables[var_names[i]] = combine_variable(var_names[i])
    combined_variables[var_names[i]].to_netcdf(combined_path / f'combined_CMIP6_{start_year}_{end_year}_{var_names[i]}.nc')
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
25-2030/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_pr_gn_2025-2030.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
31-2035/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_pr_gn_2031-2035.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
36-2040/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_pr_gn_2036-2040.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
41-2045/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_pr_gn_2041-2045.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
46-2050/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_pr_gn_2046-2050.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
25-2030/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_tas_gn_2025-2030.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
31-2035/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_tas_gn_2031-2035.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
36-2040/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_tas_gn_2036-2040.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
41-2045/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_tas_gn_2041-2045.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
46-2050/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_tas_gn_2046-2050.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
25-2030/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_rsds_gn_2025-2030.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
31-2035/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_rsds_gn_2031-2035.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
36-2040/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_rsds_gn_2036-2040.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
41-2045/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_rsds_gn_2041-2045.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
46-2050/work/diagnostic/script/CMIP6_MPI-ESM1-2-HR_day_ssp370_r1i1p1f1_rsds_gn_2046-2050.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
25-2030/work/diagnostic/script/Derived_Makkink_evspsblpot.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
31-2035/work/diagnostic/script/Derived_Makkink_evspsblpot.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
36-2040/work/diagnostic/script/Derived_Makkink_evspsblpot.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
41-2045/work/diagnostic/script/Derived_Makkink_evspsblpot.nc
Loading 
/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forcings/CMIP6/ssp370/CMIP6-20
46-2050/work/diagnostic/script/Derived_Makkink_evspsblpot.nc

Create YAML file#

# Create yaml
forcing_yaml = {
    "start_time": f"{start_year}-01-01T00:00:00Z",
    "end_time": f"{end_year}-12-31T00:00:00Z",
    "shape": str(shape_file),
    "filenames": {
        "pr": f"combined_CMIP6_{start_year}_{end_year}_pr.nc",
        "tas": f"combined_CMIP6_{start_year}_{end_year}_tas.nc",
        "rsds": f"combined_CMIP6_{start_year}_{end_year}_rsds.nc",
        "evspsblpot": f"combined_CMIP6_{start_year}_{end_year}_evspsblpot.nc"
    }}

# Save the YAML file
yaml_file_path = combined_path / "ewatercycle_forcing.yaml"

with open(yaml_file_path, "w") as yaml_file:
    yaml.dump(forcing_yaml, yaml_file, default_flow_style=False)
CMIP6_combined_forcing = ewatercycle.forcing.sources["LumpedMakkinkForcing"].load(directory=combined_path)
print(CMIP6_combined_forcing)
LumpedMakkinkForcing(
    start_time='2025-01-01T00:00:00Z',
    end_time='2050-12-31T00:00:00Z',
    directory=PosixPath('/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/forci
ngs/CMIP6/ssp370/CMIP6-2025-2050'),
    shape=PosixPath('/home/maxime/BEP-maxime/book/thesis_projects/BSc/2026_Q4_MaximedeBekker_CEG/Workyard/Shapefile
s/07DA001_basin.shp'),
    filenames={
        'evspsblpot': 'combined_CMIP6_2025_2050_evspsblpot.nc',
        'pr': 'combined_CMIP6_2025_2050_pr.nc',
        'rsds': 'combined_CMIP6_2025_2050_rsds.nc',
        'tas': 'combined_CMIP6_2025_2050_tas.nc'
    }
)