Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Annual max summer heat index

Annual Maximum JJA Heat Index at Boulder, CO and Chicago, IL (1990–2005)

ERA5 and JRA-3Q Reanalyses from NCAR GDEX

Scientific motivation

This notebook computes the annual maximum June–July–August (JJA) heat index at two mid-latitude United States cities — Boulder, Colorado (semi-arid, elevation 1655 m) and Chicago, Illinois (humid continental, Lake Michigan influence) — using two independent global reanalyses spanning 1990–2005.

The analysis follows the methodology of Romps (2024), who demonstrated using ERA5 that the annual maximum summer heat index in Texas has increased at a rate several times larger than the contemporaneous increase in dry-bulb temperature. The physical basis for this amplification lies in the nonlinearity of the heat index with respect to temperature and humidity: at near-saturated conditions, the Clausius–Clapeyron relation implies that a given increment of warming produces a disproportionately large increase in physiological heat stress. Chicago, situated in a persistently humid air mass, is expected to exhibit stronger amplification than semi-arid Boulder, where the heat index remains close to the dry-bulb temperature throughout the summer.

Datasets

ProductGDEX identifierTemporal resolutionQuantity archivedAccess format
ERA5 surface analysisd633000Hourly instantaneousVAR_2T, VAR_2DNative Zarr
JRA-3Q surface analysis (anl_surf125)d6400006-hourly instantaneousTMP_GDS0_HTGL, RH_GDS0_HTGLkerchunk reference

Because ERA5 provides 24 candidate values per day versus 4 for JRA-3Q, computing the JJA season maximum from the full hourly ERA5 stream would systematically inflate ERA5 relative to JRA-3Q through a sample-size effect alone. To remove this artifact, ERA5 is subsampled to one randomly selected observation per 6-hour window prior to computing the annual maximum, matching the JRA-3Q cadence. A random draw, rather than a fixed synoptic hour, is used to avoid introducing a systematic diurnal-phase bias into the sampled ERA5 distribution. We re-sample instead of taking a mean over a 6-hour window because the average might miss an extreme heat index value in this window.

Heat index formulation

The heat index is computed using the heatindex Python package of Lu & Romps (2025), which implements the extended formulation of Lu & Romps (2022) in python. The function takes air temperature in Kelvin and relative humidity on the interval [0, 1][0,\,1] and returns the heat index in Kelvin:

hi.heatindex(T_K, rh)  →  HI_K

ERA5 does not archive near-surface relative humidity directly. It is derived from the 2-m temperature (VAR_2T) and 2-m dew-point temperature (VAR_2D) via the August–Roche–Magnus approximation, consistent with the ECMWF IFS documentation. JRA-3Q provides relative humidity in percent as RH_GDS0_HTGL in the surface analysis product.

Required Packages

Please make sure to install the packages before moving forward

  • intake

  • intake-esm >= 2025.12.12

  • xarray

  • dask

  • zarr

  • kerchunk

  • numpy

  • pandas

  • matplotlib

  • heatindex >= 0.0.2

!pip install heatindex
Defaulting to user installation because normal site-packages is not writeable
Requirement already satisfied: heatindex in /glade/u/home/harshah/.local/lib/python3.10/site-packages (0.0.2)
Requirement already satisfied: numpy in /glade/u/home/harshah/.local/lib/python3.10/site-packages (from heatindex) (1.26.4)

[notice] A new release of pip is available: 25.3 -> 26.0.1
[notice] To update, run: pip install --upgrade pip

import matplotlib.pyplot as plt
import numpy as np
import os
import xarray as xr
import intake
import intake_esm
import pandas as pd
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import dask 
from dask_jobqueue import PBSCluster
from dask.distributed import Client
/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/pyproj/network.py:59: UserWarning: pyproj unable to set PROJ database path.
  _set_context_ca_bundle_path(ca_bundle_path)

Locate the Dataset

On the NCAR GDEX portal, go to the Data Access tab for the ERA5/ JRA-3Q dataset to find the intake-ESM catalogs needed to access data. In this notebook we will use GDEX OSDF catalog.

