⚠️ 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: Change under climate scenarios workflow [Hazard assessment]#

Step 1: Description#

Content:

  • 1.1 Sources of the Heat-related products provided by CDS and Climate adapt

  • 1.2 heatwave hazard estimation based on the Xclim methodology:

1.1 Heatwave Hazard methodologies description#

Heatwave based on the XCLIM package for the calculation of the climate indices [XCLIM]#

With this package, you can calculate heat and cold spells, heatwave frequency (years), and total length. Different methodologies for the Heatwave identification. Need to download also the minimum temperature data. The heatwave is estimated by the selected temperature thresholds for minimum and maximum temperature (these thresholds can be estimated based on the user preferences)

Advantages:

  • heatwave calculation based on the minimum and maximum temperature, better for the regions where the temperature can significantly cool down at nights

  • fast processing

  • 12x12km grid for years 1971-2000, 2011-2100n

  • With only one line of code, is easier to understand -the possibility of changing the thresholds (for maximum and minimum temperature)is given by the authorities which take into consideration the influence of heat on human health

  • This methodology is more suitable for the heat-eave estimation in the summer months because the treshold is estimated considering the impact on human health

  • With the XCLIM package you can calculate the tens of climate indices.

Disadvantages

  • We do not see the code behind it - we need to trust to XCLIM package

  • 40 GB of input data

  • heatwave estimation only for the year

Step 2: Prepare a workspace and load a libraries#

2.1 Import packages#

import zipfile        # Working with ZIP archive files (compression and extraction)
import os             # Interacting with the operating system, including handling the current working directory
import rasterio       # Reading and writing geospatial raster data
from pathlib import Path # Working with file system paths in an object-oriented way
import cdsapi         # Interacting with the Climate Data Store (CDS) API for downloading climate data
import numpy as np    # Handling large, multi-dimensional arrays and matrices, along with a collection of mathematical functions to operate on these arrays
import xarray as xr   # Handling labeled multi-dimensional arrays, commonly used for 2D and 3D array data handling
import cartopy.crs as ccrs # Cartographic projections and coordinate transformations for geospatial data
import cartopy.feature as cf # Adding common geographic features to maps
from rasterio.crs import CRS  # Handling Coordinate Reference Systems (CRS) in raster data
from pylab import pcolormesh, show, colorbar, plot, title, legend, subplot, savefig # Various plotting functions, mainly for 2D plotting
from rasterio.plot import show # Displaying raster data
import rasterio as rio # Reading and writing geospatial raster data (redundant with previous import of rasterio)
import plotly.graph_objects as go # Creating interactive plots and visualizations
from ipyleaflet import Map, DrawControl, Marker, LayersControl # Creating interactive maps in Jupyter notebooks
import ipywidgets as widgets # Creating interactive HTML widgets for Jupyter notebooks
import leafmap.leafmap as leafmap # Visualizing geospatial data and interactive mapping
from localtileserver import get_leaflet_tile_layer, TileClient # Serving local raster data as tiles and integrating with Leaflet maps

2.1.1 XCLIM packages#

# This code imports the XCLIM package
import xclim.indices as xci
import xclim

2.2 Create a directory structure#

# Define the directory for the flashflood workflow preprocess
workflow_folder = 'Heat_test2'
# 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')
r1971_2000_dir = os.path.join(workflow_folder,'r1971_2000')
r85_2011_2040_dir = os.path.join(workflow_folder,'r85_2011_2040')
r85_2041_2070_dir = os.path.join(workflow_folder,'r85_2041_2070')
r85_2071_2100_dir = os.path.join(workflow_folder,'r85_2071_2100')
r45_2011_2040_dir = os.path.join(workflow_folder,'r45_2011_2040')
r45_2041_2070_dir = os.path.join(workflow_folder,'r45_2041_2070')
r45_2071_2100_dir = os.path.join(workflow_folder,'r45_2071_2100')
# 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))
    os.makedirs(os.path.join(r1971_2000_dir))
    os.makedirs(os.path.join(r85_2011_2040_dir))
    os.makedirs(os.path.join(r85_2041_2070_dir))
    os.makedirs(os.path.join(r85_2071_2100_dir))
    os.makedirs(os.path.join(r45_2011_2040_dir))
    os.makedirs(os.path.join(r45_2041_2070_dir))
    os.makedirs(os.path.join(r45_2071_2100_dir))

