logo


Analysis of Projected versus Historical Climatology with CORDEX Data#

About#

This notebook is a practical introduction to CORDEX regional climate model data on single levels. CORDEX data are available for 14 regional domains and variable spatial resolutions, from 0.11 x 0.11 degrees up to 0.44 x 0.44 degrees. This workflow will demonstrate how to compute the difference in the air temperature climatology for 2071-2100 (according to a projected scenario) relative to the reference period 1971-2000 in Africa, with a spatial resolution of 0.44 x 0.44 degrees.

We will use the historical experiment to compute the past climatology and repeat the steps to compute the projected climatology according to the scenario RCP4.5. Finally, we will take the difference between the two in order to assess the extent of change. This is called the “delta method” which takes the climate change signal as the difference between the future and past considering that the model biases are the same and can therefore be removed with the subtraction.

This particular example uses data from the Canadian Regional Climate Model CCCma_CanRCM4, with lateral boundary conditions from the Global Climate Model (GCM) CanESM2, also Canadian. This coupled model is referred to as CanRCM4 - CanESM2. See here a full list of models included in the CDS-CORDEX subset.

Note: The workflow uses a specific model, scenario and geographical area, but you can use any model, scenario or domain, and indeed any period or variable. To have a more accurate idea of projected climate change, it is typical to calculate an ensemble of multiple models. Only for the sake of simplicity do we focus on one model.


In general, CORDEX experiments can be divided into three main categories:

  • Evaluation: model simulations for the past with imposed “perfect” lateral boundary conditions following ERA-Interim reanalysis (1979-2015).

  • Historical: Historical experiments cover the period where modern climate observations exist. These experiements show how the RCMs perform for the past climate when forced by the Global Circulation Model (GCM) and can be used as a reference period for comparison with scenario runs for the future.

  • Scenario experiments RCP2.6, RCP4.5, RCP8.5: ensemble of CORDEX climate projection experiments driven by boundary conditions from GCMs using RCP (Representative Concentration Pathways) forcing scenarios. The scenarios used here are RCP 2.6, 4.5 and 8.5, which provide different pathways of the future climate forcing.

Learn here more about CORDEX experiments in the Copernicus Climate Data Store (CDS).

The notebook has the following outline:

  1. Request data from the CDS programmatically with the CDS API

  2. Unzip the downloaded data files

  3. Compute CORDEX historical climatology for reference period 1971-2000

  4. Compute CORDEX projected climatology for target period 2071-2100 based on scenario RCP4.5

  5. Compute and visualize the climate change signal

Data#

This notebook introduces you to CORDEX regional climate model data on single levels. The data used in the notebook has the following specifications:

Data: CORDEX regional climate model data on single levels - Experiment: Historical
Temporal coverage: 1 Jan 1971 to 31 Dec 2000
Spatial coverage: Domain: Africa
Format: NetCDF in zip archives

Data: CORDEX regional climate model data on single levels - Experiment: RCP4.5
Temporal coverage: 1 Jan 2071 to 31 Dec 2100
Spatial coverage: Domain: Africa
Format: NetCDF in zip archives

Run the tutorial via free cloud platforms: Binder Kaggle Colab

Install CDS API via pip#

!pip install cdsapi

Load libraries#

# CDS API
import cdsapi

# Libraries for working with multidimensional arrays
import numpy as np
import xarray as xr

# Libraries for plotting and visualising data
import matplotlib.path as mpath
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.feature as cfeature

# Other libraries (e.g. paths, filenames, zipfile extraction)
from glob import glob
from pathlib import Path
from os.path import basename
import zipfile
import urllib3 
urllib3.disable_warnings() # Disable "InsecureRequestWarning" 
                           # for data download via API

Request data from the CDS programmatically with the CDS API#

The first step is to request data from the Climate Data Store programmatically with the help of the CDS API. Let us make use of the option to manually set the CDS API credentials. First, you have to define two variables: URL and KEY which build together your CDS API key. Below, you have to replace the ######### with your personal CDS key. Please find here your personal CDS key.

URL = 'https://cds.climate.copernicus.eu/api/v2'
KEY = '##################################'
DATADIR = './'

