⚠️ This is the first operational version of the handbook, but it is still a work in progress and will be heavily updated during 2024! ⚠️

Heatwave hazard assessment based on EURO-CORDEX data analysis#

Description of the assessment#

In this workflow we will assess the heatwave hazard by directly analyzing EURO-CORDEX climate model data, making use of the xclim package to compute heatwave indicators.

Xclim package allows to identify hot and cold spells in the data, calculate heatwave frequency (per year) and total length of heatwaves.

In this workflow the heatwave definition is based on user-defined absolute temperature thresholds for the maximum and minimum daily temperatures (i.e. day and night temperatures), and a mimimum time duration (e.g. 2 or 3 days). Using absolute thresholds for temperature requires input from the user, these limits can be made specific to the region of interest and based on expected health impacts.

The calculation is based on climate data with resolution of 12x12km (EURO-CORDEX), where data is available for 1971-2100.

Advantages of using this methodology:

  • Heatwave calculation based on the minimum and maximum temperature, suitable also for the regions where the temperature can significantly cool down at night.

  • Flexibility of the methodology: the possibility of changing the temperature thresholds according to the needs of the user.

Disadvantages:

  • Computationally intensive due to the large amount of data.

  • Heatwaves are estimated on a yearly basis - not possible to analyze seasonal heatwaves.

Step 1: Preparation work#

Import packages#

import zipfile
import os
from itertools import chain

import rasterio
import cdsapi
import numpy as np
import xarray as xr
from dask.diagnostics import ProgressBar
import cartopy.crs as ccrs
import plotly.graph_objects as go
from ipyleaflet import Map, DrawControl, Marker, LayersControl
import ipywidgets as widgets
from localtileserver import get_leaflet_tile_layer, TileClient
import xclim
from IPython.display import display
# Set host forwarding for remote jupyterlab sessions
if 'JUPYTERHUB_SERVICE_PREFIX' in os.environ:
    os.environ['LOCALTILESERVER_CLIENT_PREFIX'] = f"{os.environ['JUPYTERHUB_SERVICE_PREFIX']}/proxy/{{port}}"

Create directory structure#

# Define the directory for the heatwave workflow preprocess
workflow_folder = 'Heatwave_hazard_XCLIM'
# Define directories for data and results within the previously defined workflow directory
data_dir = os.path.join(workflow_folder,'data')
results_dir = os.path.join(workflow_folder,'results')
# Check if the workflow directory exists, if not, create it along with subdirectories for data and results
if not os.path.exists(workflow_folder):
    os.makedirs(workflow_folder)
    os.makedirs(os.path.join(data_dir))
    os.makedirs(os.path.join(results_dir))

Step 2: Retrieving data#

Download climate data from the CDS#

You can download data from the CDS with the API - look at [How to change your API KEY] and change the KEY below for authentication.

URL = "https://cds-beta.climate.copernicus.eu/api"
KEY = None # put your key here

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

We need to download the daily maximum and daily minimum air temperature at 2m height from the EURO-CORDEX dataset. The code below will download this data for the selected period and RCP scenario.

  • In the code below you can find information about selected locations, resolutions, models, periods, etc.

  • More information about the data you can find here

  • More information about RCP scenarios can be found here

Important: we will only download data for a single combination of GCM and RCM models. Below you can specify which combination you would like to use. For a more reliable result of the assessment, it is recommended to make the analysis for several combinations of GCM-RCM. For an overview of available combinations see the CDS download form and documentation on GCM-RCM models.

gcm_model = 'mpi_m_mpi_esm_lr' # 'cnrm_cerfacs_cm5'  'mohc_hadgem2_es'
rcm_model = 'clmcom_clm_cclm4_8_17'

First we download the daily maximum 2m temperature for 1971-2005.

# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'eurocordex_daily_t2m_eurmax_1971_2005.zip')

# Download the data for the selected time-period and RCP
c.retrieve(
    'projections-cordex-domains-single-levels',
    {
        'format': 'zip',
        'domain': 'europe',
        'experiment': 'historical',
        'horizontal_resolution': '0_11_degree_x_0_11_degree',
        'temporal_resolution': 'daily_mean',
        'variable': 'maximum_2m_temperature_in_the_last_24_hours',
        'gcm_model': gcm_model,
        'rcm_model': rcm_model,
        'ensemble_member': 'r1i1p1',
        'start_year': [
            '1971', '1976', '1981', '1986', '1991', '1996', '2001', 
        ],
        'end_year': [
            '1975', '1980', '1985', '1990', '1995', '2000', '2005',
        ],
    },
    zip_path)

# Unzip the downloaded files so they are ready for computing 
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

# Remove the downloaded zip file after unpacking
os.remove(zip_path)

Now we can download the daily minimum 2m temperature for 1971-2005.

# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'eurocordex_daily_t2m_eurmin_1971_2005.zip')

# Download the data for the selected time-period and RCP
c.retrieve(
    'projections-cordex-domains-single-levels',
    {
        'format': 'zip',
        'domain': 'europe',
        'experiment': 'historical',
        'horizontal_resolution': '0_11_degree_x_0_11_degree',
        'temporal_resolution': 'daily_mean',
        'variable': 'minimum_2m_temperature_in_the_last_24_hours',
        'gcm_model': gcm_model,
        'rcm_model': rcm_model,
        'ensemble_member': 'r1i1p1',
        'start_year': [
            '1971', '1976', '1981', '1986', '1991', '1996','2001',
        ],
        'end_year': [
            '1975', '1980', '1985', '1990', '1995', '2000','2005',
        ],
    },
    zip_path)