Step 3: Data download and preparation#

3.1 Download a data from the CDS#

URL = "https://cds.climate.copernicus.eu/api/v2"
KEY = "PUT-HERE-YOUR-KEY!!!" ### put here your key!!!
  • download of the maximum 2m air temperature (for PESETA IV), if you want to use also XCLIM download also minimum 2m air temperature

  • This code will automatically download the EuroCordex maximum air temperature data for the selected period and RCP

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

    • More information about the data you can find here: [CDS data source]

    • More information about what is the RCP can find here: [Info about RCP]

    • For PESETA IV download only maximum 2m temperature (minimum 2m temperature download only if you want to use XCLIM package)

1971-2000 maximum 2m temperature in the last 24h

# This code downloads the data for the selected time-period and rcp
## Data for the 1971-2000
c = cdsapi.Client(url=URL, key=KEY)
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'era5_daily_t2m_eurminmax_1971_2000.zip')
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': 'cnrm_cerfacs_cm5',
        'rcm_model': 'clmcom_clm_cclm4_8_17',
        'ensemble_member': 'r1i1p1',
        'start_year': [
            '1971', '1976', '1981', '1986', '1991', '1996', 
        ],
        'end_year': [
            '1975', '1980', '1985', '1990', '1995', '2000', 
        ],
    },
   f"{data_dir}/era5_daily_t2m_eurminmax_1971_2000.zip")
# This code unzips the downloaded files in your working directory, so they will be ready for computing 
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'era5_daily_t2m_eurminmax_1971_2000.zip')
# Extract from zip file
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

1971-2000 minimum 2m temperature in the last 24h (only for XCLIM methodology)

# This code downloads the data for the selected time and rcp
## Data for the 1971-2000
c = cdsapi.Client(url=URL, key=KEY)
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'era5_daily_t2m_eurmin_1971_2000.zip')
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': 'cnrm_cerfacs_cm5',
        'rcm_model': 'clmcom_clm_cclm4_8_17',
        'ensemble_member': 'r1i1p1',
        'start_year': [
            '1971', '1976', '1981', '1986', '1991', '1996', 
        ],
        'end_year': [
            '1975', '1980', '1985', '1990', '1995', '2000', 
        ],
    },
   f"{data_dir}/era5_daily_t2m_eurmin_1971_2000.zip")
# This code unzips the downloaded files in your working directory, so they will be ready for computing 
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'era5_daily_t2m_eurmin_1971_2000.zip')
# Extract from zip file
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

2011-2100 rcp_8_5 maximum 2m temperature in the last 24h (or you can select other suggested period 2041-2070, 2071-2100 or rcp 4.5 or 8.5)

# This code  downloads the data for the selected time period and rcp
c = cdsapi.Client(url=URL, key=KEY)
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'era5_daily_t2m_eurminmax_2011_2100.zip')
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': 'cnrm_cerfacs_cm5',
        'rcm_model': 'clmcom_clm_cclm4_8_17',
        'ensemble_member': 'r1i1p1',
        'start_year': [
            '2011', '2016', '2021', '2026', '2031', '2036', '2041', '2046', '2051', '2056', '2061', '2066', '2071', '2076', '2081', '2086', '2091', '2096',
        ],
        'end_year': [
            '2015', '2020', '2025', '2030', '2035', '2040', '2045', '2050', '2055', '2060', '2065', '2070', '2075', '2080', '2085', '2090', '2095', '2100',
        ],
    },
   f"{data_dir}/era5_daily_t2m_eurminmax_2011_2100.zip")