# Set up your scratch folder path - Ignore these lines if you are not on NCAR's cluster
username       = os.environ["USER"]
glade_scratch  = "/glade/derecho/scratch/" + username
print(glade_scratch)
/glade/derecho/scratch/harshah

# OSDF based data access 
era5_catalog  = 'https://data.gdex.ucar.edu/d633000/catalogs/d633000-osdf.json'
jra3q_catalog = 'https://data.gdex.ucar.edu/d640000/catalogs/d640000-osdf.json'

Step 2 - Set up cluster

  • Please set these two flags to False if you are not on NCAR’s HPC cluster or not using a dask gateway.

  • Setting these flags to False immediately selects a local cluster which can run on your personal device

USE_PBS_SCHEDULER = True # Use NCAR's HPC cluster 
USE_DASK_GATEWAY  = False 
# Create a PBS cluster object
def get_pbs_cluster():
    """ Create cluster through dask_jobqueue.   
    """
    from dask_jobqueue import PBSCluster
    cluster = PBSCluster(
        job_name = 'dask-osdf-24',
        cores = 1,
        memory = '8GiB',
        processes = 1,
        local_directory = glade_scratch + '/dask/spill',
        log_directory = glade_scratch + '/dask/logs/',
        resource_spec = 'select=1:ncpus=1:mem=8GB',
        queue = 'casper',
        walltime = '3:00:00',
        #interface = 'ib0'
        interface = 'ext'
    )
    return cluster

def get_gateway_cluster():
    """ Create cluster through dask_gateway
    """
    from dask_gateway import Gateway

    gateway = Gateway()
    cluster = gateway.new_cluster()
    cluster.adapt(minimum=2, maximum=4)
    return cluster

def get_local_cluster():
    """ Create cluster using the Jupyter server's resources
    """
    from distributed import LocalCluster, performance_report
    cluster = LocalCluster()    

    cluster.scale(5)
    return cluster
# Obtain dask cluster in one of three ways
if USE_PBS_SCHEDULER:
    cluster = get_pbs_cluster()
elif USE_DASK_GATEWAY:
    cluster = get_gateway_cluster()
else:
    cluster = get_local_cluster()

