Heat Index Indicators
Introduction
Heatwaves have become one significant weather-related hazards, affecting human health, ecosystems, and various economic sectors.
To better quantify the magnitude and the effects, several indices have been developed to represent thermal comfort and heat stress. Among them, the Heat Index computed by the U.S. National Weather Service (NWS) will be used to implement this indicator.
Heat Index computes the perceived temperature resulting from the combined effects of air temperature and relative humidity. It is used to indicate discomfort or potential health risks due to heat. The computation is based on a formula mainly developed by Lans P. Rothfusz, which estimates the apparent temperature under warm and humid conditions. The index increases with higher humidity levels because the body’s ability to cool itself through evaporation (sweating) becomes less effective.
The computed indicator is:
- HeatIndex
The Heat Index can be computed using both multidimensional arrays and s2dv_cube objects. Functions prefixed with CST_ work with s2dv_cube objects, while the non-prefixed versions operate on multidimensional arrays.
Some supplementary functions are useful to compute the Heat Index and allow additional operations such as basic aggregation (mean, max, min), quantile calculation, or application of other CSIndicators functions.
Here, two different functions are used to illustrate the computation procedure:
-
CST_HeatIndex (Mean) computes the mean Heat Index over a period, here for August, representing typical thermal conditions.
-
CST_HeatIndex (TotalSpellTimeExceedingThreshold) counts consecutive days where the Heat Index exceeds 26.7°C (80°F), corresponding to the first level of health risk where fatigue is possible with prolonged exposure or physical activity.
When period selection is required, start and end parameters must be provided to subset the time dimension. The Heat Index indicator uses helper functions to support period and data selection when needed. Otherwise, the entire period is used.
Data Aquisition
Load required libraries.
library(CSIndicators)
library(CSTools)
library(s2dv)
library(multiApply)
source("https://gitlab.earth.bsc.es/es/csindicators/-/raw/master/R/HeatIndex.R?ref_type=heads")
Define spatial domain (West Europe) and forecast start dates (August 1st from 1993 to 2016).
lat_min = 30
lat_max = 60
lon_min = -10
lon_max = 20
sdates <- paste0(1993:2016, '08', '01')
To compute the Heat Index, we load daily mean temperature (tas) and relative humidity (hurs) from ECMWF system51c3s forecasts for August 1993–2016 and from ERA5 reanalysis for the same period. ERA5 data are regridded onto the forecast grid using conservative remapping to ensure consistency.
PATH_exp_system51 <- paste0("/esarchive/exp/ecmwf/system51c3s/daily_mean/$var$_f6h/$var$_$sdate$.nc")
tas_exp <- CST_Start(dataset = PATH_exp_system51,
var = "tas",
member = startR::indices(1:3),
sdate = sdates,
ftime = startR::indices(1:31),
lat = startR::values(list(lat_min, lat_max)),
lat_reorder = startR::Sort(decreasing = TRUE),
lon = startR::values(list(lon_min, lon_max)),
lon_reorder = startR::CircularSort(0, 360),
synonims = list(lon = c('lon', 'longitude'),
lat = c('lat', 'latitude'),
member = c('member', 'ensemble'),
ftime = c('ftime', 'time')),
transform_vars = c('lat', 'lon'),
return_vars = list(lat = NULL,
lon = NULL, ftime = 'sdate'),
retrieve = TRUE)
hurs_exp <- CST_Start(dataset = PATH_exp_system51,
var = "hurs",
member = startR::indices(1:3),
sdate = sdates,
ftime = startR::indices(1:31),
lat = startR::values(list(lat_min, lat_max)),
lat_reorder = startR::Sort(decreasing = TRUE),
lon = startR::values(list(lon_min, lon_max)),
lon_reorder = startR::CircularSort(0, 360),
synonims = list(lon = c('lon', 'longitude'),
lat = c('lat', 'latitude'),
member = c('member', 'ensemble'),
ftime = c('ftime', 'time')),
transform_vars = c('lat', 'lon'),
return_vars = list(lat = NULL,
lon = NULL, ftime = 'sdate'),
retrieve = TRUE)
dates_exp <- tas_exp$attrs$Dates
PATH_obs_era5 <- paste0("/esarchive/recon/ecmwf/era5/daily_mean/$var$_f1h-r1440x721cds/$var$_$date$.nc")
tas_obs <- CST_Start(dataset = PATH_obs_era5,
var = "tas",
date = unique(format(dates_exp, '%Y%m')),
ftime = startR::values(dates_exp),
ftime_across = 'date',
ftime_var = 'ftime',
merge_across_dims = TRUE,
split_multiselected_dims = TRUE,
lat = startR::values(list(lat_min, lat_max)),
lat_reorder = startR::Sort(decreasing = TRUE),
lon = startR::values(list(lon_min, lon_max)),
lon_reorder = startR::CircularSort(0, 360),
synonims = list(lon = c('lon', 'longitude'),
lat = c('lat', 'latitude'),
ftime = c('ftime', 'time')),
transform = startR::CDORemapper,
transform_extra_cells = 2,
transform_params = list(grid = "/esarchive/scratch/abojaly/GitLab/sunset/conf/grid_description/griddes_system51c3s.txt",
method = "conservative"),
transform_vars = c('lat', 'lon'),
return_vars = list(lon = NULL,
lat = NULL,
ftime = 'date'),
retrieve = TRUE)
hurs_obs <- CST_Start(dataset = PATH_obs_era5,
var = "hurs",
date = unique(format(dates_exp, '%Y%m')),
ftime = startR::values(dates_exp),
ftime_across = 'date',
ftime_var = 'ftime',
merge_across_dims = TRUE,
split_multiselected_dims = TRUE,
lat = startR::values(list(lat_min, lat_max)),
lat_reorder = startR::Sort(decreasing = TRUE),
lon = startR::values(list(lon_min, lon_max)),
lon_reorder = startR::CircularSort(0, 360),
synonims = list(lon = c('lon', 'longitude'),
lat = c('lat', 'latitude'),
ftime = c('ftime', 'time')),
transform = startR::CDORemapper,
transform_extra_cells = 2,
transform_params = list(grid = "/esarchive/scratch/abojaly/GitLab/sunset/conf/grid_description/griddes_system51c3s.txt",
method = "conservative"),
transform_vars = c('lat', 'lon'),
return_vars = list(lon = NULL,
lat = NULL,
ftime = 'date'),
retrieve = TRUE)
The output contains data and metadata for both the forecast and observational datasets. This shows the order and size of dimensions:
The summary of the data values to inspect ranges, units, and possible anomalies. Useful for checking for missing values, extreme values, or unexpected distributions.
dim(tas_exp$data)
# dataset var member sdate ftime lat lon
# 1 1 3 24 31 30 30
dim(tas_obs$data)
# dataset var sdate ftime lat lon
# 1 1 24 31 30 30
dim(hurs_exp$data)
# dataset var member sdate ftime lat lon
# 1 1 3 24 31 30 30
dim(hurs_obs$data)
# dataset var sdate ftime lat lon
# 1 1 24 31 30 30
summary(tas_exp$data)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 273.2 289.5 293.5 294.2 298.8 313.9
summary(tas_obs$data)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 273.5 289.7 293.7 294.4 298.7 314.8
summary(hurs_exp$data)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 6.14 62.17 74.37 68.38 81.41 99.76
summary(hurs_obs$data)
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 6.833 64.328 75.221 69.893 82.109 99.338
Quick temperature overview
This section provides a quick overview of the climatology, showing the spatial distribution of temperature and relative humidity from ERA5 observations. The long-term mean values for August 1993 - 2016 serve as a reference for interpreting the Heat Index and identifying areas prone to heat stress.
map_tas <- apply(tas_obs$data, c(5,6), mean, na.rm = TRUE) - 273.15
map_hurs <- apply(hurs_obs$data, c(5,6), mean, na.rm = TRUE)
brks <- seq(ceiling(min(map_tas)), ceiling(max(map_tas)), by = 1)
map_plot <- t(map_tas)
names(dim(map_plot)) <- c('longitude','latitude')
PlotEquiMap(map_plot,
lat=tas_obs$coords$lat,
lon=tas_obs$coords$lon,
filled.continents = F,
country.borders = TRUE,
triangle_ends = c(TRUE, TRUE),
toptitle = paste("ERA5 / Mean Temperature in August / Reference 1993 - 2016"),
title_scale = 0.6,
units = "deg C",
brks = brks)
brks <- seq(ceiling(min(map_hurs)), ceiling(max(map_hurs)), length.out = 50)
cols <- colorRampPalette(c("#8c510a", "#d8b365", "#f6e8c3", "#c7eae5", "#016c9a"))(length(brks)-1)
map_plot <- t(map_hurs)
names(dim(map_plot)) <- c('longitude','latitude')
PlotEquiMap(map_plot,
lat=hurs_obs$coords$lat,
lon=hurs_obs$coords$lon,
filled.continents = F,
country.borders = TRUE,
triangle_ends = c(TRUE, TRUE),
toptitle = paste("ERA5 / Mean Humidity in August / Reference 1993 - 2016"),
title_scale = 0.6,
cols = cols,
col_inf = cols[1],
col_sup = cols[length(cols)],
units = "%",
brks = brks)
Heat Index computation
This section presents the computation of the Heat Index from daily temperature and relative humidity for both forecast (system51c3s) and observational (ERA5) datasets. The mean Heat Index represents thermal conditions in August, while periods exceeding a critical threshold of 26.7°C (80°F), corresponding to the first level of health risk where fatigue is possible with prolonged exposure or physical activity, highlight heat stress events and areas at risk.
# EXP
HI_exp_mean <- CST_HeatIndex(temp = tas_exp, rh = hurs_exp, temp_units_output = "C", time_dim= "ftime")
HI_exp_thresh <- CST_HeatIndex(temp = tas_exp, rh = hurs_exp, temp_units_output = "C", time_dim = "ftime",
fun = function(x, na.rm = TRUE) {data <- array(x, dim = c(ftime = length(x)))
TotalSpellTimeExceedingThreshold(data = data, threshold = 26.7, spell = 3,
op = ">", time_dim = "ftime")})
# OBS
HI_obs_mean <- CST_HeatIndex(temp = tas_obs, rh = hurs_obs, temp_units_output = "C", time_dim= "ftime")
HI_obs_thresh <- CST_HeatIndex(temp = tas_obs, rh = hurs_obs, temp_units_output = "C",
time_dim = "ftime",
fun = function(x, na.rm = TRUE) {data <- array(x, dim = c(time = length(x)))
TotalSpellTimeExceedingThreshold(data = data, threshold = 26.7, spell = 3,
op = ">", time_dim = "time")})
Visualisation
This section presents visualizations of Heat Index anomalies and forecast bias. Anomalies for August 2003 are calculated relative to the mean of the reference period 1993 - 2016 (excluding 2003), highlighting the heatwave over France during this period. Forecast bias is computed as the difference between observed and forecasted mean Heat Index, indicating regions where the forecast overestimates or underestimates thermal stress.
# Anomalies
map_obs2003 <- t(HI_obs_mean$data[,,11,,,])
names(dim(map_obs2003)) <- c("longitude", "latitude")
map_obs <- t(apply(HI_obs_mean$data[,,-11,,,], c(2,3), mean, na.rm = TRUE))
names(dim(map_obs)) <- c("longitude", "latitude")
max_abs <- ceiling(max(abs(map_obs2003 - map_obs), na.rm = TRUE))
brks <- seq(-max_abs, max_abs, length.out = 21)
PlotEquiMap(map_obs2003 - map_obs,
lat = HI_obs_mean$coords$lat,
lon = HI_obs_mean$coords$lon,
filled.continents = FALSE, country.borders = TRUE,
toptitle = "ERA5 / Heat Index / Anomalies / mean August - 2003 \n Reference Years: 1993 - 2016 (excl. 2003)",
title_scale = 0.5,
brks = brks,
bar_label_scale = 0.8,
bar_extra_margin = c(2,0,1,0),
triangle_ends = c(F, F),
draw_bar_ticks = TRUE,
units = "HI [degree C]")
# Bias
map_obs_mean <- t(apply(HI_obs_mean$data[,,,,,], c(2,3), mean, na.rm = TRUE))
names(dim(map_obs_mean)) <- c("longitude", "latitude")
map_exp_mean <- t(apply(HI_exp_mean$data[,,,,,,], c(3,4), mean, na.rm = TRUE))
names(dim(map_exp_mean)) <- c("longitude", "latitude")
max_abs <- ceiling(max(abs(map_obs_mean - map_exp_mean), na.rm = TRUE))
brks <- seq(-max_abs, max_abs, length.out = 21)
PlotEquiMap(
lat = HI_exp_mean$coords$lat,
lon = HI_exp_mean$coords$lon,
var = map_obs_mean - map_exp_mean,
filled.continents = FALSE,
country.borders = TRUE,
toptitle = "Heat Index / BIAS (obs - hcst) / mean August, 1993 - 2016",
title_scale = 0.6,
brks = brks,
bar_label_scale = 0.8,
triangle_ends = c(FALSE, FALSE),
draw_bar_ticks = TRUE,
units = "HI [°C]")
This section visualizes the mean duration of Heat Index spells exceeding the critical threshold (26.7°C) for both observed and forecasted datasets over the reference period 1993 - 2016. These maps show regions where prolonged heat stress events are most frequent and provide a spatial overview of areas at higher risk for fatigue or heat-related health impacts.
map_thresholds <- MeanDims(data = HI_obs_thresh$data, dims = 'sdate', drop = TRUE)
map_thresholds <- drop(map_thresholds)
brks = seq(ceiling(min(HI_obs_thresh$data)), ceiling(max(HI_obs_thresh$data)), by = 3)
PlotEquiMap(map_thresholds,
lat = HI_obs_mean$coords$lat,
lon = HI_obs_mean$coords$lon,
filled.continents = FALSE, country.borders = TRUE,
toptitle = "ERA5 / Heat Index / Mean Spell Time Exceeding Threshold / Reference 1993 - 2016",
title_scale = 0.5,
bar_label_scale = 0.8,
brks = brks,
units = "Time [days]")
map_thresholds <- MeanDims(data = HI_exp_thresh$data, dims = c('member','sdate'), drop = TRUE)
map_thresholds <- drop(map_thresholds)
brks = seq(ceiling(min(HI_obs_thresh$data)), ceiling(max(HI_obs_thresh$data)), by = 3)
PlotEquiMap(map_thresholds,
lat = HI_obs_mean$coords$lat,
lon = HI_obs_mean$coords$lon,
filled.continents = FALSE, country.borders = TRUE,
toptitle = "SYSTEM5.1 / Heat Index / Mean Spell Time Exceeding Threshold / Reference 1993 - 2016",
title_scale = 0.5,
bar_label_scale = 0.8,
brks = brks,
units = "Time [days]")