# This code unzips the downloaded files in your working directory, data will be ready for computing 
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'era5_daily_t2m_eurminmax_2011_2100.zip')
# Extract from zip file
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

2011-2100 rcp_8_5 minimum 2m temperature in the last 24h (only for XCLIM methodology) (or you can select other suggested period 2041-2070, 2071-2100 or rcp 4.5 or 8.5 in ‘experiment’: rcp_8_5 to rcp_4_5)

# This code downloads the data for the selected time period and rcp
c = cdsapi.Client(url=URL, key=KEY)
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'era5_daily_t2m_eurmin_2011_2100.zip')
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': 'cnrm_cerfacs_cm5',
        'rcm_model': 'clmcom_clm_cclm4_8_17',
        'ensemble_member': 'r1i1p1',
        'start_year': [
            '2011', '2016', '2021', '2026', '2031', '2036', '2041', '2046', '2051', '2056', '2061', '2066', '2071', '2076', '2081', '2086', '2091', '2096',
        ],
        'end_year': [
            '2015', '2020', '2025', '2030', '2035', '2040', '2045', '2050', '2055', '2060', '2065', '2070', '2075', '2080', '2085', '2090', '2095', '2100',
        ],
    },
   f"{data_dir}/era5_daily_t2m_eurmin_2011_2100.zip")
# This code unzips the downloaded files in your working directory, so they will be ready for computing 
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'era5_daily_t2m_eurmin_2011_2100.zip')
# Extract from zip file
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

2011-2100 rcp_4_5 maximum 2m temperature in the last 24h

# This code will download the data for the selected time-period and rcp
c = cdsapi.Client(url=URL, key=KEY)
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'era5_daily_t2m_eurmax_rcp45_2011_2100.zip')
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': 'cnrm_cerfacs_cm5',
        'rcm_model': 'clmcom_clm_cclm4_8_17',
        'ensemble_member': 'r1i1p1',
        'start_year': [
            '2011', '2016', '2021', '2026', '2031', '2036', '2041', '2046', '2051', '2056', '2061', '2066', '2071', '2076', '2081', '2086', '2091', '2096',
        ],
        'end_year': [
            '2015', '2020', '2025', '2030', '2035', '2040', '2045', '2050', '2055', '2060', '2065', '2070', '2075', '2080', '2085', '2090', '2095', '2100',
        ],
    },
   f"{data_dir}/era5_daily_t2m_eurmax_rcp45_2011_2100.zip")
# This code unzips the downloaded files in your working directory, so they will be ready for computing 
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'era5_daily_t2m_eurmax_rcp45_2011_2100.zip')
# Extract from zip file
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

2011-2100 rcp_4_5 minimum 2m temperature in the last 24h (only for XCLIM methodology)

# This code downloads the data for the selected time period and rcp
c = cdsapi.Client(url=URL, key=KEY)
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'era5_daily_t2m_eurmin_rcp45_2011_2100.zip')
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': 'cnrm_cerfacs_cm5',
        'rcm_model': 'clmcom_clm_cclm4_8_17',
        'ensemble_member': 'r1i1p1',
        'start_year': [
            '2011', '2016', '2021', '2026', '2031', '2036', '2041', '2046', '2051', '2056', '2061', '2066', '2071', '2076', '2081', '2086', '2091', '2096',
        ],
        'end_year': [
            '2015', '2020', '2025', '2030', '2035', '2040', '2045', '2050', '2055', '2060', '2065', '2070', '2075', '2080', '2085', '2090', '2095', '2100',
        ],
    },
   f"{data_dir}/era5_daily_t2m_eurmin_rcp45_2011_2100.zip")
# This code will unzip the downloaded files in your working directory, so they will be ready for computing 
# Define zip file's absolute path
zip_path = os.path.join(data_dir, 'era5_daily_t2m_eurmin_rcp45_2011_2100.zip')
# Extract from zip file
with zipfile.ZipFile(zip_path, 'r') as zObject:
    zObject.extractall(path=data_dir)