# Unzip the downloaded files so they are ready for computing 
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

# Remove the downloaded zip file after unpacking
os.remove(zip_path)

Next we download the daily maximum temperature at 2m for RCP8.5.

You can select the period for download: 2006-2100 by default but you can select other shorter period (e.g. 2041-2070, 2071-2100) to limit the amount of data.

# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'eurocordex_daily_t2m_eurminmax_2006_2100.zip')

# Download the data for the selected time-period and RCP
c.retrieve(
    'projections-cordex-domains-single-levels',
    {
        'format': 'zip',
        'domain': 'europe',
        'experiment': 'rcp_8_5',
        'horizontal_resolution': '0_11_degree_x_0_11_degree',
        'temporal_resolution': 'daily_mean',
        'variable': 'maximum_2m_temperature_in_the_last_24_hours',
        'gcm_model': gcm_model,
        'rcm_model': rcm_model,
        'ensemble_member': 'r1i1p1',
        'start_year': [
            '2006','2011', '2016', '2021', '2026', '2031', '2036', '2041', '2046'],#, '2051', '2056', '2061', '2066', '2071', '2076', '2081', '2086', '2091', '2096'],
        'end_year': [
            '2010','2015', '2020', '2025', '2030', '2035', '2040', '2045', '2050'],#, '2055', '2060', '2065', '2070', '2075', '2080', '2085', '2090', '2095', '2100'],
    },
    zip_path)

# Unzip the downloaded files so they are ready for computing 
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

# Remove the downloaded zip file after unpacking
os.remove(zip_path)

Next we download the daily minimum temperature at 2m for RCP8.5.

You can select the period for download: 2006-2100 by default but you can select other shorter period (e.g. 2041-2070, 2071-2100) to limit the amount of data.

# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'eurocordex_daily_t2m_eurmin_2006_2100.zip')

# Download the data for the selected time-period and RCP
c.retrieve(
    'projections-cordex-domains-single-levels',
    {
        'format': 'zip',
        'domain': 'europe',
        'experiment': 'rcp_8_5',
        'horizontal_resolution': '0_11_degree_x_0_11_degree',
        'temporal_resolution': 'daily_mean',
        'variable': 'minimum_2m_temperature_in_the_last_24_hours',
        'gcm_model': gcm_model,
        'rcm_model': rcm_model,
        'ensemble_member': 'r1i1p1',
        'start_year': [
            '2006','2011', '2016', '2021', '2026', '2031', '2036', '2041', '2046'],#, '2051', '2056', '2061', '2066', '2071', '2076', '2081', '2086', '2091', '2096'],
        'end_year': [
            '2010','2015', '2020', '2025', '2030', '2035', '2040', '2045', '2050'],#, '2055', '2060', '2065', '2070', '2075', '2080', '2085', '2090', '2095', '2100'],
    },
    zip_path)

# Unzip the downloaded files so they are ready for computing 
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

# Remove the downloaded zip file after unpacking
os.remove(zip_path)

Next we download the daily maximum temperature at 2m for RCP4.5.

You can select the period for download: 2006-2100 by default but you can select other shorter period (e.g. 2041-2070, 2071-2100) to limit the amount of data.

# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'eurocordex_daily_t2m_eurmax_rcp45_2006_2100.zip')

# Download the data for the selected time-period and RCP
c.retrieve(
    'projections-cordex-domains-single-levels',
    {
        'format': 'zip',
        'domain': 'europe',
        'experiment': 'rcp_4_5',
        'horizontal_resolution': '0_11_degree_x_0_11_degree',
        'temporal_resolution': 'daily_mean',
        'variable': 'maximum_2m_temperature_in_the_last_24_hours',
        'gcm_model': gcm_model,
        'rcm_model': rcm_model,
        'ensemble_member': 'r1i1p1',
        'start_year': [
            '2006','2011', '2016', '2021', '2026', '2031', '2036', '2041', '2046'],#, '2051', '2056', '2061', '2066', '2071', '2076', '2081', '2086', '2091', '2096'],
        'end_year': [
            '2010','2015', '2020', '2025', '2030', '2035', '2040', '2045', '2050'],#, '2055', '2060', '2065', '2070', '2075', '2080', '2085', '2090', '2095', '2100'],
    },
    zip_path)

# Unzip the downloaded files so they are ready for computing 
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

# Remove the downloaded zip file after unpacking
os.remove(zip_path)

Next we download the daily minimum temperature at 2m for RCP4.5.

You can select the period for download: 2006-2100 by default but you can select other shorter period (e.g. 2041-2070, 2071-2100) to limit the amount of data.

# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'eurocordex_daily_t2m_eurmin_rcp45_2006_2100.zip')

# Download the data for the selected time-period and RCP
c.retrieve(
    'projections-cordex-domains-single-levels',
    {
        'format': 'zip',
        'domain': 'europe',
        'experiment': 'rcp_4_5',
        'horizontal_resolution': '0_11_degree_x_0_11_degree',
        'temporal_resolution': 'daily_mean',
        'variable': 'minimum_2m_temperature_in_the_last_24_hours',
        'gcm_model': gcm_model,
        'rcm_model': rcm_model,
        'ensemble_member': 'r1i1p1',
        'start_year': [
            '2006','2011', '2016', '2021', '2026', '2031', '2036', '2041', '2046'],#, '2051', '2056', '2061', '2066', '2071', '2076', '2081', '2086', '2091', '2096'],
        'end_year': [
            '2010','2015', '2020', '2025', '2030', '2035', '2040', '2045', '2050'],#, '2055', '2060', '2065', '2070', '2075', '2080', '2085', '2090', '2095', '2100'],
    },
    zip_path)