The next step is then to request the data with the help of the CDS API. Below are two blocks of code for the following data requests:

  • Historical experiment: Daily aggregated historical 2m air temperature (from the CanRCM4 - CanESM2 model) from 1971 to 2000 for Africa.

  • RCP4.5 experiment: Daily aggregated RCP4.5 projections of 2m air temperature (from the CanRCM4 - CanESM2 model) from 2071 to 2100 for Africa

Note: Before you run the cells below, the terms and conditions on the use of the data need to have been accepted in the CDS. You can view and accept these conditions by logging into the CDS, searching for the dataset, then scrolling to the end of the Download data section.

c = cdsapi.Client(url=URL, key=KEY)

c.retrieve(
    'projections-cordex-domains-single-levels',
    {
        'format': 'zip',
        'domain': 'africa',
        'experiment': 'historical',
        'horizontal_resolution': '0_44_degree_x_0_44_degree',
        'temporal_resolution': 'daily_mean',
        'variable': '2m_air_temperature',
        'gcm_model': 'cccma_canesm2',
        'rcm_model': 'cccma_canrcm4',
        'ensemble_member': 'r1i1p1',
        'start_year': ['1971', '1976', '1981', '1986', '1991', '1996'],
        'end_year': ['1975', '1980', '1985', '1990', '1995', '2000'],
    },
    f'{DATADIR}1971-2000_cordex_historical_africa.zip')