3.2 Load data and convert Kelvin to Celsius#

# This code loads a data from your working directory, and set the coordinates system (CRS)
# Historical data Maximum temperature
dmaxh = xr.open_mfdataset(f'{data_dir}/tasmax*historical*.nc', decode_coords='all')
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*historical*.nc', decode_coords='all')
rotated_pole = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
dminh.rio.write_crs(rotated_pole, inplace=True)
# Maximum temperature rcp 4.5
dmax45 = xr.open_mfdataset(f'{data_dir}/tasmax*rcp45*.nc', decode_coords='all')
rotated_pole = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
dmax45.rio.write_crs(rotated_pole, inplace=True)
# Minimum temperature rcp 4.5
dmin45 = xr.open_mfdataset(f'{data_dir}/tasmin*rcp45*.nc', decode_coords='all')
rotated_pole = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
dmin45.rio.write_crs(rotated_pole, inplace=True)
# Maximum temperature rcp 8.5
dmax85 = xr.open_mfdataset(f'{data_dir}/tasmax*rcp85*.nc', decode_coords='all')
rotated_pole = ccrs.RotatedPole(pole_latitude=39.25, pole_longitude=-162)
dmax85.rio.write_crs(rotated_pole, inplace=True)
# Minimum temperature rcp 8.5
dmin85 = xr.open_mfdataset(f'{data_dir}/tasmin*rcp85*.nc', decode_coords='all')
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!!!

3.3 Selection of the bbox (boundary box)#

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 EuroCordex 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(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)
# Arrange the map and text widget vertically
display(widgets.VBox([m, coord_output]))

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: [[18.199646, 48.917587], [19.40311, 49.566692]]
Rotated Pole Coordinates: [(0.13126311280544092, -1.832235902681797), (0.9101431184540201, -1.1746779340800233)]
# 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)
(18.199646, 48.917587, 19.40311, 49.566692)
# 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)
(0.13126311280544092, -1.832235902681797, 0.9101431184540201, -1.1746779340800233)

3.4 Clip area 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])  

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

# This code will convert the maximum temperature data from Kelvin to Celsius. The first part will also remove unnecessary data from the disk
# Remove data from the memory
del dmaxh
del dminh
del dmax45
del dmin45
del dmax85
del dmin85
# 2. Convert from K to °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 an extensible framework for constructing custom climate indicators, statistical downscaling and bias adjustment of climate model simulations, as well as climate model ensemble analysis tools. XCLIMim 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 services data treatment workflows. 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)

The thresholds of 22° and 25°C for night temperatures and 30° and 35°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].

In Robinson [2001], the parameters would be: thresh_tasmin=27.22, thresh_tasmax=39.44, window=2 (81F, 103F). Casati, Yagouti, and Chaumont [2013], Robinson [2001]

  • The biggest difference vs PESETA IV is that in XCLIM is the treshold for the heatwaves calculated directly from the 30 years (rough estimate, identifies heatwave mainly in summer months), while in PESETA IV is calculated for monthly values from 30 years and we can detect also seasonal heatwave changes.

4.1 Data load#

  • For XCLIM we are using the same data as for PESETA IV + data about minimum air temperature

  • These data are downloaded at the beginning of this notebook

  • We can select our own bbox or use the selected one at the begginig of the notebook

  • Before proceeding this code you need to run whole section 1

# 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", "2000-12-31"))
tasmaxp_45=tasmax45.sel(time=slice("2011-01-01", "2100-12-31"))
tasmaxp_85=tasmax85.sel(time=slice("2011-01-01", "2100-12-31"))
# Min T
tasminh=tasminh.sel(time=slice("1971-01-01", "2000-12-31"))
tasminp_45=tasmin45.sel(time=slice("2011-01-01", "2100-12-31"))
tasminp_85=tasmin85.sel(time=slice("2011-01-01", "2100-12-31"))