# Unzip the downloaded files so they are ready for computing 
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

# Remove the downloaded zip file after unpacking
os.remove(zip_path)
# rename all files to lowercase to enable easier subsetting of files per model
for file in os.listdir(data_dir):
    os.rename(os.path.join(data_dir,file),os.path.join(data_dir,file.lower()))

Load data and select region of interest#

Load data from your working directory, and set the coordinates system (CRS).

In the code below we will load the data that was downloaded in the previous step, filtering the filenames to select only the data for the specific GCM that was selected (needed in case several models were already downloaded to the same folder, and assuming the RCM stays the same).

# General configuration for data loaders
open_kwargs = {
    "decode_coords": "all",
    "chunks": { "time": 32 }
}

# Historical data maximum temperature
dmaxh = xr.open_mfdataset(f'{data_dir}/tasmax*{gcm_model.replace("_","-")}*historical*.nc', **open_kwargs)
rotated_pole = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
dmaxh.rio.write_crs(rotated_pole, inplace=True)

# Historical data minimum temperature
dminh = xr.open_mfdataset(f'{data_dir}/tasmin*{gcm_model.replace("_","-")}*historical*.nc', **open_kwargs)
rotated_pole = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
dminh.rio.write_crs(rotated_pole, inplace=True)

# Maximum temperature RCP4.5
dmax45 = xr.open_mfdataset(f'{data_dir}/tasmax*{gcm_model.replace("_","-")}*rcp45*.nc', **open_kwargs)
rotated_pole = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
dmax45.rio.write_crs(rotated_pole, inplace=True)

# Minimum temperature RCP4.5
dmin45 = xr.open_mfdataset(f'{data_dir}/tasmin*{gcm_model.replace("_","-")}*rcp45*.nc', **open_kwargs)
rotated_pole = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
dmin45.rio.write_crs(rotated_pole, inplace=True)

# Maximum temperature RCP8.5
dmax85 = xr.open_mfdataset(f'{data_dir}/tasmax*{gcm_model.replace("_","-")}*rcp85*.nc', **open_kwargs)
rotated_pole = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
dmax85.rio.write_crs(rotated_pole, inplace=True)

# Minimum temperature RCP8.5
dmin85 = xr.open_mfdataset(f'{data_dir}/tasmin*{gcm_model.replace("_","-")}*rcp85*.nc', **open_kwargs)
rotated_pole = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
dmin85.rio.write_crs(rotated_pole, inplace=True)

print("Data loaded and crs has been set!")
Data loaded and crs has been set!

Selection of the region of interest#

Tip for the selection of the area: Select at least the area of the region or state where you want to compute (bigger area is better for graphical comparison, for EURO-CORDEX resolution of 12x12km)

  • β€œZoom” to your area with [+] (top left corner)

  • β€œClick” on the rectangle on the left panel

  • β€œSelect your area” by left click and drag

  • And run a code below which β€œautomatically” takes the coordinates of your area and transforms it to rotated pole CRS which is used by the EuroCordex data

# This code plots a map where you need to select a bounding box, as an area for the heatwave estimation
# Create a map centered at a specific location
m = Map(center=(0, 0), zoom=2)
# Create lists to store rectangle coordinates
min_lon_list = []
min_lat_list = []
max_lon_list = []
max_lat_list = []
# Create a DrawControl with rectangle drawing enabled
draw_control = DrawControl(rectangle={'shapeOptions': {'color': '#ff0000'}})
# Add the DrawControl to the map
m.add_control(draw_control)
# Create a text widget to display coordinates
coord_output = widgets.Text(placeholder='Coordinates will appear here', disabled=True)
# Define a function to handle draw events
def handle_draw_polygon(self, action, geo_json):
    if action == 'created':
        if geo_json['geometry']['type'] == 'Polygon':
            # Extract coordinates of the rectangle
            coords = geo_json['geometry']['coordinates'][0]
            # Compute rectangle coordinates (min_lon, min_lat, max_lon, max_lat)
            min_lon, min_lat = min(coord[0] for coord in coords), min(coord[1] for coord in coords)
            max_lon, max_lat = max(coord[0] for coord in coords), max(coord[1] for coord in coords)
            # Update text widget with coordinates
            coord_output.value = f'Coordinates: ({min_lon:.2f}, {min_lat:.2f}), ({max_lon:.2f}, {max_lat:.2f})'
            # Append coordinates to lists
            min_lon_list.append(min_lon)
            min_lat_list.append(min_lat)
            max_lon_list.append(max_lon)
            max_lat_list.append(max_lat)
# Attach the event handler to the DrawControl
draw_control.on_draw(handle_draw_polygon)
# Arrange the map and text widget vertically
display(widgets.VBox([m, coord_output]))

Note: the coordinates will be automatically taken to next steps, you do not need to copy them!

# This code extracts the coordinates from the bounding box selected above and transform it to Rotated pole coordinates
a=min_lon_list + min_lat_list
b=max_lon_list + max_lat_list
# Define the source and destination coordinate reference systems
source_crs = ccrs.PlateCarree()  # WGS84 CRS
#dest_crs = ccrs.RotatedPole(pole_longitude=190, pole_latitude=39.25, central_rotated_longitude=0)
dest_crs = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
# Example coordinates in WGS 84
#wgs84_coords = [(18.528168, 49.076529), (18.965041, 49.349247)]
wgs84_coords = [(a), (b)]
# Perform the transformation for each pair of coordinates
rotated_pole_coords = []
for lon, lat in wgs84_coords:
    rotated_pole_coords.append(dest_crs.transform_point(lon, lat, source_crs))