2022-03-03 22:13:47,034 INFO Welcome to the CDS
2022-03-03 22:13:47,049 INFO Sending request to https://cds.climate.copernicus.eu/api/v2/resources/projections-cordex-domains-single-levels
2022-03-03 22:13:47,142 INFO Request is queued
2022-03-03 22:13:48,199 INFO Request is running
2022-03-03 22:18:05,397 INFO Request is completed
2022-03-03 22:18:05,400 INFO Downloading https://download-0012.copernicus-climate.eu/cache-compute-0012/cache/data6/dataset-projections-cordex-domains-single-levels-19a7525a-dc2e-4e68-bba7-609414960047.zip to ./cordex/1971-2000_cordex_historical_africa.zip (908.2M)
2022-03-03 22:32:25,835 INFO Download rate 1.1M/s                                                                      
Result(content_length=952308195,content_type=application/zip,location=https://download-0012.copernicus-climate.eu/cache-compute-0012/cache/data6/dataset-projections-cordex-domains-single-levels-19a7525a-dc2e-4e68-bba7-609414960047.zip)
c.retrieve(
    'projections-cordex-domains-single-levels',
    {
        'format': 'zip',
        'domain': 'africa',
        'experiment': 'rcp_4_5',
        'horizontal_resolution': '0_44_degree_x_0_44_degree',
        'temporal_resolution': 'daily_mean',
        'variable': '2m_air_temperature',
        'gcm_model': 'cccma_canesm2',
        'rcm_model': 'cccma_canrcm4',
        'ensemble_member': 'r1i1p1',
        'start_year': ['2071', '2076', '2081', '2086', '2091', '2096'],
        'end_year': ['2075', '2080', '2085', '2090', '2095', '2100'],
    },
    f'{DATADIR}2071-2100_cordex_rcp_4_5_africa.zip')
2022-03-03 22:03:00,724 INFO Welcome to the CDS
2022-03-03 22:03:00,730 INFO Sending request to https://cds.climate.copernicus.eu/api/v2/resources/projections-cordex-domains-single-levels
2022-03-03 22:03:00,863 INFO Request is queued
2022-03-03 22:03:01,925 INFO Request is running
2022-03-03 22:09:19,471 INFO Request is completed
2022-03-03 22:09:19,478 INFO Downloading https://download-0000.copernicus-climate.eu/cache-compute-0000/cache/data3/dataset-projections-cordex-domains-single-levels-b19f6b3d-9012-41cd-98fa-c4820a3daeae.zip to ./cordex/2071-2100_cordex_rcp_4_5_africa.zip (908.1M)
2022-03-03 22:12:09,215 INFO Download rate 5.4M/s                                                                      
Result(content_length=952229444,content_type=application/zip,location=https://download-0000.copernicus-climate.eu/cache-compute-0000/cache/data3/dataset-projections-cordex-domains-single-levels-b19f6b3d-9012-41cd-98fa-c4820a3daeae.zip)

Unzip the downloaded data files#

From the Copernicus Climate Data Store, CORDEX data are available in NetCDF format in compressed archives, either zip or tar.gz. For this reason, before we can load any data, we have to unzip or untar the files. Having downloaded the two experiments historical and RCP4.5 as seperate zip files, we can use the functions from the zipfile Python package to extract the content of a zip file. First, we construct a ZipFile() object and second, we apply the function extractall() to extract the content of the zip file.

Below, the same process is repeated for the two zip files downloaded. We see that the actual data files are disseminated in NetCDF. CORDEX files are generally organised in five year periods. For this reason, for each zip file, six NetCDF files in five year periods are extracted.

cordex_zip_paths = glob(f'{DATADIR}*.zip')
for j in cordex_zip_paths:
    with zipfile.ZipFile(j, 'r') as zip_ref:
        zip_ref.extractall(f'{DATADIR}')

Compute CORDEX historical climatology for reference period 1971-2000#

The first step is to load the data of the historical experiment.

We can use the Python library xarray and its function open_mfdataset to read and concatenate multiple NetCDF files at once.

The result is a xarray.Dataset object with four dimensions: bnds, rlat, rlon, time, of which the dimension bnds is not callable.

ds_1971_2000_historical = xr.open_mfdataset(f'{DATADIR}*CanESM2_historical*.nc')
ds_1971_2000_historical
<xarray.Dataset>
Dimensions:       (time: 10950, rlon: 194, rlat: 201, bnds: 2)
Coordinates:
  * time          (time) object 1971-01-01 12:00:00 ... 2000-12-31 12:00:00
  * rlon          (rlon) float64 -24.64 -24.2 -23.76 -23.32 ... 59.4 59.84 60.28
  * rlat          (rlat) float64 -45.76 -45.32 -44.88 ... 41.36 41.8 42.24
    lon           (rlat, rlon) float64 335.4 335.8 336.2 ... 59.4 59.84 60.28
    lat           (rlat, rlon) float64 -45.76 -45.76 -45.76 ... 42.24 42.24
    height        float64 2.0
Dimensions without coordinates: bnds
Data variables:
    rotated_pole  (time) |S1 b'' b'' b'' b'' b'' b'' ... b'' b'' b'' b'' b'' b''
    tas           (time, rlat, rlon) float32 279.3 279.3 279.3 ... 269.3 268.8
    time_bnds     (time, bnds) object 1971-01-01 00:00:00 ... 2001-01-01 00:0...
Attributes: (12/26)
    title:                          CanRCM4 model output prepared for CORDEX ...
    institution:                    CCCma (Canadian Centre for Climate Modell...
    institute_id:                   CCCma
    contact:                        cccma_info@ec.gc.ca
    Conventions:                    CF-1.4
    experiment:                     Historical run driven by CCCma-CanESM2
    ...                             ...
    history:                        created: 2012-06-29 23:52:08 by rcm2nc
    data_licence:                   1) GRANT OF LICENCE - The Government of C...
    creation_date:                  2012-06-29-T23:47:03Z
    c3s_comment:                    This data has been published at ESGF with...
    tracking_id:                    hdl:21.14103/3d35a8d3-105d-451c-88fa-4d46...
    c3s_disclaimer:                 This data has been curated and prepared i...

The next step is to load the data variable tas from the Dataset above. We can load data variables by adding the name of the variable (tas) in square brackets. The resulting object is a xarray.DataArray, which offers additional metadata information about the data, e.g. the unit or long_name of the variable.

tas_1971_2000_historical = ds_1971_2000_historical['tas']
tas_1971_2000_historical
<xarray.DataArray 'tas' (time: 10950, rlat: 201, rlon: 194)>
array([[[279.2568 , 279.26056, 279.26654, ..., 277.77548, 277.74847,
         277.72147],
        [279.64325, 279.64853, 279.6545 , ..., 278.1537 , 278.12518,
         278.10565],
        [280.0147 , 280.02823, 280.03497, ..., 278.54614, 278.5334 ,
         278.52063],
        ...,
        [286.1471 , 286.06378, 286.12534, ..., 265.70358, 265.35767,
         265.00046],
        [285.68558, 285.6721 , 285.67133, ..., 263.67523, 263.27975,
         263.11618],
        [285.18805, 285.2241 , 285.2526 , ..., 261.59732, 261.50504,
         261.1576 ]],

       [[278.82474, 278.78436, 278.7532 , ..., 277.96725, 277.92688,
         277.89856],
        [279.1356 , 279.1129 , 279.0931 , ..., 278.36023, 278.34042,
         278.31915],
        [279.44855, 279.4174 , 279.39758, ..., 278.7688 , 278.7674 ,
         278.75464],
...
        [288.48715, 288.50494, 288.51172, ..., 272.33804, 271.80722,
         271.60403],
        [288.1138 , 288.12567, 288.1392 , ..., 272.37784, 272.03494,
         271.43893],
        [287.74554, 287.77008, 287.80225, ..., 272.35327, 271.90204,
         270.8226 ]],

       [[281.28778, 281.2156 , 281.14505, ..., 278.70392, 278.6908 ,
         278.6793 ],
        [281.68234, 281.60852, 281.54044, ..., 279.0599 , 279.0583 ,
         279.05008],
        [282.05804, 281.99896, 281.93417, ..., 279.44708, 279.4233 ,
         279.40854],
        ...,
        [289.0427 , 289.0788 , 289.12308, ..., 271.27466, 270.85794,
         270.69144],
        [288.75888, 288.79825, 288.85077, ..., 269.81375, 269.76205,
         269.11487],
        [288.44717, 288.49805, 288.5538 , ..., 269.28876, 269.27728,
         268.813  ]]], dtype=float32)
Coordinates:
  * time     (time) object 1971-01-01 12:00:00 ... 2000-12-31 12:00:00
  * rlon     (rlon) float64 -24.64 -24.2 -23.76 -23.32 ... 59.4 59.84 60.28
  * rlat     (rlat) float64 -45.76 -45.32 -44.88 -44.44 ... 41.36 41.8 42.24
    lon      (rlat, rlon) float64 335.4 335.8 336.2 336.7 ... 59.4 59.84 60.28
    lat      (rlat, rlon) float64 -45.76 -45.76 -45.76 ... 42.24 42.24 42.24
    height   float64 2.0
Attributes:
    long_name:      Near-Surface Air Temperature
    standard_name:  air_temperature
    units:          K
    grid_mapping:   rotated_pole
    cell_methods:   time: mean

We will now compute the monthly climatology for the reference period. The monthly climatology represents for each month, and each grid point, the average value over the reference period of 30 years, from 1971 to 2000. We can compute the climatology for each grid point with xarray in two steps: first, all the values for one specific month have to be grouped with the function groupby() and second, we can create the average of each group with the function mean().

The result is a three-dimensional DataArray with the dimensions month, rlat and rlon. The operation has changed the name of the time dimension to month, which has the average near-surface air temperature for each month. Please note that this operation also drops all the attributes (metadata) from the DataArray.

tas_1971_2000_climatology = tas_1971_2000_historical.groupby('time.month').mean()

In a last step, we want to convert the near-surface air temperature values in Kelvin to degrees Celsisus. We can do this by subtracting 273.15 from air temperature values in Kelvin.

tas_1971_2000_climatology_degC = tas_1971_2000_climatology - 273.15

Compute CORDEX future climatology for target period 2071-2100 based on scenario RCP4.5#

Let us now repeat the steps executed for the historical experiment data with the scenario RCP4.5 data. With the projection data for the target period 2071 to 2100, we can compute the projected climatology according to this particular scenario, which can be compared with the historical climatology from the reference period.

Again, we begin with the function open_mfdataset() to load the multiple NetCDF files. The scenario data is also organised into separate files for each 5-year period. These will be concatenated along the time dimension.

ds_2071_2100_projection = xr.open_mfdataset(f'{DATADIR}*CanESM2_rcp45*.nc')

Again we load the data variable tas from the Dataset above into an xarray.DataArray object to facilitate further processing:

tas_2071_2100_projection = ds_2071_2100_projection['tas']

# uncomment the line below before running this cell to view the data characteristics
# tas_2071_2100_projection

Now we will compute the monthly climatology for the target period 2071 to 2100, in the same way as we did for the historical period:

tas_2071_2100_climatology = tas_2071_2100_projection.groupby('time.month').mean()

And finally we convert from degrees Kelvin to Celsius:

tas_2071_2100_climatology_degC = tas_2071_2100_climatology - 273.15

Compute and visualize the climate change signal between the target (2071-2100) and reference (1971-2000) period#

Now, we can compute the difference between the two climatologies. Since both DataArrays have the same structure, we can simply subtract the historical climatology from the projected climatology. The resulting data array has the projected averaged temperature difference for each month.

tas_difference_climatology = tas_2071_2100_climatology - tas_1971_2000_climatology

Now, let us visualize the climate change signal between the target (2071-2100) and reference (1971-2000) period for Africa for June. The plotting function below has five main parts:

  • Initiate the plot: initiate a matplotlib plot with plt.figure and plt.subplot()

  • Plot the 2D map with pcolormesh: plot the xarray data array as 2D map with the function pcolormesh

  • Add coast- and gridlines: add additional items to the plot, such as coast- and gridlines

  • Add title and legend: Set title and legend of the map

  • Save the figure as png file

We can change the value for the variable month and we see that the near-surface air temperature difference is for some regions more pronounced for some months, e.g. June, compared to others.

month=5

# Initiate the plot
fig = plt.figure(figsize=(10,10))
ax = plt.subplot(1,1,1, projection=ccrs.PlateCarree())

# Plot the 2D map with pcolormesh
im = plt.pcolormesh(tas_difference_climatology.rlon, tas_difference_climatology.rlat,
                   tas_difference_climatology[month,:,:], cmap='coolwarm', transform=ccrs.PlateCarree())

# Add coast- and gridlines
ax.coastlines(color='black')
ax.set_extent([-20,60,-40,40], crs=ccrs.PlateCarree())
gl = ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.right_labels = False

# Add title and legend
ax.set_title('CORDEX near surface air temp difference in June: (2071 to 2100 climatology) minus (1971 to 2000 climatology)\n', fontsize=12)
cbar = plt.colorbar(im,fraction=0.046, pad=0.04)
cbar.set_label('\nTemperature difference (Kelvin/Celsius)\n')

# Save the figure
fig.savefig(f'{DATADIR}TAS_1971-2000_2071-2100_june.png', dpi=300)
_images/177fb9ead26ebc499ea8095e20702cdcefb0ed3d232318802d34982c7f4094a5.png

As a last step, we repeat the visualization of climate change signal between the target (2071-2100) and reference (1971-2000) period. But this time, we create the average signal for the whole period, not only one month. For this, we have to compute the average over all months with the function mean().

fig = plt.figure(figsize=(10,10))
ax = plt.subplot(1,1,1, projection=ccrs.PlateCarree())

im = plt.pcolormesh(tas_difference_climatology.rlon, tas_difference_climatology.rlat,
                   tas_difference_climatology.mean('month'), cmap='coolwarm', transform=ccrs.PlateCarree())

ax.coastlines(color='black')
ax.set_extent([-20,60,-40,40], crs=ccrs.PlateCarree())
gl = ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')
gl.top_labels = False
gl.right_labels = False
ax.set_title('CORDEX near surface air temp difference: (2071 to 2100 climatology) minus (1971 to 2000 climatology)\n', fontsize=12)

cbar = plt.colorbar(im,fraction=0.046, pad=0.04)
cbar.set_label('\nTemperature difference (Kelvin/Celsius)\n')

fig.savefig(f'{DATADIR}TAS_1971-2000_2071-2100.png', dpi=300)
_images/b6ea186079c49c772d8ebd5bc3f6d39fb4792f7dda82ce39333491ebf0d5ad5b.png

This project is licensed under APACHE License 2.0. | View on GitHub