# Connect to cluster
from distributed import Client
client = Client(cluster)
BokehUserWarning: reference already known 'p1011'
BokehUserWarning: reference already known 'p1051'
BokehUserWarning: reference already known 'p1289'
BokehUserWarning: reference already known 'p1295'
BokehUserWarning: reference already known 'p1353'
2026-04-23 18:27:46,037 - tornado.application - ERROR - Uncaught exception GET /status/ws (127.0.0.1)
HTTPServerRequest(protocol='http', host='jupyterhub.hpc.ucar.edu', method='GET', uri='/status/ws', version='HTTP/1.1', remote_ip='127.0.0.1')
Traceback (most recent call last):
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/websocket.py", line 965, in _accept_connection
    open_result = handler.open(*handler.open_args, **handler.open_kwargs)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/web.py", line 3388, in wrapper
    return method(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/bokeh/server/views/ws.py", line 149, in open
    raise ProtocolError("Token is expired. Configure the app with a larger value for --session-token-expiration if necessary")
bokeh.protocol.exceptions.ProtocolError: Token is expired. Configure the app with a larger value for --session-token-expiration if necessary
2026-04-23 18:27:50,727 - tornado.application - ERROR - Uncaught exception GET /status/ws (127.0.0.1)
HTTPServerRequest(protocol='http', host='jupyterhub.hpc.ucar.edu', method='GET', uri='/status/ws', version='HTTP/1.1', remote_ip='127.0.0.1')
Traceback (most recent call last):
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/websocket.py", line 965, in _accept_connection
    open_result = handler.open(*handler.open_args, **handler.open_kwargs)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/web.py", line 3388, in wrapper
    return method(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/bokeh/server/views/ws.py", line 149, in open
    raise ProtocolError("Token is expired. Configure the app with a larger value for --session-token-expiration if necessary")
bokeh.protocol.exceptions.ProtocolError: Token is expired. Configure the app with a larger value for --session-token-expiration if necessary
2026-04-23 18:39:05,909 - tornado.application - ERROR - Uncaught exception GET /status/ws (127.0.0.1)
HTTPServerRequest(protocol='http', host='jupyterhub.hpc.ucar.edu', method='GET', uri='/status/ws', version='HTTP/1.1', remote_ip='127.0.0.1')
Traceback (most recent call last):
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/websocket.py", line 965, in _accept_connection
    open_result = handler.open(*handler.open_args, **handler.open_kwargs)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/web.py", line 3388, in wrapper
    return method(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/bokeh/server/views/ws.py", line 149, in open
    raise ProtocolError("Token is expired. Configure the app with a larger value for --session-token-expiration if necessary")
bokeh.protocol.exceptions.ProtocolError: Token is expired. Configure the app with a larger value for --session-token-expiration if necessary
2026-04-23 18:39:14,910 - tornado.application - ERROR - Uncaught exception GET /status/ws (127.0.0.1)
HTTPServerRequest(protocol='http', host='jupyterhub.hpc.ucar.edu', method='GET', uri='/status/ws', version='HTTP/1.1', remote_ip='127.0.0.1')
Traceback (most recent call last):
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/websocket.py", line 965, in _accept_connection
    open_result = handler.open(*handler.open_args, **handler.open_kwargs)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/web.py", line 3388, in wrapper
    return method(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/bokeh/server/views/ws.py", line 149, in open
    raise ProtocolError("Token is expired. Configure the app with a larger value for --session-token-expiration if necessary")
bokeh.protocol.exceptions.ProtocolError: Token is expired. Configure the app with a larger value for --session-token-expiration if necessary
2026-04-23 18:39:31,643 - tornado.application - ERROR - Uncaught exception GET /status/ws (127.0.0.1)
HTTPServerRequest(protocol='http', host='jupyterhub.hpc.ucar.edu', method='GET', uri='/status/ws', version='HTTP/1.1', remote_ip='127.0.0.1')
Traceback (most recent call last):
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/websocket.py", line 965, in _accept_connection
    open_result = handler.open(*handler.open_args, **handler.open_kwargs)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/web.py", line 3388, in wrapper
    return method(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/bokeh/server/views/ws.py", line 149, in open
    raise ProtocolError("Token is expired. Configure the app with a larger value for --session-token-expiration if necessary")
bokeh.protocol.exceptions.ProtocolError: Token is expired. Configure the app with a larger value for --session-token-expiration if necessary
2026-04-23 18:57:35,899 - tornado.application - ERROR - Uncaught exception GET /status/ws (127.0.0.1)
HTTPServerRequest(protocol='http', host='jupyterhub.hpc.ucar.edu', method='GET', uri='/status/ws', version='HTTP/1.1', remote_ip='127.0.0.1')
Traceback (most recent call last):
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/websocket.py", line 965, in _accept_connection
    open_result = handler.open(*handler.open_args, **handler.open_kwargs)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/web.py", line 3388, in wrapper
    return method(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/bokeh/server/views/ws.py", line 149, in open
    raise ProtocolError("Token is expired. Configure the app with a larger value for --session-token-expiration if necessary")
bokeh.protocol.exceptions.ProtocolError: Token is expired. Configure the app with a larger value for --session-token-expiration if necessary
2026-04-23 18:57:38,467 - tornado.application - ERROR - Uncaught exception GET /status/ws (127.0.0.1)
HTTPServerRequest(protocol='http', host='jupyterhub.hpc.ucar.edu', method='GET', uri='/status/ws', version='HTTP/1.1', remote_ip='127.0.0.1')
Traceback (most recent call last):
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/websocket.py", line 965, in _accept_connection
    open_result = handler.open(*handler.open_args, **handler.open_kwargs)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/web.py", line 3388, in wrapper
    return method(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/bokeh/server/views/ws.py", line 149, in open
    raise ProtocolError("Token is expired. Configure the app with a larger value for --session-token-expiration if necessary")
bokeh.protocol.exceptions.ProtocolError: Token is expired. Configure the app with a larger value for --session-token-expiration if necessary
2026-04-23 19:00:18,540 - tornado.application - ERROR - Uncaught exception GET /status/ws (127.0.0.1)
HTTPServerRequest(protocol='http', host='jupyterhub.hpc.ucar.edu', method='GET', uri='/status/ws', version='HTTP/1.1', remote_ip='127.0.0.1')
Traceback (most recent call last):
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/websocket.py", line 965, in _accept_connection
    open_result = handler.open(*handler.open_args, **handler.open_kwargs)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/web.py", line 3388, in wrapper
    return method(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/bokeh/server/views/ws.py", line 149, in open
    raise ProtocolError("Token is expired. Configure the app with a larger value for --session-token-expiration if necessary")
bokeh.protocol.exceptions.ProtocolError: Token is expired. Configure the app with a larger value for --session-token-expiration if necessary
2026-04-23 19:00:27,137 - tornado.application - ERROR - Uncaught exception GET /status/ws (127.0.0.1)
HTTPServerRequest(protocol='http', host='jupyterhub.hpc.ucar.edu', method='GET', uri='/status/ws', version='HTTP/1.1', remote_ip='127.0.0.1')
Traceback (most recent call last):
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/websocket.py", line 965, in _accept_connection
    open_result = handler.open(*handler.open_args, **handler.open_kwargs)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/web.py", line 3388, in wrapper
    return method(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/bokeh/server/views/ws.py", line 149, in open
    raise ProtocolError("Token is expired. Configure the app with a larger value for --session-token-expiration if necessary")
bokeh.protocol.exceptions.ProtocolError: Token is expired. Configure the app with a larger value for --session-token-expiration if necessary
2026-04-23 19:00:44,304 - tornado.application - ERROR - Uncaught exception GET /status/ws (127.0.0.1)
HTTPServerRequest(protocol='http', host='jupyterhub.hpc.ucar.edu', method='GET', uri='/status/ws', version='HTTP/1.1', remote_ip='127.0.0.1')
Traceback (most recent call last):
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/websocket.py", line 965, in _accept_connection
    open_result = handler.open(*handler.open_args, **handler.open_kwargs)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/tornado/web.py", line 3388, in wrapper
    return method(self, *args, **kwargs)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/glade/u/home/harshah/.conda/envs/osdf/lib/python3.11/site-packages/bokeh/server/views/ws.py", line 149, in open
    raise ProtocolError("Token is expired. Configure the app with a larger value for --session-token-expiration if necessary")
bokeh.protocol.exceptions.ProtocolError: Token is expired. Configure the app with a larger value for --session-token-expiration if necessary

# Scale the cluster and display cluster dashboard URL
n_workers =5
cluster.scale(n_workers)
client.wait_for_workers(n_workers = n_workers)
cluster
Loading...

Open the catalog, find and load the variable of interest

%%time
era5_cat = intake.open_esm_datastore(era5_catalog)
jra3q_cat = intake.open_esm_datastore(jra3q_catalog)
era5_cat
CPU times: user 234 ms, sys: 34.7 ms, total: 268 ms
Wall time: 1.76 s
Loading...
# ── Time range ────────────────────────────────────────────────────────────────
YEAR_START_ERA5  = 1990
YEAR_START_JRA3Q = 1990
YEAR_END         = 2005

JJA_MONTHS = [6, 7, 8]

# ── Site coordinates (degrees north, degrees east 0–360) ──────────────────────
SITES = {
    "Boulder_CO" : {"lat":  40.01, "lon": 254.73},
    "Chicago_IL" : {"lat":  41.88, "lon": 272.37},
}

# ── Reproducibility ───────────────────────────────────────────────────────────
RANDOM_SEED = 42

  • Search for humidity and temperature variables in the catalog using the long_name and/or short_name columns. We only show the JRA-3Q example in the code below

jra3q_cat.df[['variable', 'short_name', 'long_name']].drop_duplicates().loc[
    jra3q_cat.df['long_name'].str.contains('2m temp|humid', case=False, na=False)]
Loading...
# ERA5: 2-m temperature and 2-m dew-point (GDEX convention)
ERA5_T_VAR  = "VAR_2T"
ERA5_TD_VAR = "VAR_2D"

era5_search = era5_cat.search(variable=[ERA5_T_VAR, ERA5_TD_VAR])
print(f"ERA5 matched: {len(era5_search.df):,} entries")

# JRA-3Q: 2-m temperature and 2-m relative humidity
# Variable names confirmed from the discovery step in the previous cell:
#   tmp2m-hgt-an-ll125  →  2-metre temperature, height level, analysis
#   rh2m-hgt-an-ll125   →  2-metre relative humidity, height level, analysis
JRA3Q_T_VAR  = "tmp2m-hgt-an-ll125"
JRA3Q_RH_VAR = "rh2m-hgt-an-ll125"

jra3q_search = jra3q_cat.search(variable=[JRA3Q_T_VAR, JRA3Q_RH_VAR])
print(f"JRA-3Q matched: {len(jra3q_search.df):,} entries")
ERA5 matched: 2 entries
JRA-3Q matched: 2 entries

jra3q_dsets = jra3q_search.to_dataset_dict(xarray_open_kwargs={'engine':'kerchunk',"chunks": {}})
# Select only required variables from JRA-3Q immediately
jra3q_dsets = {k: v[[JRA3Q_T_VAR, JRA3Q_RH_VAR]] for k, v in jra3q_dsets.items()}
# jra3q_dsets

--> The keys in the returned dictionary of datasets are constructed as follows:
	'variable.short_name'
Loading...
Loading...
%%time
era5_dsets = era5_search.to_dataset_dict()

--> The keys in the returned dictionary of datasets are constructed as follows:
	'variable.short_name'
Loading...
CPU times: user 749 ms, sys: 72.6 ms, total: 822 ms
Wall time: 18 s

#Inspect keys 
print(era5_dsets.keys())
jra3q_dsets.keys()
dict_keys(['VAR_2T.2t', 'VAR_2D.2d'])
dict_keys(['tmp2m-hgt-an-ll125.tmp2m-hgt-an-ll125', 'rh2m-hgt-an-ll125.rh2m-hgt-an-ll125'])

Extract data

  • Resample ERA5 data to 6hr frequency

  • Subset the data in space (Colorado and Chicago) and time

# Subset to analysis period and JJA
era5_points  = {}
jra3q_points = {}

for site, coords in SITES.items():
    lat, lon = coords["lat"], coords["lon"]

    # ERA5: merge variables, subset time and JJA
    era5_pt = xr.merge(
        [era5_dsets['VAR_2T.2t'][[ERA5_T_VAR]],
         era5_dsets['VAR_2D.2d'][[ERA5_TD_VAR]]]
    ).sel(time=slice(str(YEAR_START_ERA5), str(YEAR_END)))
    era5_pt = era5_pt.sel(time=era5_pt.time.dt.month.isin(JJA_MONTHS))

    # JRA-3Q: select variables, subset time and JJA
    jra3q_pt = jra3q_dsets[list(jra3q_dsets.keys())[0]][[JRA3Q_T_VAR, JRA3Q_RH_VAR]]
    jra3q_pt = jra3q_pt.sel(time=slice(str(YEAR_START_JRA3Q), str(YEAR_END)))
    jra3q_pt = jra3q_pt.sel(time=jra3q_pt.time.dt.month.isin(JJA_MONTHS))

    # Extract nearest grid cell for each site
    era5_pt   = era5_pt.sel(latitude=lat,  longitude=lon, method="nearest")
    jra3q_pt  = jra3q_pt.sel(lat=lat,      lon=lon,       method="nearest")

    print(f"ERA5   / {site}: requested ({lat:.2f}°N, {lon:.2f}°E) → "
          f"nearest ({float(era5_pt.latitude):.3f}°N, {float(era5_pt.longitude):.3f}°E)")
    print(f"JRA-3Q / {site}: requested ({lat:.2f}°N, {lon:.2f}°E) → "
          f"nearest ({float(jra3q_pt.lat):.3f}°N, {float(jra3q_pt.lon):.3f}°E)")

    era5_points[site]  = era5_pt
    jra3q_points[site] = jra3q_pt
ERA5   / Boulder_CO: requested (40.01°N, 254.73°E) → nearest (40.000°N, 254.750°E)
JRA-3Q / Boulder_CO: requested (40.01°N, 254.73°E) → nearest (40.000°N, 255.000°E)
ERA5   / Chicago_IL: requested (41.88°N, 272.37°E) → nearest (42.000°N, 272.250°E)
JRA-3Q / Chicago_IL: requested (41.88°N, 272.37°E) → nearest (42.500°N, 272.500°E)

# Randomly subsample ERA5 to one observation per 6-hour window
# to match the JRA-3Q temporal cadence prior to computing the annual maximum

def subsample_6h_random(ds, seed=RANDOM_SEED):
    """
    Retain one randomly selected observation per 6-hour window.
    The selection mask is built from the time coordinate only,
    leaving data variables Dask-lazy.
    """
    times        = pd.DatetimeIndex(ds["time"].values)
    window_start = times.floor("6h")
    unique_wins  = np.unique(window_start)

    rng           = np.random.default_rng(seed)
    offsets_h     = rng.integers(0, 6, size=len(unique_wins)).astype(int)
    win_to_offset = dict(zip(unique_wins, offsets_h))

    pos_h = ((times - window_start).total_seconds() / 3600).astype(int)
    keep  = np.fromiter(
        (pos_h[i] == win_to_offset[window_start[i]] for i in range(len(times))),
        dtype=bool,
        count=len(times),
    )

    print(f"  {ds.sizes['time']:,} → {keep.sum():,} time steps retained "
          f"({keep.sum() / len(times) * 100:.1f}%, expected ≈16.7%)")
    return ds.isel(time=keep)


for site in SITES:
    era5_points[site] = subsample_6h_random(era5_points[site])
    print(f"ERA5 / {site}: subsampling complete")
  35,328 → 5,888 time steps retained (16.7%, expected ≈16.7%)
ERA5 / Boulder_CO: subsampling complete
  35,328 → 5,888 time steps retained (16.7%, expected ≈16.7%)
ERA5 / Chicago_IL: subsampling complete

# Derive relative humidity from ERA5 2-m temperature and dew-point temperature.
# Following Romps (2024, supplementary section 1.5), relative humidity is computed
# as the ratio of the saturation vapor pressure evaluated at the dew-point temperature
# to that evaluated at the air temperature. The saturation vapor pressure is
# approximated using the August-Roche-Magnus formula (Alduchov and Eskridge, 1996).
#
# Alduchov, O.A. and Eskridge, R.E., 1996. Improved Magnus form approximation of
# saturation vapor pressure. Journal of Applied Meteorology and Climatology, 35(4),
# pp.601-609. DOI: 10.1175/1520-0450(1996)035<0601:IMFAOS>2.0.CO;2
#
# JRA-3Q provides relative humidity directly in %; convert to [0, 1].

def rh_from_T_Td(T_K, Td_K):
    a, b = 17.625, 243.04
    T_C  = T_K  - 273.15
    Td_C = Td_K - 273.15
    rh   = np.exp(a * Td_C / (b + Td_C)) / np.exp(a * T_C  / (b + T_C))
    return rh.clip(0.0, 1.0)


for site in SITES:
    era5_points[site] = era5_points[site].assign(
        rh=rh_from_T_Td(era5_points[site][ERA5_T_VAR],era5_points[site][ERA5_TD_VAR]))
    
    jra3q_points[site] = jra3q_points[site].chunk({"time": -1})
    jra3q_points[site] = jra3q_points[site].compute()
    
    jra3q_points[site] = jra3q_points[site].assign(rh=(jra3q_points[site][JRA3Q_RH_VAR] / 100.0).clip(0.0, 1.0))

print("Relative humidity prepared for all sites.")

Data Analysis

%%time
from heatindex import heatindex

def compute_heatindex_xr(T_K, rh):
    """
    Compute the extended heat index (Lu and Romps, 2022) using the
    Python implementation of Romps (2024).
    
    Parameters
    ----------
    T_K : xr.DataArray  — 2-m temperature in Kelvin
    rh  : xr.DataArray  — relative humidity in [0, 1]

    Returns
    -------
    xr.DataArray — heat index in Kelvin
    """
    hi_K = xr.apply_ufunc(
        heatindex,
        T_K,
        rh,
        dask="parallelized",
        output_dtypes=[float],
    )
    hi_K.attrs = {"units": "K", "long_name": "Heat Index (Lu and Romps, 2022)"}
    return hi_K


# Compute heat index and annual maximum JJA value for each site
results = {}

for site in SITES:
    results[site] = {}
    for label, points in [("ERA5", era5_points), ("JRA-3Q", jra3q_points)]:
        ds   = points[site]
        T_K  = ds[ERA5_T_VAR]  if label == "ERA5" else ds[JRA3Q_T_VAR]
        rh   = ds["rh"]

        hi_K     = compute_heatindex_xr(T_K, rh)
        ann_max  = hi_K.groupby("time.year").max(dim="time")

        print(f"Computing {label} / {site} ...")
        results[site][label] = (ann_max - 273.15).compute()
        print(f"  Done. Years: {int(results[site][label].year[0])}–"
              f"{int(results[site][label].year[-1])}")
%%time
fig, ax = plt.subplots(figsize=(10, 5))

colors = {"Boulder_CO": "tab:blue", "Chicago_IL": "tab:red"}
ls     = {"ERA5": "-", "JRA-3Q": "--"}

for site in SITES:
    for label, da in results[site].items():
        ax.plot(
            da.year,
            da.values,
            color     = colors[site],
            linestyle = ls[label],
            linewidth = 1.5,
            label     = f"{site.replace('_', ', ')} — {label}",
        )

ax.set_xlabel("Year")
ax.set_ylabel("Annual maximum JJA heat index (°C)")
ax.set_title("Annual Maximum JJA Heat Index\nBoulder, CO and Chicago, IL (1990–2005)")
ax.legend(framealpha=0.9)
ax.grid(linestyle=":", alpha=0.5)
plt.tight_layout()
plt.show()

  • The plot confirms our intuition: Chicago, a more humid place, has a higher JJA maximum heat index than Boulder

  • We see the infamous 1995 Chicago Heat Wave in the plot, but it is not captured by JRA-3Q!

  • For both Chicago and Boulder, the predictions from both the ERA5 and JRA-3Q reanalyses seem to broadly agree!

# Close the cluster
cluster.close()

References

  • Romps, D.M., 2024. Heat index extremes increasing several times faster than the air temperature, ERL, 2024 Environmental Research Letters

  • Lu, Y.-C. and Romps, D.M., 2022. Extending the heat index. Journal of Applied Meteorology and Climatology, 61(10), 1367–1383. DOI: Lu & Romps (2022)

  • Lu, Y.-C. and Romps, D.M., 2025. heatindex: Tools for Calculating Heat Stress, version 0.0.2. https://heatindex.org

  • Hersbach, H., et al., 2020. The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society, 146(730), 1999–2049. DOI: Hersbach et al. (2020)

  • Kosaka, Y., et al., 2024. The JRA-3Q Reanalysis. Journal of the Meteorological Society of Japan, 102, 49–109. DOI: KOSAKA et al. (2024)

References
  1. Lu, Y.-C., & Romps, D. M. (2022). Extending the Heat Index. Journal of Applied Meteorology and Climatology, 61(10), 1367–1383. 10.1175/jamc-d-22-0021.1
  2. Hersbach, H., Bell, B., Berrisford, P., Hirahara, S., Horányi, A., Muñoz‐Sabater, J., Nicolas, J., Peubey, C., Radu, R., Schepers, D., Simmons, A., Soci, C., Abdalla, S., Abellan, X., Balsamo, G., Bechtold, P., Biavati, G., Bidlot, J., Bonavita, M., … Thépaut, J. (2020). The ERA5 global reanalysis. Quarterly Journal of the Royal Meteorological Society, 146(730), 1999–2049. 10.1002/qj.3803
  3. KOSAKA, Y., KOBAYASHI, S., HARADA, Y., KOBAYASHI, C., NAOE, H., YOSHIMOTO, K., HARADA, M., GOTO, N., CHIBA, J., MIYAOKA, K., SEKIGUCHI, R., DEUSHI, M., KAMAHORI, H., NAKAEGAWA, T., TANAKA, T. Y., TOKUHIRO, T., SATO, Y., MATSUSHITA, Y., & ONOGI, K. (2024). The JRA-3Q Reanalysis. Journal of the Meteorological Society of Japan. Ser. II, 102(1), 49–109. 10.2151/jmsj.2024-004
NCAR data origin
Nino indices from HadISST data
Other data origins
Load data