print("WGS 84 Coordinates:", wgs84_coords)
print("Rotated Pole Coordinates:", rotated_pole_coords)
WGS 84 Coordinates: [[4.196777, 51.99841], [4.43161, 52.12843]]
Rotated Pole Coordinates: [(-8.452396779175887, 2.037647168587727), (-8.286523211884221, 2.1389936593755525)]
# This code will create a list from the coordinates in WGS84
# Original list of tuples
list_of_tuples = wgs84_coords
# Convert the list of tuples to a single tuple
bbox84 = tuple(chain.from_iterable(list_of_tuples))
print(bbox84)
(4.196777, 51.99841, 4.43161, 52.12843)
# This code will create a list from the coordinates in rotated pole
# Original list of tuples
list_of_tuples2 = rotated_pole_coords
# Convert the list of tuples to a single tuple
bbox = tuple(chain.from_iterable(list_of_tuples2))
print(bbox)
(-8.452396779175887, 2.037647168587727, -8.286523211884221, 2.1389936593755525)

Clip data to the region of interest#

The bouding box coordinates will be automatically copied after a selection of the bbox

# This code crops the downloaded data from CDS to your selected bounding box (selected area)
dmaxch=dmaxh.rio.clip_box(minx=bbox[0], miny=bbox[1], maxx=bbox[2], maxy=bbox[3])  
dminch=dminh.rio.clip_box(minx=bbox[0], miny=bbox[1], maxx=bbox[2], maxy=bbox[3])
dmaxc45=dmax45.rio.clip_box(minx=bbox[0], miny=bbox[1], maxx=bbox[2], maxy=bbox[3])  
dminc45=dmin45.rio.clip_box(minx=bbox[0], miny=bbox[1], maxx=bbox[2], maxy=bbox[3])  
dmaxc85=dmax85.rio.clip_box(minx=bbox[0], miny=bbox[1], maxx=bbox[2], maxy=bbox[3])  
dminc85=dmin85.rio.clip_box(minx=bbox[0], miny=bbox[1], maxx=bbox[2], maxy=bbox[3])  
# If the size of the area is relatively small (city-scale) we can load the data into memory for faster calculation
dmaxch.load()
dminch.load()
dmaxc45.load()
dminc45.load()
dmaxc85.load()
dminc85.load()
<xarray.Dataset> Size: 526kB
Dimensions:                     (time: 16436, rlat: 2, rlon: 3)
Coordinates:
  * time                        (time) datetime64[ns] 131kB 2006-01-01T12:00:...
  * rlat                        (rlat) float64 16B 2.035 2.145
  * rlon                        (rlon) float64 24B -8.465 -8.355 -8.245
    lat                         (rlat, rlon) float32 24B 51.99 52.01 ... 52.14
    lon                         (rlat, rlon) float32 24B 4.177 4.353 ... 4.496
    height                      float64 8B 2.0
    rotated_latitude_longitude  int32 4B 0
Data variables:
    tasmin                      (time, rlat, rlon) float32 394kB 273.4 ... 268.3
Attributes: (12/30)
    institution:                    Climate Limited-area Modelling Community ...
    institute_id:                   CLMcom
    experiment_id:                  rcp85
    source:                         CLMcom-CCLM4-8-17
    model_id:                       CLMcom-CCLM4-8-17
    contact:                        cordex-cclm@dkrz.de
    ...                             ...
    table_id:                       Table day (Sept 2013) 0cf1782745489246c9f...
    modeling_realm:                 atmos
    realization:                    1
    cmor_version:                   2.9.1
    tracking_id:                    hdl:21.14103/9f40a965-de36-42c4-9f19-cf4d...
    c3s_disclaimer:                 This data has been produced in the contex...

Step 3: Calculation of the heatwave occurrence based on the XCLIM methodology#

Convert the maximum temperature data from Kelvin (K) to Celsius (Β°C)

# Max Temperature historical
dmaxch = dmaxch['tasmax']
tasmaxh = dmaxch - 273.15
tasmaxh = tasmaxh.assign_attrs(dmaxch.attrs)
tasmaxh.attrs['units']='Β°C'

# Min Temperature historical
dminch = dminch['tasmin']
tasminh = dminch - 273.15
tasminh = tasminh.assign_attrs(dminch.attrs)
tasminh.attrs['units']='Β°C'

# Max Temperature rcp 45
dmaxc45 = dmaxc45['tasmax']
tasmax45 = dmaxc45 - 273.15
tasmax45 = tasmax45.assign_attrs(dmaxc45.attrs)
tasmax45.attrs['units']='Β°C'

# Min Temperature rcp 45
dminc45 = dminc45['tasmin']
tasmin45 = dminc45 - 273.15
tasmin45 = tasmin45.assign_attrs(dminc45.attrs)
tasmin45.attrs['units']='Β°C'

# Max Temperature rcp 85
dmaxc85 = dmaxc85['tasmax']
tasmax85 = dmaxc85 - 273.15
tasmax85 = tasmax85.assign_attrs(dmaxc85.attrs)
tasmax85.attrs['units']='Β°C'

