⚠️ 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 using EuroHEAT methodology#

Description of the assessment#

Heatwave based on the EuroHEAT project approach#

The heatwave hazard assessment presented in this page is directly based on an existing dataset of heat wave days on CDS. You can also view these results directly on the CDS using an application that shows aggregated results for the NUTS3 administrative regions. This notebook provides an example on how to extract this data for any grid cell of the dataset (resolution of 12 x 12km) in the EU.

Methodology is based on the EuroHEAT project (source) with two options for defining heatwaves:

  • Health-related EU-wide definition: for the summer period of June to August, heat waves were defined as days in which the maximum apparent temperature (Tappmax) exceeds the threshold (90th percentile of Tappmax for each month) and the minimum temperature (Tmin) exceeds its threshold (90th percentile of Tmin for each month) for at least two days. The apparent temperature is a measure of relative discomfort due to combined heat and high humidity, developed based on physiological studies on evaporative skin cooling. It can be calculated as a combination of air and dew point temperature.

  • National heat-wave definition: each nation uses a different methodology e.g. for Belgium, a heatwave-day is a day on which both the three-day running daily minimum and maximum temperature exceeds a threshold of 18.2 °C and 29.6 °C between April and September. This condition has to be met for three consecutive days. For this example, the source of the heat wave definition is the Belgian federal public health agency. See this page for more national definitions.

Advantages of using this methodology and dataset:

  • Pre-computed indicators - less data that needs to be downloaded and processed (in comparison with the Peseta-IV and Xclim hazard workflows).

  • Data available on a 12x12km grid for years 1986-2085 for the whole EU.

  • The heatwave definition is based on both maximum and minimum daily temperatures, so it accounts for the effect of minimum temperature on the severity of a heat wave (e.g. when temperature significantly drops during the nights thus making the heatwave less severe).

Disadvantages:

  • No possibility of changing thresholds for temperature and number of days in the heatwave definition.

  • Only contains values aggregated per year.

  • National heat wave definitions are not available for all EU states (only a selection).

Step 1. Preparation work#

Import packages#

import zipfile
import os
from pathlib import Path
import cdsapi
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
from rasterio.crs import CRS
from pylab import pcolormesh, show, colorbar, plot, title, legend, subplot, savefig
from matplotlib import pyplot
from rasterio.plot import show
import rasterio as rio
import plotly.graph_objects as go
from sklearn.linear_model import LinearRegression
from ipyleaflet import Map, DrawControl, Marker, LayersControl
import ipywidgets as widgets
import leafmap.leafmap as leafmap
from localtileserver import get_leaflet_tile_layer, TileClient

Create a directory structure#

# Define the directory for storing data
workflow_folder = 'Heatwave_hazard_EuroHEAT'

# 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 on heat wave occurence#

Download data from the CDS#

  • You can download data from the CDS with the API. If this is the first time you do this, you will need to register on CDS and set your API key. Check [How to change your API KEY] and change the KEY in the code cell below.

URL = "https://cds.climate.copernicus.eu/api/v2"
KEY = "key" ### put here your key!!!

Heat wave data based on National definitions from EuroHEAT [source]#

# Heat waves and cold spells in Europe derived from climate projections
c = cdsapi.Client(url=URL, key=KEY)
c.retrieve(
    'sis-heat-and-cold-spells',
    {
        'variable': 'heat_wave_days',
        'definition': 'country_related',
        'experiment': [
            'rcp4_5', 'rcp8_5',
        ],
        'ensemble_statistic': [
            'ensemble_members_average', 'ensemble_members_standard_deviation',
        ],
        'format': 'zip',
    },
   f"{data_dir}/heat_spells_country_1986_2085.zip")
