Reproducing results from Hart et al., 2025

Here, we use climepi to reproduce results from Hart et al., PNAS, 2025 (https://doi.org/10.1073/pnas.2507311122).

Note: Running the code in step 1 will trigger the download of ~130 MB of climate projection data.

import climepi  # noqa
from climepi import climdata, epimod

2025-12-29 12:24:28,256 [INFO]: rcparams.py(_validate_backend:241) >> Found 'auto' as default backend, checking available backends
2025-12-29 12:24:28,256 [INFO]: rcparams.py(_validate_backend:241) >> Found 'auto' as default backend, checking available backends
2025-12-29 12:24:28,257 [INFO]: rcparams.py(_validate_backend:247) >> Matplotlib is available, defining as default backend
2025-12-29 12:24:28,257 [INFO]: rcparams.py(_validate_backend:247) >> Matplotlib is available, defining as default backend
2025-12-29 12:24:29,355 [INFO]: preview.py(<module>:14) >> arviz_base available, exposing its functions as part of arviz.preview
2025-12-29 12:24:29,355 [INFO]: preview.py(<module>:14) >> arviz_base available, exposing its functions as part of arviz.preview
2025-12-29 12:24:29,399 [INFO]: preview.py(<module>:33) >> arviz_stats available, exposing its functions as part of arviz.preview
2025-12-29 12:24:29,399 [INFO]: preview.py(<module>:33) >> arviz_stats available, exposing its functions as part of arviz.preview
2025-12-29 12:24:29,463 [INFO]: preview.py(<module>:47) >> arviz_plots available, exposing its functions as part of arviz.preview
2025-12-29 12:24:29,463 [INFO]: preview.py(<module>:47) >> arviz_plots available, exposing its functions as part of arviz.preview

1. Loading climate data

We load daily global climate projections for 2030-2100 in five cities from the ISIMIP project. The data are included as an example dataset (stored on the climepi GitHub repository) accessible via the get_example_dataset() method of the climdata subpackage, but can also be downloaded from the original source using climdata.get_climate_data() as follows:

ds_clim = climdata.get_climate_data(
    data_source="isimip",
    frequency="daily",
    subset={
            "years": list(range(2030, 2101)),
            "locations": ["London", "Paris", "Istanbul", "Cape Town", "Los Angeles"],
            "lon": [-0.08, 2.35, 28.98, 18.42, -118.42],
            "lat": [51.51, 48.86, 41.01, -33.93, 33.94],
        },
    save_dir="some/directory",
)

By default, the data are downloaded to the OS cache directory; change the ‘base_dir’ argument below to use a different file path.

ds_clim = climdata.get_example_dataset("isimip_cities_daily", base_dir=None)
ds_clim

<xarray.Dataset> Size: 32MB
Dimensions:        (time: 25932, bnds: 2, location: 5, scenario: 3, model: 10,
                    realization: 1)
Coordinates:
  * time           (time) datetime64[ns] 207kB 2030-01-01T12:00:00 ... 2100-1...
  * location       (location) <U11 220B 'London' 'Paris' ... 'Los Angeles'
  * scenario       (scenario) <U6 72B 'ssp126' 'ssp370' 'ssp585'
  * model          (model) <U13 520B 'gfdl-esm4' 'ipsl-cm6a-lr' ... 'miroc6'
  * realization    (realization) int64 8B 0
    lat            (location) float64 40B dask.array<chunksize=(1,), meta=np.ndarray>
    lon            (location) float64 40B dask.array<chunksize=(1,), meta=np.ndarray>
Dimensions without coordinates: bnds
Data variables:
    time_bnds      (time, bnds) datetime64[ns] 415kB 2030-01-01 ... 2101-01-01
    temperature    (location, scenario, model, realization, time) float32 16MB dask.array<chunksize=(1, 1, 1, 1, 25932), meta=np.ndarray>
    precipitation  (location, scenario, model, realization, time) float32 16MB dask.array<chunksize=(1, 1, 1, 1, 25932), meta=np.ndarray>
    lon_bnds       (location, bnds) float64 80B dask.array<chunksize=(1, 2), meta=np.ndarray>
    lat_bnds       (location, bnds) float64 80B dask.array<chunksize=(1, 2), meta=np.ndarray>
Attributes:
    title:        ISIMIP3b bias-adjusted climate input data
    institution:  Potsdam Institute for Climate Impact Research (PIK)
    project:      Inter-Sectoral Impact Model Intercomparison Project phase 3...
    contact:      ISIMIP cross-sectoral science team <info@isimip.org> <https...
    summary:      CMIP6 daily output data bias-adjusted and statistically dow...
    references:   Lange (2019) <https://doi.org/10.5194/gmd-12-3055-2019> and...

xarray.Dataset

• Dimensions:
  • time: 25932
  • bnds: 2
  • location: 5
  • scenario: 3
  • model: 10
  • realization: 1
• Coordinates: (7)
  • time
    (time)
    datetime64[ns]
    2030-01-01T12:00:00 ... 2100-12-...

    standard_name :

    time

    long_name :

    Time

    axis :

    T

    bounds :

    time_bnds

    array(['2030-01-01T12:00:00.000000000', '2030-01-02T12:00:00.000000000',
           '2030-01-03T12:00:00.000000000', ..., '2100-12-29T12:00:00.000000000',
           '2100-12-30T12:00:00.000000000', '2100-12-31T12:00:00.000000000'],
          shape=(25932,), dtype='datetime64[ns]')

  • location
    (location)
    <U11
    'London' 'Paris' ... 'Los Angeles'

    array(['London', 'Paris', 'Istanbul', 'Cape Town', 'Los Angeles'], dtype='<U11')

  • scenario
    (scenario)
    <U6
    'ssp126' 'ssp370' 'ssp585'

    array(['ssp126', 'ssp370', 'ssp585'], dtype='<U6')

  • model
    (model)
    <U13
    'gfdl-esm4' ... 'miroc6'

    array(['gfdl-esm4', 'ipsl-cm6a-lr', 'mpi-esm1-2-hr', 'mri-esm2-0',
           'ukesm1-0-ll', 'canesm5', 'cnrm-cm6-1', 'cnrm-esm2-1', 'ec-earth3',
           'miroc6'], dtype='<U13')

  • realization
    (realization)
    int64
    0

    array([0])

  • lat
    (location)
    float64
    dask.array<chunksize=(1,), meta=np.ndarray>

    standard_name :

    latitude

    long_name :

    Latitude

    units :

    °N

    axis :

    Y

    bounds :

    lat_bnds

    Array Chunk Bytes 40 B 8 B Shape (5,) (1,) Dask graph 5 chunks in 12 graph layers Data type float64 numpy.ndarray

  • lon
    (location)
    float64
    dask.array<chunksize=(1,), meta=np.ndarray>

    standard_name :

    longitude

    long_name :

    Longitude

    units :

    °E

    axis :

    X

    bounds :

    lon_bnds

    Array Chunk Bytes 40 B 8 B Shape (5,) (1,) Dask graph 5 chunks in 12 graph layers Data type float64 numpy.ndarray

• Data variables: (5)
  • time_bnds
    (time, bnds)
    datetime64[ns]
    2030-01-01 ... 2101-01-01

    xcdat_bounds :

    True

    array([['2030-01-01T00:00:00.000000000', '2030-01-02T00:00:00.000000000'],
           ['2030-01-02T00:00:00.000000000', '2030-01-03T00:00:00.000000000'],
           ['2030-01-03T00:00:00.000000000', '2030-01-04T00:00:00.000000000'],
           ...,
           ['2100-12-29T00:00:00.000000000', '2100-12-30T00:00:00.000000000'],
           ['2100-12-30T00:00:00.000000000', '2100-12-31T00:00:00.000000000'],
           ['2100-12-31T00:00:00.000000000', '2101-01-01T00:00:00.000000000']],
          shape=(25932, 2), dtype='datetime64[ns]')

  • temperature
    (location, scenario, model, realization, time)
    float32
    dask.array<chunksize=(1, 1, 1, 1, 25932), meta=np.ndarray>

    long_name :

    Temperature

    units :

    °C

    Array Chunk Bytes 14.84 MiB 101.30 kiB Shape (5, 3, 10, 1, 25932) (1, 1, 1, 1, 25932) Dask graph 150 chunks in 345 graph layers Data type float32 numpy.ndarray

  • precipitation
    (location, scenario, model, realization, time)
    float32
    dask.array<chunksize=(1, 1, 1, 1, 25932), meta=np.ndarray>

    long_name :

    Precipitation

    units :

    mm/day

    Array Chunk Bytes 14.84 MiB 101.30 kiB Shape (5, 3, 10, 1, 25932) (1, 1, 1, 1, 25932) Dask graph 150 chunks in 345 graph layers Data type float32 numpy.ndarray

  • lon_bnds
    (location, bnds)
    float64
    dask.array<chunksize=(1, 2), meta=np.ndarray>

    Array Chunk Bytes 80 B 16 B Shape (5, 2) (1, 2) Dask graph 5 chunks in 12 graph layers Data type float64 numpy.ndarray

  • lat_bnds
    (location, bnds)
    float64
    dask.array<chunksize=(1, 2), meta=np.ndarray>

    Array Chunk Bytes 80 B 16 B Shape (5, 2) (1, 2) Dask graph 5 chunks in 12 graph layers Data type float64 numpy.ndarray

• Attributes: (6)

  title :

  ISIMIP3b bias-adjusted climate input data

  institution :

  Potsdam Institute for Climate Impact Research (PIK)

  project :

  Inter-Sectoral Impact Model Intercomparison Project phase 3b (ISIMIP3b)

  contact :

  ISIMIP cross-sectoral science team <info@isimip.org> <https://www.isimip.org>

  summary :

  CMIP6 daily output data bias-adjusted and statistically downscaled to 0.5 degree horizontal resolution using ISIMIP3BASD v2.5.0 and W5E5 v2.0

  references :

  Lange (2019) <https://doi.org/10.5194/gmd-12-3055-2019> and Lange (2021) <https://doi.org/10.5281/zenodo.4686991> for ISIMIP3BASD; Cucchi et al. (2020) <https://doi.org/10.5194/essd-2020-28> and Lange et al. (2021) <https://doi.org/10.48364/ISIMIP.342217> for W5E5

2. Defining the epidemiological model

We use the temperature niche model for dengue transmission Aedes albopictus from Mordecai et al., PLOS Negl Trop Dis, 2017 (https://doi.org/10.1371/journal.pntd.0005568), which is included as an example model accessible via the get_example_model() method of the epimod subpackage.

suitability_model = epimod.get_example_model("mordecai_ae_albopictus_niche")

We can use the plot_suitability() method to visualize the model, which here simply comprises a range of temperatures within which transmission is assumed to be possible.

suitability_model.plot_suitability()

3. Running the model

We use the climepi accessor for xarray Datasets to run the epidemiological model on on the climate data, obtaining projections of the number of days suitable for transmission each year.

ds_months_suitable = ds_clim.climepi.run_epi_model(
    suitability_model, return_yearly_portion_suitable=True
)
ds_months_suitable

<xarray.Dataset> Size: 88kB
Dimensions:           (time: 71, location: 5, scenario: 3, model: 10,
                       realization: 1, bnds: 2)
Coordinates:
  * time              (time) object 568B 2030-07-02 12:00:00 ... 2100-07-02 1...
  * location          (location) <U11 220B 'London' 'Paris' ... 'Los Angeles'
  * scenario          (scenario) <U6 72B 'ssp126' 'ssp370' 'ssp585'
  * model             (model) <U13 520B 'gfdl-esm4' 'ipsl-cm6a-lr' ... 'miroc6'
  * realization       (realization) int64 8B 0
    lat               (location) float64 40B dask.array<chunksize=(1,), meta=np.ndarray>
    lon               (location) float64 40B dask.array<chunksize=(1,), meta=np.ndarray>
Dimensions without coordinates: bnds
Data variables:
    portion_suitable  (time, location, scenario, model, realization) int64 85kB dask.array<chunksize=(71, 1, 1, 1, 1), meta=np.ndarray>
    lat_bnds          (location, bnds) float64 80B dask.array<chunksize=(1, 2), meta=np.ndarray>
    lon_bnds          (location, bnds) float64 80B dask.array<chunksize=(1, 2), meta=np.ndarray>
    time_bnds         (time, bnds) object 1kB 2030-01-01 00:00:00 ... 2101-01...

xarray.Dataset

• Dimensions:
  • time: 71
  • location: 5
  • scenario: 3
  • model: 10
  • realization: 1
  • bnds: 2
• Coordinates: (7)
  • time
    (time)
    object
    2030-07-02 12:00:00 ... 2100-07-...

    standard_name :

    time

    long_name :

    Time

    axis :

    T

    bounds :

    time_bnds

    array([cftime.DatetimeProlepticGregorian(2030, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2031, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2032, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2033, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2034, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2035, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2036, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2037, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2038, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2039, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2040, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2041, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2042, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2043, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2044, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2045, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2046, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2047, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2048, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2049, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2050, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2051, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2052, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2053, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2054, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2055, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2056, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2057, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2058, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2059, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2060, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2061, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2062, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2063, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2064, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2065, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2066, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2067, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2068, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2069, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2070, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2071, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2072, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2073, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2074, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2075, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2076, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2077, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2078, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2079, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2080, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2081, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2082, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2083, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2084, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2085, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2086, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2087, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2088, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2089, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2090, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2091, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2092, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2093, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2094, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2095, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2096, 7, 2, 0, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2097, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2098, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2099, 7, 2, 12, 0, 0, 0, has_year_zero=True),
           cftime.DatetimeProlepticGregorian(2100, 7, 2, 12, 0, 0, 0, has_year_zero=True)],
          dtype=object)

  • location
    (location)
    <U11
    'London' 'Paris' ... 'Los Angeles'

    array(['London', 'Paris', 'Istanbul', 'Cape Town', 'Los Angeles'], dtype='<U11')

  • scenario
    (scenario)
    <U6
    'ssp126' 'ssp370' 'ssp585'

    array(['ssp126', 'ssp370', 'ssp585'], dtype='<U6')

  • model
    (model)
    <U13
    'gfdl-esm4' ... 'miroc6'

    array(['gfdl-esm4', 'ipsl-cm6a-lr', 'mpi-esm1-2-hr', 'mri-esm2-0',
           'ukesm1-0-ll', 'canesm5', 'cnrm-cm6-1', 'cnrm-esm2-1', 'ec-earth3',
           'miroc6'], dtype='<U13')

  • realization
    (realization)
    int64
    0

    array([0])

  • lat
    (location)
    float64
    dask.array<chunksize=(1,), meta=np.ndarray>

    standard_name :

    latitude

    long_name :

    Latitude

    units :

    °N

    axis :

    Y

    bounds :

    lat_bnds

    Array Chunk Bytes 40 B 8 B Shape (5,) (1,) Dask graph 5 chunks in 12 graph layers Data type float64 numpy.ndarray

  • lon
    (location)
    float64
    dask.array<chunksize=(1,), meta=np.ndarray>

    standard_name :

    longitude

    long_name :

    Longitude

    units :

    °E

    axis :

    X

    bounds :

    lon_bnds

    Array Chunk Bytes 40 B 8 B Shape (5,) (1,) Dask graph 5 chunks in 12 graph layers Data type float64 numpy.ndarray

• Data variables: (4)
  • portion_suitable
    (time, location, scenario, model, realization)
    int64
    dask.array<chunksize=(71, 1, 1, 1, 1), meta=np.ndarray>

    long_name :

    Days where suitability > 0

    Array Chunk Bytes 83.20 kiB 568 B Shape (71, 5, 3, 10, 1) (71, 1, 1, 1, 1) Dask graph 150 chunks in 354 graph layers Data type int64 numpy.ndarray

  • lat_bnds
    (location, bnds)
    float64
    dask.array<chunksize=(1, 2), meta=np.ndarray>

    Array Chunk Bytes 80 B 16 B Shape (5, 2) (1, 2) Dask graph 5 chunks in 12 graph layers Data type float64 numpy.ndarray

  • lon_bnds
    (location, bnds)
    float64
    dask.array<chunksize=(1, 2), meta=np.ndarray>

    Array Chunk Bytes 80 B 16 B Shape (5, 2) (1, 2) Dask graph 5 chunks in 12 graph layers Data type float64 numpy.ndarray

  • time_bnds
    (time, bnds)
    object
    2030-01-01 00:00:00 ... 2101-01-...

    xcdat_bounds :

    True

    array([[cftime.DatetimeProlepticGregorian(2030, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2031, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2031, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2032, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2032, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2033, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2033, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2034, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2034, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2035, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2035, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2036, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2036, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2037, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2037, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2038, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2038, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2039, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2039, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2040, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
    ...
            cftime.DatetimeProlepticGregorian(2092, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2092, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2093, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2093, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2094, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2094, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2095, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2095, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2096, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2096, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2097, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2097, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2098, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2098, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2099, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2099, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2100, 1, 1, 0, 0, 0, 0, has_year_zero=True)],
           [cftime.DatetimeProlepticGregorian(2100, 1, 1, 0, 0, 0, 0, has_year_zero=True),
            cftime.DatetimeProlepticGregorian(2101, 1, 1, 0, 0, 0, 0, has_year_zero=True)]],
          dtype=object)

4. Visualizing the results

We now reproduce Figure 1 from Hart et al., 2025, which shows uncertainty in the future number of months suitable, decomposed into contributions from internal climate variability, climate model uncertainty and scenario uncertainty. This is achieved by using the plot_uncertainty_interval_decomposition method from the climepi package. Since only one simulation is available for each model/scenario pair, we estimate the extent of internal variability for each pair based on year-to-year deviations around a cubic polynomial fit (keyword arguments internal_variability_method='polyfit' and deg=3; note these are the default settings when a single simulation is available per model/scenario pair, but are explicitly set here for clarity). For methodological details, see Hart et al., 2025, as well as Hawkins and Sutton, Bull Am Meteorol Soc, 2009 (https://doi.org/10.1175/2009BAMS2607.1), from which the uncertainty decomposition methodology is adapted.

The return value of plot_uncertainty_interval_decomposition() method is a HoloViews Overlay object. Customization options can be applied using the opts() method.

ds_months_suitable.sel(
    location="London"
).climepi.plot_uncertainty_interval_decomposition(
    uncertainty_level=90,
    internal_variability_method="polyfit",
    deg=3,
).opts(
    ylim=(0, 220),
    show_title=False,
    legend_position="top_left",
    legend_opts={
        "location": (5, 150),
        "label_text_font_size": "9pt",
        "padding": 4,
        "spacing": 1,
    },
)