# Min Temperature rcp 85
dminc85 = dminc85['tasmin']
tasmin85 = dminc85 - 273.15
tasmin85 = tasmin85.assign_attrs(dminc85.attrs)
tasmin85.attrs['units']='Β°C'

XCLIM is an operational Python library for climate services, providing numerous climate-related indicator tools with a framework for constructing custom climate indicators, statistical downscaling and bias adjustment of climate model simulations, as well as climate model ensemble analysis tools. XCLIM is built using xarray and can seamlessly benefit from the parallelization handling provided by dask. Its objective is to make it as simple as possible for users to perform typical climate data analyses. Leveraging xarray and dask, users can easily bias-adjust climate simulations over large spatial domains or compute indices from large climate datast [XCLIM]

Examples of the code for the computing of the climate indices.

xclim.indices.heat_wave_max_length(tasmin, tasmax, thresh_tasmin='22.0 degC', thresh_tasmax='30 degC', window=3, freq='YS', op='>', resample_before_rl=True)

In the example above (from XCLIM documentation) the thresholds of 22Β°C for night temperatures and 30Β°C for day temperatures were selected by Health Canada professionals, following a temperature–mortality analysis. These absolute temperature thresholds characterize the occurrence of hot weather events that can result in adverse health outcomes for Canadian communities [Casati et al., 2013]. For other regions these thresholds may need to be selected differently.

Loading data#

# This code selects the periods and RCPs for the computation of the climate indices by the XCLIM
# Max T
tasmaxh = tasmaxh.sel(time=slice("1971-01-01", "2005-12-31"))
tasmaxp_45 = tasmax45.sel(time=slice("2006-01-01", "2100-12-31"))
tasmaxp_85 = tasmax85.sel(time=slice("2006-01-01", "2100-12-31"))
# Min T
tasminh = tasminh.sel(time=slice("1971-01-01", "2005-12-31"))
tasminp_45 = tasmin45.sel(time=slice("2006-01-01", "2100-12-31"))
tasminp_85 = tasmin85.sel(time=slice("2006-01-01", "2100-12-31"))

Defining the temperature and duration tresholds#

Treshold temperatures based on the impact on human health:

  1. Based on the preferred studies e.g.: A multi-city epidemiologic study in Europe found that the mortality threshold for heat effects was 29.4 Β°C for Mediterranean cities and 23.3 Β°C for the north European cities [Maximum temperature threshold for EU].

  2. Or use the thresholds for minimum and maximum temperature given by your national meteorological or health service.

We select three climate indices from the Xclim package:

  1. Heatwave index

  2. Heatwave frequency

  3. Heatwave total length

But you can also choose others, all climate indices are listed on Xclim website [source]

Calculating the heatwave index [source]#

We need to specify the temperature threshold for maximum daily air temperature (in degC) and the duration threshold (in days).

thresh_tempmax = '25.0 degC' # maximum daily temperature threshold
thresh_window = 5 #duration threshold in days
# This code uses the XCLIM package for the computation of the climate indices
# heatwave index
# RCP 4.5
HWIp45 = xclim.indices.heat_wave_index(tasmaxp_45, thresh=thresh_tempmax, window=thresh_window, freq='YS', op='>', resample_before_rl=True)
# RCP 8.5
HWIp85 = xclim.indices.heat_wave_index(tasmaxp_85, thresh=thresh_tempmax, window=thresh_window, freq='YS', op='>', resample_before_rl=True)
# Historical
HWIh = xclim.indices.heat_wave_index(tasmaxh, thresh=thresh_tempmax, window=thresh_window, freq='YS', op='>', resample_before_rl=True)

Calculating the heatwave frequency [source]#

We need to specify:

  • Temperature threshold for minimum daily air temperature (in degC)

  • Temperature threshold for maximum daily air temperature (in degC)

  • Duration threshold (in days)

thresh_tasmin = '22.0 degC' # threshold for minimum daily air temperature
thresh_tasmax = '28.0 degC' # threshold for maximum daily air temperature
thresh_window_hfreq = 3 #duration threshold in days
# Heatwave frequency
# RCP 4.5
HWFp45 = xclim.indices.heat_wave_frequency(tasminp_45, tasmaxp_45, thresh_tasmin=thresh_tasmin, thresh_tasmax=thresh_tasmax, window=thresh_window_hfreq, freq='YS', op='>', resample_before_rl=True)
# RCP 8.5
HWFp85 = xclim.indices.heat_wave_frequency(tasminp_85, tasmaxp_85, thresh_tasmin=thresh_tasmin, thresh_tasmax=thresh_tasmax, window=thresh_window_hfreq, freq='YS', op='>', resample_before_rl=True)
#Historical
HWFh = xclim.indices.heat_wave_frequency(tasminh, tasmaxh, thresh_tasmin=thresh_tasmin, thresh_tasmax=thresh_tasmax, window=thresh_window_hfreq, freq='YS', op='>', resample_before_rl=True)

# Data correction: set the right type
HWFp45 = HWFp45.astype(np.float64)
HWFp85 = HWFp85.astype(np.float64)
HWFh = HWFh.astype(np.float64)

Calculating heatwave total length [source]#

For this indicator we use the same thresholds as for the heatwave frequency indicator (above).

# heatwave total length
# RCP 4.5
HWTLp45=xclim.indices.heat_wave_total_length(tasminp_45, tasmaxp_45, thresh_tasmin=thresh_tasmin, thresh_tasmax=thresh_tasmax, window=thresh_window, freq='YS', op='>', resample_before_rl=True)
# RCP 8.5
HWTLp85=xclim.indices.heat_wave_total_length(tasminp_85, tasmaxp_85, thresh_tasmin=thresh_tasmin, thresh_tasmax=thresh_tasmax, window=thresh_window, freq='YS', op='>', resample_before_rl=True)
# Historical
HWTLh=xclim.indices.heat_wave_total_length(tasminh, tasmaxh, thresh_tasmin=thresh_tasmin, thresh_tasmax=thresh_tasmax, window=thresh_window, freq='YS', op='>', resample_before_rl=True)
# This code sets the rotated pole CRS to the calculated data
## CRS for Eurocordex data
rotated_pole = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
HWIp45.rio.write_crs(rotated_pole, inplace=True)
HWIp85.rio.write_crs(rotated_pole, inplace=True)
HWIh.rio.write_crs(rotated_pole, inplace=True)
HWFp45.rio.write_crs(rotated_pole, inplace=True)
HWFp85.rio.write_crs(rotated_pole, inplace=True)
HWFh.rio.write_crs(rotated_pole, inplace=True)
HWTLp45.rio.write_crs(rotated_pole, inplace=True)
HWTLp85.rio.write_crs(rotated_pole, inplace=True)
HWTLh.rio.write_crs(rotated_pole, inplace=True)
print("CRS written!")
CRS written!

The next step will perform the calculations and save the result to the disk. Depending on your machine’s computational capacity, this step may take some time.

# This code saves the computed indices in the disk in raster tif format 
with ProgressBar():
    print("Heatwave index")
    HWIp45.rio.to_raster(raster_path=f'{data_dir}/HWIp45.tif')
    HWIp85.rio.to_raster(raster_path=f'{data_dir}/HWIp85.tif')
    HWIh.rio.to_raster(raster_path=f'{data_dir}/HWIh.tif')
    print("Heatwave frequency")
    HWFp45.rio.to_raster(raster_path=f'{data_dir}/HWFp45.tif')
    HWFp85.rio.to_raster(raster_path=f'{data_dir}/HWFp85.tif')
    HWFh.rio.to_raster(raster_path=f'{data_dir}/HWFh.tif')
    print("Heatwave total length")
    HWTLp45.rio.to_raster(raster_path=f'{data_dir}/HWTLp45.tif')
    HWTLp85.rio.to_raster(raster_path=f'{data_dir}/HWTLp85.tif')
    HWTLh.rio.to_raster(raster_path=f'{data_dir}/HWTLh.tif')
Heatwave index
Heatwave frequency
Heatwave total length

Select location for plotting#

# This code reprojects the computed indices to the WGS84 (EPSG:4326)
# Define the files and their corresponding variables
file_variables = [
    ("HWIp45.tif", "HWIp45"),
    ("HWIp85.tif", "HWIp85"),
    ("HWIh.tif", "HWIh"),
    ("HWFp45.tif", "HWFp45"),
    ("HWFp85.tif", "HWFp85"),
    ("HWFh.tif", "HWFh"),
    ("HWTLp45.tif", "HWTLp45"),
    ("HWTLp85.tif", "HWTLp85"),
    ("HWTLh.tif", "HWTLh")
]
# Define the target CRS
target_crs = 'EPSG:4326'  # WGS84
# Iterate over each file and variable
for filename, var_name in file_variables:
    # Open the raster dataset
    data = xr.open_dataset(f'{data_dir}/{filename}')['band_data']   
    # Write the CRS information to the dataset
    data.rio.write_crs(rotated_pole, inplace=True)  
    # Reproject the clipped raster dataset to the target CRS
    reprojected_data = data.rio.reproject(target_crs)
    # Replace NaN values
    reprojected_data = reprojected_data.where(reprojected_data != 1.7976931348623157e+308, np.nan)   
    # Optionally, you can save the reprojected raster to a new file
    reprojected_data.rio.to_raster(raster_path=f'{data_dir}/{var_name}_reprojected.tif') 
    # Calculate mean and save to file
    mean_data = reprojected_data.mean(dim='band', skipna=True, keep_attrs=True)
    mean_data.rio.to_raster(raster_path=f'{data_dir}/{var_name}_m.tif')

With this code, you select the point on the map, for which you want to plot the bar and line graphs (select only one point, if you accidentally select multiple points, run this code (below) again and select the point again).

# This code plots the map where we need to select the point for the plot of the graphs for computed indices
# First, create a tile server from local raster file
client1 = TileClient(f'{data_dir}/HWIp85_m.tif')
# Create ipyleaflet tile layer from that server
t1 = get_leaflet_tile_layer(client1, colormap='Reds', opacity=0.5, nodata=0, name='heatwave index')
m = Map(center=client1.center(), zoom=client1.default_zoom)
m.add(t1)
control = LayersControl(position='topright')
m.add_control(control)
# Create lists to store point coordinates
point_lon_list = []
point_lat_list = []
# Create a DrawControl with point drawing enabled
draw_control = DrawControl(marker={'shapeOptions': {'color': '#FF0000'}})
# Add the DrawControl to the map
m.add_control(draw_control)
# Create a text widget to display coordinates
coord_output = widgets.Text(placeholder='Coordinates will appear here', disabled=True)
# Define a function to handle draw events
def handle_draw_point(self, action, geo_json):
    if action == 'created':
        if geo_json['geometry']['type'] == 'Point':
            # Extract coordinates of the point
            lon, lat = geo_json['geometry']['coordinates']
            # Update text widget with coordinates
            coord_output.value = f'Coordinates: ({lon:.2f}, {lat:.2f})'
            # Append coordinates to lists
            point_lon_list.append(lon)
            point_lat_list.append(lat)
            # Create and add a marker to the map
            marker = Marker(location=(lat, lon))
            m.add_layer(marker)