2024-06-28 13:19:03,563 INFO Welcome to the CDS
2024-06-28 13:19:03,564 INFO Sending request to https://cds.climate.copernicus.eu/api/v2/resources/sis-heat-and-cold-spells
2024-06-28 13:19:03,656 INFO Request is queued
2024-06-28 13:19:04,714 INFO Request is running
2024-06-28 13:19:17,147 INFO Request is completed
2024-06-28 13:19:17,148 INFO Downloading https://download-0006-clone.copernicus-climate.eu/cache-compute-0006/cache/data1/dataset-sis-heat-and-cold-spells-32ed2bf0-d6e7-4d5e-8cc0-edff0c5d0322.zip to Heatwave_hazard_EuroHEAT\data/heat_spells_country_1986_2085.zip (32M)
2024-06-28 13:19:21,749 INFO Download rate 7M/s     
Result(content_length=33605837,content_type=application/zip,location=https://download-0006-clone.copernicus-climate.eu/cache-compute-0006/cache/data1/dataset-sis-heat-and-cold-spells-32ed2bf0-d6e7-4d5e-8cc0-edff0c5d0322.zip)
zip_path = os.path.join(data_dir, 'heat_spells_country_1986_2085.zip')
# Extract from zip file
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

Plotting the data#

Prepare data for plotting#

  • In this step we prepare the data for plotting. Plotting will later be done for a selected pixel on the EU map.

# This code loads a data form data dir and sets the CRS
hwd45h = xr.open_dataset(f'{data_dir}/HWD_EU_health_rcp45_mean_v1.0.nc', decode_coords='all')
hwd85h = xr.open_dataset(f'{data_dir}/HWD_EU_health_rcp85_mean_v1.0.nc', decode_coords='all')
hwd45n = xr.open_dataset(f'{data_dir}/HWD_national_rcp45_mean_v1.0.nc', decode_coords='all')
hwd85n = xr.open_dataset(f'{data_dir}/HWD_national_rcp85_mean_v1.0.nc', decode_coords='all')
hwd45h.rio.write_crs("epsg:4326", inplace=True)
hwd85h.rio.write_crs("epsg:4326", inplace=True)
hwd45n.rio.write_crs("epsg:4326", inplace=True)
hwd85n.rio.write_crs("epsg:4326", inplace=True)
<xarray.Dataset> Size: 102MB
Dimensions:      (lat: 425, lon: 599, time: 100)
Coordinates:
    height       float64 8B ...
  * lat          (lat) float64 3kB 30.1 30.2 30.3 30.4 ... 72.2 72.3 72.4 72.5
  * lon          (lon) float64 5kB -24.9 -24.8 -24.7 -24.6 ... 34.7 34.8 34.9
  * time         (time) datetime64[ns] 800B 1986-01-01 1987-01-01 ... 2085-01-01
    spatial_ref  int64 8B 0
Data variables:
    HWD_merged   (time, lat, lon) float32 102MB ...
Attributes:
    title:          Processed EURO-CORDEX future climate data for the health ...
    conventions:    CF-1.6
    project:        Copernicus Climate Change Service Sectoral Information Sy...
    source:         Processing of bias-corrected EURO-CORDEX data by VITO
    contact:        bd_rma@vito.be
    creation_date:  Thu May  9 14:08:15 2019
    institution:    VITO (https://vito.be/en)
    history:        Mon Jul 15 13:48:58 2019: ncatted -O -a units,HWD_merged,...
    NCO:            netCDF Operators version 4.7.5 (Homepage = http://nco.sf....
# This code selects the variable for plotting within each dataset
hwd45h=hwd45h['HWD_EU_health']
hwd85h=hwd85h['HWD_EU_health']
hwd45n=hwd45n['HWD_merged']
hwd85n=hwd85n['HWD_merged']
hwd45health= hwd45h.mean(dim='time', skipna=True, keep_attrs=True)
hwd45health.rio.to_raster(raster_path=f'{data_dir}/hwd45_health_mean.tif')
hwd45h.rio.to_raster(raster_path=f'{data_dir}/hwd45_health.tif')
#
hwd85health= hwd85h.mean(dim='time', skipna=True, keep_attrs=True)
hwd85health.rio.to_raster(raster_path=f'{data_dir}/hwd85_health_mean.tif')
hwd85h.rio.to_raster(raster_path=f'{data_dir}/hwd85_health.tif')
#
hwd45national= hwd45n.mean(dim='time', skipna=True, keep_attrs=True)
hwd45national.rio.to_raster(raster_path=f'{data_dir}/hwd45_national_mean.tif')
hwd45n.rio.to_raster(raster_path=f'{data_dir}/hwd45_national.tif')
#
hwd85national= hwd85n.mean(dim='time', skipna=True, keep_attrs=True)
hwd85national.rio.to_raster(raster_path=f'{data_dir}/hwd85_national_mean.tif')
hwd85n.rio.to_raster(raster_path=f'{data_dir}/hwd85_national.tif')

Select the location for plotting#

We will plot a 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}/hwd45_health_mean.tif')
# Create ipyleaflet tile layer from that server
t1 = get_leaflet_tile_layer(client1, cmap='Reds', opacity=0.5, nodata=0, name='Heat-wave health rcp 4.5')
#m = Map(center=client1.center(), zoom=client1.default_zoom)
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(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)
# Arrange the map and text widget vertically
display(widgets.VBox([m, coord_output]))

On the map above you can see the mean values of the heat-wave occurrence for the years 1986-2085 for the RCP4.5 based on the health-related definitions. The total number of heat occurrences for this map is not important, important is to see which areas can be most influenced by the heat.

To select the location for further analysis:

  • Zoom with [+]/[-] to your area,

  • Select the location for plotting the values, with “Draw a marker” from the left panel.

# This code extracts the data values for the selected pixel
hwd45h = xr.open_dataset(f'{data_dir}/hwd45_health.tif')
hwd85h = xr.open_dataset(f'{data_dir}/hwd85_health.tif')
hwd45n = xr.open_dataset(f'{data_dir}/hwd45_national.tif')
hwd85n = xr.open_dataset(f'{data_dir}/hwd85_national.tif')

lat = point_lat_list[0] 
lon = point_lon_list[0]

extracted_data1 = hwd45h.sel(y=lat, x=lon, method='nearest')
extracted_data2 = hwd85h.sel(y=lat, x=lon, method='nearest')
extracted_data3 = hwd45n.sel(y=lat, x=lon, method='nearest')
extracted_data4 = hwd85n.sel(y=lat, x=lon, method='nearest')

d1=extracted_data1['band_data']
hwd45h = d1.values.tolist()
d2=extracted_data2['band_data']
hwd85h = d2.values.tolist()
d3=extracted_data3['band_data']
hwd45n = d3.values.tolist()
d4=extracted_data4['band_data']
hwd85n = d4.values.tolist()

Plot results for heatwaves defined using national thresholds#

The heatwave data calculated using national heatwave definitions is only available for a few states in the EU (Belgium, Hungary, Italy, United Kingdom, Sweden, Lithuania, Latvia, Estonia, Finland).

if not any(np.isnan(hwd45n)): # check if there is national data for the previously selected location
    # Define years
    year = list(range(1986, 2086))
    # Fit polynomial regression curves (degree 2 for simplicity)
    b_poly = np.polyfit(np.arange(len(b)), b, 2)
    t_poly = np.polyfit(np.arange(len(t)), t, 2)
    # Calculate polynomial values
    b_poly_values = np.polyval(b_poly, np.arange(len(b)))
    t_poly_values = np.polyval(t_poly, np.arange(len(t)))
    # Calculate the residuals and standard error for b
    b_residuals = b - b_poly_values
    b_std_error = np.std(b_residuals)
    # Calculate the residuals and standard error for t
    t_residuals = t - t_poly_values
    t_std_error = np.std(t_residuals)
    # Create figure
    fig = go.Figure()
    # Add bar trace for 'b'
    fig.add_trace(go.Bar(x=year, y=b, name='RCP 4.5', marker_color='blue'))
    # Add bar trace for 't'
    fig.add_trace(go.Bar(x=year, y=t, name='RCP 8.5', opacity=0.5, marker_color='red'))
    # Update layout
    fig.update_layout(
        title='Heatwave occurrence in selected pixel by year for the period 1986-2086 rcp4.5 and 8.5 [national-related]',
        xaxis_title='Year',
        yaxis_title='Number of heatwaves',
        legend=dict(x=0, y=1, traceorder='normal', orientation='h'),
        width=1100,
        height=500
    )
    fig.show()
else:
    print("No data is available for this location based on national heatwave definition.")
No data is available for this location based on national heatwave definition.

The graph above shows the heatwave occurrence for the selected location for 1986-2086 for RCP4.5 and RCP8.5. These results are based on the national thresholds for defining a heatwave, based on the EuroHEAT project [source].

Conclusion for the hazard assessment result based on EuroHEAT#

  • The results show the projected heatwave occurrence per year. We can plot the results based on the health-related thresholds or national thresholds for the heatwave definition.

  • The health-related data are available for the whole EU, but the national-related data are available only for a few EU states.

  • The workflow allows to evaluate the projected trend in the heat wave occurence over time under two climate change scenarios at a local resolution of 12 x 12 km. Please note that this data does not include the urban heat-island effect and therefore is not fully representative for urban environments. Extra corrections for urban locations are necessary to include the increased temperatures within built environment.

Authors#

Main contributors: