Projected climate change signals and uncertainty under global warming levels

Projected climate change signals and uncertainty under global warming levels#

This notebook is a reproducibility example of the IPCC-WGI AR6 Interactive Atlas products. This work is licensed under a Creative Commons Attribution 4.0 International License.

E. Cimadevilla and M. Iturbide (Santander Meteorology Group. Instituto de Física de Cantabria, CSIC-UC, Santander, Spain).

This notebook is an example of the calculation and visualization of the IPCC-WGI AR6 uncertainty methods (simple and advanced) for projected delta changes. Please refer to the AR6 WGI Cross-Chapter Box Atlas 1 (Gutiérrez et al., 2021) for more information. We also introduce the Global Warming Level dimension for the analysis of the climate change signal.

See also hatching-uncertainty_R.ipynb at the IPCC-WGI/Atlas GitHub repository (also available in the hub: shared/repositories/IPCC-WGI-Atlas/notebooks).

This notebook works with the data available in this Hub, this is the dataset underlying the IPCC-WGI AR6 Interactive Atlas, originally published at DIGITAL.CSIC for the long-term archival, and also available through the Copernicus Data Store (CDS). Open the Getting_started.ipynb for a quick description of the available data.

Contents in this notebook#

Libraries
The Global Warming Level analysis dimension
Data loading for the different GWLs
Data loading for the historical reference
Uncertainty calculation and representation

1. Libraries#

import re
import math

import xarray
import numpy as np
import pandas as pd
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
import requests

2. The Global Warming Level analysis dimension#

Instead of calculating climate change anomalies for a fixed period, we will calculate the anomaly for a given level of global warming (GWL). To do so, we need the information on the time windows where the global surface temperature reaches the different levels of warming. This information is available at the IPCC-WGI/Atlas GitHub repository.

url = "https://github.com/SantanderMetGroup/ATLAS/raw/refs/heads/main/warming-levels/CMIP6_Atlas_WarmingLevels.csv"
gwls_file = url.split("/")[-1]
with open(gwls_file, "w") as f:
    response = requests.get(url)
    response.raise_for_status()
    f.write(response.text)
gwls = pd.read_csv(gwls_file)

gwls.head()

	model_run	1.5_ssp126	2_ssp126	3_ssp126	4_ssp126	1.5_ssp245	2_ssp245	3_ssp245	4_ssp245	1.5_ssp370	2_ssp370	3_ssp370	4_ssp370	1.5_ssp585	2_ssp585	3_ssp585	4_ssp585
0	ACCESS-CM2_r1i1p1f1	2027.0	2042.0	NaN	NaN	2028	2040	2070.0	NaN	2027	2039	2062.0	2082.0	2025	2038	2055	2071.0
1	ACCESS-ESM1-5_r1i1p1f1	2030.0	2072.0	NaN	NaN	2029	2045	NaN	NaN	2033	2048	2069.0	NaN	2027	2039	2060	2078.0
2	AWI-CM-1-1-MR_r1i1p1f1	2022.0	2050.0	NaN	NaN	2020	2039	NaN	NaN	2021	2037	2064.0	NaN	2019	2036	2059	2079.0
3	BCC-CSM2-MR_r1i1p1f1	2041.0	NaN	NaN	NaN	2035	2057	NaN	NaN	2032	2046	2074.0	NaN	2031	2043	2065	NaN
4	CAMS-CSM1-0_r2i1p1f1	NaN	NaN	NaN	NaN	2055	2088	NaN	NaN	2046	2065	NaN	NaN	2047	2064	2088	NaN

Note that we have pointed to the CMIP6_Atlas_WarmingLevels.csv file, as we are going to work with CMIP6 data. In this example, we will focus on the +3ºC GWL. We will consider the ssp585 scenario, however, any other scenario can be considered, as the anomalies do not vary significantly across scenarios when the GWL dimension is analyzed.

gwls3 = gwls[["model_run", "3_ssp585"]]
gwls3.head()

	model_run	3_ssp585
0	ACCESS-CM2_r1i1p1f1	2055
1	ACCESS-ESM1-5_r1i1p1f1	2060
2	AWI-CM-1-1-MR_r1i1p1f1	2059
3	BCC-CSM2-MR_r1i1p1f1	2065
4	CAMS-CSM1-0_r2i1p1f1	2088

Object gwl3 is a data frame containing the central year in a time window of 20 years where the +3ºC GWL is reached for each model (Find more information in the IPCC-WGI/Atlas repository).

5. Uncertainty calculation and representation#

climate4R implements the “simple” and “advanced” methods for the uncertainty characterization defined in the IPCC Sixth Assessment report. The function to apply is computeUncertainty. Please refer to the AR6 WGI Cross-Chapter Box Atlas 1 (Gutiérrez et al., 2021) for more information.

def model_agreement(da, axis, th = 80):
    nmembers, nlat, nlon = da.shape
    mask = np.array([
        (da[:,i,j] > 0).sum() > int(nmembers * th / 100)
        if da[:,i,j].mean() > 0
        else (da[:,i,j] < 0).sum() > int(nmembers * th / 100)
        for i in range(nlat)
        for j in range(nlon)]).reshape((nlat,nlon))
    
    return mask

def hatching(plot, mask, data):
    rows, cols = mask.shape
    for i in range(rows):
        for j in range(cols):
            lat, lon = data["lat"][i].item(), data["lon"][j].item()
            if mask[i,j]:
                plot.axes.plot([lon-.5,lon+.5],[lat+.5,lat-.5],'-',c="black", linewidth=.5)

mask_simple = delta.reduce(model_agreement, "member")

plot = ens_mean.plot(
    add_colorbar=True,
    cmap="BrBG",
    cbar_kwargs={"shrink": .5},
    vmin=-50, vmax=50,
    subplot_kws=dict(projection=ccrs.PlateCarree(central_longitude=0)),
    transform=ccrs.PlateCarree())

hatching(plot, ~mask_simple, ens_mean)

plot.axes.coastlines()

plt.savefig("uncer-europe.svg", bbox_inches="tight")

../_images/95917779dd2421cfd3c78d6c1d04ce7aba39de2a4ab3c1524a2793b4d0404dea.png

Now we compute the “advanced” method:

def advanced_method(historical, delta):
    vth = historical.reduce(
        lambda da, axis: ((math.sqrt(2) * 1.645 * np.std(da, axis=axis, ddof=1)) / math.sqrt(20)),
        "time")
    anom_abs = np.abs(delta)
    sig = anom_abs - vth
    si = xarray.where(sig > 0, 1, 0)
    uncer3 = delta.reduce(model_agreement, "member", th=80)
    uncer1 = si.reduce(
        lambda da, axis, nmembers, th: (da.sum(axis)/nmembers)*100 > th,
        "member",
        nmembers=si["member"].size,
        th=66)
    uncer2 = si.reduce(
        lambda da, axis, nmembers, th: (da.sum(axis)/nmembers)*100 < th,
        "member",
        nmembers=si["member"].size,
        th=66)
    uncer23 = (uncer2.astype(int) + uncer3.astype(int)) > 0
    uncer_a_aux1 = ((uncer1.astype(int) - 1) * -1).values
    uncer_a_aux2 = ((uncer23.astype(int) - 1) * -1).values
    uncer_a_aux2[uncer_a_aux2 == 1] = 2
    uncer_a_aux = uncer_a_aux1 + uncer_a_aux2
    uncer_a_aux[uncer_a_aux>1] = 2

    return uncer_a_aux

a = cmip6_hist["pr"].resample(time="QS-DEC").mean()
historical = a.sel(time=a.time.dt.month==12)

mask_advanced = advanced_method(historical, delta)

plot = ens_mean.plot(
    add_colorbar=True,
    cmap="BrBG",
    cbar_kwargs={"shrink": .5},
    vmin=-50, vmax=50,
    subplot_kws=dict(projection=ccrs.PlateCarree(central_longitude=0)),
    transform=ccrs.PlateCarree())

hatching(plot, mask_advanced, ens_mean)

plot.axes.coastlines()

<cartopy.mpl.feature_artist.FeatureArtist at 0x7adac02a2f90>

../_images/1688fb0bf5fa14c04c61067e91e9ea9d34a737d60ae922073e37756294fe930c.png

fig, axes = plt.subplots(1,2, figsize=(12,6), subplot_kw=dict(projection=ccrs.PlateCarree(central_longitude=0)))

plot_simple = ens_mean.plot(
    ax=axes[0],
    add_colorbar=True,
    cmap="BrBG",
    cbar_kwargs={"shrink": .5},
    vmin=-50, vmax=50,
    transform=ccrs.PlateCarree())

plot_advanced = ens_mean.plot(
    ax=axes[1],
    add_colorbar=True,
    cmap="BrBG",
    cbar_kwargs={"shrink": .5},
    vmin=-50, vmax=50,
    transform=ccrs.PlateCarree())

hatching(plot_simple, ~mask_simple, ens_mean)
hatching(plot_advanced, mask_advanced, ens_mean)

for ax in axes.flatten():
    ax.coastlines()
    ax.set_extent((lons.start, lons.stop, lats.start, lats.stop), ccrs.PlateCarree())

plt.savefig("uncertainty.svg")

../_images/b8c446c029932f7df07cdf4bd0e3bbeb7bac5a77339bd838c03a70ddf63f6c19.png

Regional stripes#

Here is an example of generating regional stripes.

cmip6_ssp585_y = cmip6_ssp585["pr"].resample(time="QS-DEC").mean()
cmip6_ssp585_y = cmip6_ssp585_y.sel(time=cmip6_ssp585_y["time"].dt.month==12)

cmip6_hist_c = cmip6_hist["pr"].mean(["time"])

year_delta = cmip6_ssp585_y - cmip6_hist_c
year_delta_rel = year_delta / cmip6_hist_c * 100

regional_mean = year_delta_rel.mean(["lat", "lon"])

regional_mean.plot.imshow(
    figsize=(14,7),
    add_colorbar=True,
    cmap="BrBG",
    vmin=-50,vmax=50)
plt.savefig("regional-stripes-pr.pdf", bbox_inches="tight")

../_images/39f16297e7e812f2a2e6a759a8631a759ed80e9f85e6101d38ee2b06553edbd3.png

cmip6_hist_y = cmip6_hist["pr"].resample(time="QS-DEC").mean()
cmip6_hist_y = cmip6_hist_y.sel(time=cmip6_hist_y["time"].dt.month==12)

cmip6_hist_y.mean(["lat", "lon"]).plot.imshow(
    figsize=(14,7),
    add_colorbar=True,
    cmap="YlGnBu",
    vmin=0,vmax=5)

<matplotlib.image.AxesImage at 0x7adaaa10cfd0>

../_images/9bb5e2dfa2309539b233f8749e02b052bd3dada15bfcfe55dc85d7d20333f1dd.png

We may omit the color bar to display only the spawn of each model time series.

regional_mean.where(np.isnan(regional_mean), 0).plot.imshow(
    figsize=(14,7),
    add_colorbar=False)

<matplotlib.image.AxesImage at 0x7adaaa141b90>

../_images/f46b03a6ac8653694bc0cc6740298fd3d0a87ce64b9aa27bbc4e00f09bcbbb60.png