# Attach the event handler to the DrawControl
draw_control.on_draw(handle_draw_point)
# Arrange the map and text widget vertically
display(widgets.VBox([m, coord_output]))
c:\Users\aleksand\AppData\Local\anaconda3\envs\climaax_heatwaves\Lib\site-packages\localtileserver\widgets.py:90: UserWarning: The `cmap` keyword argument is deprecated. Please use `colormap` instead.
  warnings.warn(
c:\Users\aleksand\AppData\Local\anaconda3\envs\climaax_heatwaves\Lib\site-packages\rio_tiler\models.py:140: RuntimeWarning: invalid value encountered in cast
  return array.astype(out_dtype)
c:\Users\aleksand\AppData\Local\anaconda3\envs\climaax_heatwaves\Lib\site-packages\rio_tiler\models.py:140: RuntimeWarning: invalid value encountered in cast
  return array.astype(out_dtype)
c:\Users\aleksand\AppData\Local\anaconda3\envs\climaax_heatwaves\Lib\site-packages\rio_tiler\models.py:140: RuntimeWarning: invalid value encountered in cast
  return array.astype(out_dtype)
c:\Users\aleksand\AppData\Local\anaconda3\envs\climaax_heatwaves\Lib\site-packages\rio_tiler\models.py:140: RuntimeWarning: invalid value encountered in cast
  return array.astype(out_dtype)
c:\Users\aleksand\AppData\Local\anaconda3\envs\climaax_heatwaves\Lib\site-packages\rio_tiler\models.py:140: RuntimeWarning: invalid value encountered in cast
  return array.astype(out_dtype)
c:\Users\aleksand\AppData\Local\anaconda3\envs\climaax_heatwaves\Lib\site-packages\rio_tiler\models.py:140: RuntimeWarning: invalid value encountered in cast
  return array.astype(out_dtype)
c:\Users\aleksand\AppData\Local\anaconda3\envs\climaax_heatwaves\Lib\site-packages\rio_tiler\models.py:140: RuntimeWarning: invalid value encountered in cast
  return array.astype(out_dtype)

In the picture above, we can find the plotted raster data in the area that we select by bbox at the beginning of the code.

The shown data is based on the calculated heatwave index. This visual comparison is intended to help understand which locations in your selected area will suffer more from heat (dark red will be most affected).

For the visual comparison of the number of heat occurrences you need to select a point (flag) from the left panel and select a pixel that you prefer. The skewness of the rasters was caused by the transformation of the coordinates system.

# This code extracts the values from the computed indices for the selected pixel 
def get_pixel_values(raster_file, lat, lon):
    with rasterio.open(raster_file) as src:
        # Transform lat/lon coordinates to pixel coordinates
        row, col = src.index(lon, lat)
        
        # Read pixel values from each band
        pixel_values = [band[row, col] for band in src.read()]
        
        # Return pixel values for each band
        return pixel_values

# Example usage
raster_file_hwi45 = f'{data_dir}/HWIp45_reprojected.tif'
raster_file_hwi85 = f'{data_dir}/HWIp85_reprojected.tif'
raster_file_hwih = f'{data_dir}/HWIh_reprojected.tif'
raster_file_hwf45 = f'{data_dir}/HWFp45_reprojected.tif'
raster_file_hwf85 = f'{data_dir}/HWFp85_reprojected.tif'
raster_file_hwfh = f'{data_dir}/HWFh_reprojected.tif'
raster_file_hwtl45 = f'{data_dir}/HWTLp45_reprojected.tif'
raster_file_hwtl85 = f'{data_dir}/HWTLp85_reprojected.tif'
raster_file_hwtlh = f'{data_dir}/HWTLh_reprojected.tif'
# extract lat lon from map
lat = point_lat_list[0]
lon = point_lon_list[0]
# extract values for selected pixel 
pixel_values_hwi45 = get_pixel_values(raster_file_hwi45, lat, lon)
pixel_values_hwi85 = get_pixel_values(raster_file_hwi85, lat, lon)
pixel_values_hwih = get_pixel_values(raster_file_hwih, lat, lon)
pixel_values_hwf45 = get_pixel_values(raster_file_hwf45, lat, lon)
pixel_values_hwf85 = get_pixel_values(raster_file_hwf85, lat, lon)
pixel_values_hwfh = get_pixel_values(raster_file_hwfh, lat, lon)
pixel_values_hwtl45 = get_pixel_values(raster_file_hwtl45, lat, lon)
pixel_values_hwtl85 = get_pixel_values(raster_file_hwtl85, lat, lon)
pixel_values_hwtlh = get_pixel_values(raster_file_hwtlh, lat, lon)

Plot the heatwave index (HWI)#

Number of days that are part of a heatwave, defined as five or more consecutive days over a threshold of 25 degrees (or with other thresholds if specified differently above!).

# This code plots the data with heatwave index for the selected pixel 
# Sample data
b = pixel_values_hwih + pixel_values_hwi45
t = pixel_values_hwih + pixel_values_hwi85 
# Define years
yearh = list(range(1971, 2005))
yearp1 = list(range(2006, 2101))
year = yearh + yearp1
# Create figure
fig = go.Figure()
# Add bar trace for 'b'
fig.add_trace(go.Bar(x=year, y=b, name='RCP4.5', marker_color='blue'))
# Add bar trace for 't'
fig.add_trace(go.Bar(x=year, y=t, name='RCP8.5', opacity=0.5, marker_color='red'))

# Add information on models
fig.add_annotation(xref='paper',yref='paper', y=-0.2,x=1, showarrow=False,
            text=f'GCM-RCM: {gcm_model} - {rcm_model}')

# Update layout
fig.update_layout(
    title={'text':f'HEATWAVE INDEX<br>Number of days that are part of a heatwave defined as {thresh_window} or more consequtive days over a threshold of {thresh_tempmax}<br>Location: lat {lat:.02f}\u00B0, lon {lon:.02f}\u00B0','x':0.5,'xanchor': 'center'},
    xaxis_title='Year',
    yaxis_title='Number of days',
    legend=dict(x=0, y=1, traceorder='normal', orientation='h'),
    width=1100,
    height=500)

Plot the heatwave frequency#

Number of heatwaves over a given period. A heatwave is defined as an event where the minimum and maximum daily temperature both exceed specific thresholds over a minimum number of days. The number of the heatwaves depends strictly on the threshold which we selected above.

# This code plots the data with heatwave frequency for the selected pixel 
# Sample data
b = pixel_values_hwfh + pixel_values_hwf45
t = pixel_values_hwfh + pixel_values_hwf85 
# Define years
yearh = list(range(1971, 2001))
yearp1 = list(range(2011, 2101))
year = yearh + yearp1
# Create figure
fig = go.Figure()
# Add bar trace for 'b'
fig.add_trace(go.Bar(x=year, y=b, name='RCP4.5', marker_color='blue'))
# Add bar trace for 't'
fig.add_trace(go.Bar(x=year, y=t, name='RCP8.5', opacity=0.5, marker_color='red'))
# Update layout

# Add information on models
fig.add_annotation(xref='paper',yref='paper', y=-0.2,x=1, showarrow=False,
            text=f'GCM-RCM: {gcm_model} - {rcm_model}')

fig.update_layout(
    title={'text':f'HEATWAVE FREQUENCY<br>based on temperature thresholds for day and night: {thresh_tasmax} and {thresh_tasmin}, for at least {thresh_window_hfreq} consequtive days<br>Location: lat {lat:.02f}\u00B0, lon {lon:.02f}\u00B0','x':0.5,'xanchor': 'center'},
    xaxis_title='Year',
    yaxis_title='Number of heatwaves per year',
    legend=dict(x=0, y=1, traceorder='normal', orientation='h'),
    width=1100,
    height=500)

Plot the heatwave total length#

Total length of heatwaves over a given period. A heatwave is defined as an event where the minimum and maximum daily temperature both exceeds specific thresholds over a minimum number of days. This the sum of all days in such events.

# This code plots the data with heatwave total length for the selected pixel 
# Sample data
b = pixel_values_hwtlh + pixel_values_hwtl45
t = pixel_values_hwtlh + pixel_values_hwtl85 
# Define years
yearh = list(range(1971, 2001))
yearp1 = list(range(2011, 2101))
year = yearh + yearp1
# Create figure
fig = go.Figure()
# Add bar trace for 'b'
fig.add_trace(go.Bar(x=year, y=b, name='RCP4.5', marker_color='blue'))
# Add bar trace for 't'
fig.add_trace(go.Bar(x=year, y=t, name='RCP8.5', opacity=0.5, marker_color='red'))

# Add information on models
fig.add_annotation(xref='paper',yref='paper', y=-0.2,x=1, showarrow=False,
            text=f'GCM-RCM: {gcm_model} - {rcm_model}')

# Update layout
fig.update_layout(
    title={'text':f'Total number of heatwave days<br>based on temperature thresholds for day and night: {thresh_tasmax[:4]}\u00B0C and {thresh_tasmin[:4]}\u00B0C, for at least {thresh_window_hfreq} consequtive days<br>Location: lat {lat:.02f}\u00B0, lon {lon:.02f}\u00B0','x':0.5,'xanchor': 'center'},
    xaxis_title='Year',
    yaxis_title='Total number of days during heatwave events',
    legend=dict(x=0, y=1, traceorder='normal', orientation='h'),
    width=1100,
    height=500)

Conclusions for the XCLIM results#

We calculated three climate indices connected to extreme temperature: heatwave index, heatwave frequency, and heatwave total length for two climate scenarios (RCP4.5 and RCP8.5) based on the EURO-CORDEX data.

To interpret the plotted results, it is important to remember that the result is dependent on the GCM-RCM combination of climate models that is used in the analysis. This workflow is based on a single model combination, and using a different combination will lead to a different distribution of temperature for specific years. However, the general trend over the entire analysis period will be similar. It is recommended to perform this analysis for several GCM-RCM combinations to check how much the result is dependent on the selection of the models.

The plotted data with indicator values per year should not be used as a forecast of heatwaves for specific years, as this is strongly dependent on the model realization. For policy and planning it may be useful to calculate average indices (e.g. average heatwave occurence) over periods of 20-30 years and compare to historical periods of the same length (e.g. 2040-2060 vs. 1990-2010).

The plotted results are obtained for a single location within the EURO-CORDEX dataset, which corresponds to a 12x12 km cell. The data within this cell does not take into account the finer resolution effects such as the urban heat island effect, and does not take into account future land use changes.

References#

Authors#