4.2 Defining of the temperature and duration tresholds#

Treshold temperatures based on the percentile:

  1. Compute the percentile for your selected area

  2. Calculate the mean value for your selected area

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]

Temperature treshold for maximum daily air temperature (in degC) = thresh

Duration threshold (in days) = window

# 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='25.0 degC', window=5, freq='YS', op='>', resample_before_rl=True)
# RCP 8.5
HWIp85=xclim.indices.heat_wave_index(tasmaxp_85, thresh='25.0 degC', window=5, freq='YS', op='>', resample_before_rl=True)
# Historical
HWIh=xclim.indices.heat_wave_index(tasmaxh, thresh='25.0 degC', window=5, freq='YS', op='>', resample_before_rl=True)

Temperature treshold for minimum daily air temperature (in degC) = thresh_tasmin

Temperature treshold for maximum daily air temperature (in degC) = thresh_tasmax

Duration threshold (in days) = window

# heatwave frequency
# RCP 4.5
HWFp45=xclim.indices.heat_wave_frequency(tasminp_45, tasmaxp_45, thresh_tasmin='17.0 degC', thresh_tasmax='28.0 degC', window=3, freq='YS', op='>', resample_before_rl=True)
# RCP 8.5
HWFp85=xclim.indices.heat_wave_frequency(tasminp_85, tasmaxp_85, thresh_tasmin='17.0 degC', thresh_tasmax='28.0 degC', window=3, freq='YS', op='>', resample_before_rl=True)
#Historical
HWFh=xclim.indices.heat_wave_frequency(tasminh, tasmaxh, thresh_tasmin='17.0 degC', thresh_tasmax='28.0 degC', window=3, freq='YS', op='>', resample_before_rl=True)
# Data correction
# Set a right type
HWFp45=HWFp45.astype(np.float64)
HWFp85=HWFp85.astype(np.float64)
HWFh=HWFh.astype(np.float64)
  • Heatwave total length [source]

Temperature treshold for minimum daily air temperature (in degC) = thresh_tasmin

Temperature treshold for minimum daily air temperature (in degC) = thresh_tasmax

Duration threshold (in days) = window

# heatwave total length
# RCP 4.5
HWTLp45=xclim.indices.heat_wave_total_length(tasminp_45, tasmaxp_45, thresh_tasmin='17.0 degC', thresh_tasmax='28.0 degC', window=3, freq='YS', op='>', resample_before_rl=True)
# RCP 8.5
HWTLp85=xclim.indices.heat_wave_total_length(tasminp_85, tasmaxp_85, thresh_tasmin='17.0 degC', thresh_tasmax='28.0 degC', window=3, freq='YS', op='>', resample_before_rl=True)
# Historical
HWTLh=xclim.indices.heat_wave_total_length(tasminh, tasmaxh, thresh_tasmin='17.0 degC', thresh_tasmax='28.0 degC', window=3, 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 writed!!!")
  • This step will take a few minutes!!!

# This code saves the computed indices in the disk in raster tif format 
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')
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')
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')

4.3 Select the pixel 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.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')
# 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, cmap='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(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]))
  • 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

  • These data are scaled for heatwave Index

  • This graphical comparison is only for the better imagination of which place 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)

4.4 Plots 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.

# 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, 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='HWI RCP 4.5', marker_color='blue'))
# Add bar trace for 't'
fig.add_trace(go.Bar(x=year, y=t, name='HWI RCP 8.5', opacity=0.5, marker_color='red'))
# Update layout
fig.update_layout(
    title='Number of days that are part of a heatwave defined as five or more consecutive days over a threshold of 25 degrees',
    xaxis_title='Year',
    yaxis_title='Number of heatwaves',
    legend=dict(x=0, y=1, traceorder='normal', orientation='h'),
    width=1100,
    height=500)