Open In Colab   Open in Kaggle

Tutorial 1: Creating Maps of CMIP6 Earth System Model (ESM) Projections#

Week 2, Day 1, Future Climate: The Physical Basis

Content creators: Brodie Pearson, Julius Busecke, Tom Nicholas

Content reviewers: Younkap Nina Duplex, Zahra Khodakaramimaghsoud, Sloane Garelick, Peter Ohue, Jenna Pearson, Agustina Pesce, Derick Temfack, Peizhen Yang, Cheng Zhang, Chi Zhang, Ohad Zivan

Content editors: Jenna Pearson, Ohad Zivan, Chi Zhang

Production editors: Wesley Banfield, Jenna Pearson, Chi Zhang, Ohad Zivan

Our 2023 Sponsors: NASA TOPS, Google DeepMind, and CMIP

Tutorial Objectives#

Earth System Models (ESMs) provide physically-based projections of how Earth’s climate could change in the coming years, decades, and centuries at both global and local scales. In the following tutorial you will:

  • Revisit how to load ESM data from the CMIP6 experiments, and

  • Create maps showing projected future changes in sea surface temperature (SST).


# installations ( uncomment and run this cell ONLY when using google colab or kaggle )

# !pip install condacolab &> /dev/null
# import condacolab
# condacolab.install()

# # Install all packages in one call (+ use mamba instead of conda), this must in one line or code will fail
# !mamba install xarray-datatree intake-esm gcsfs xmip aiohttp cartopy nc-time-axis cf_xarray xarrayutils "esmf<=8.3.1" xesmf &> /dev/null
# # For xesmf install we need to pin "esmf<=8.3.1". More context here:
# imports
import time

tic = time.time()

import intake
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import xesmf as xe

from xmip.preprocessing import combined_preprocessing
from xarrayutils.plotting import shaded_line_plot

from datatree import DataTree
from xmip.postprocessing import _parse_metric

import as ccrs

Figure settings#

Hide code cell source
# @title Figure settings
import ipywidgets as widgets  # interactive display

%matplotlib inline

Video 1: Recap of Earth System Models#

Section 1.1: Loading CMIP6 SST Data with xarray#

As a reminder, these ESMs simulate several systems (ocean, atmosphere, cryosphere, land) that are coupled to each other, and each system has its own variables, physics, and discretizations (grid & timestep).


Atmospheric Model Schematic (Credit: Wikipedia)

Let’s repeat the CMIP6 loading method that we learned in Tutorial 6 on last week’s Climate Modelling day (Day 5).

Although we will only work with monthly SST (ocean) data today, the methods introduced can easily be applied/extended to load and analyze other CMIP6 variables, including from other components of the Earth system.

To learn more about CMIP, including additional methods to access CMIP data, please see our CMIP Resource Bank and the CMIP website.

As a reminder, the facets that have to be specified for CMIP6, along with the facet choice(s) we make in this tutorial, are:

  1. variable_id: tos = sea surface temperature

  2. source_id: The CMIP6 model(s) that we want data from

  3. table_id: Omon (ocean monthly output)

  4. grid_id: gn = data on the model’s native grid

  5. experiment_id: ssp585 (we’ll discuss experiments later today)

  6. member_id: r1i1p1f1 for now

Now, let’s repeat our CMIP6 loading method from the previous tutorial.

Note: today we will start by using only use one model, TaiESM1, which stands for Taiwan Earth System Model version 1, and a single experiment, ssp585 which is a high-emissions future scenario. In later tutorials you will work with 5 distinct CMIP6 models (source_id), and two additional experiments (experiment_id). TaiESM1 was developed by modifying the Community Earth System Model (CESM) version 1.2.2 to include different parameterizations (i.e., physics). As a result, the TaiESM1 model output is distinct from the CESM output you used in previous tutorials/days.

# open an intake catalog containing the Pangeo CMIP cloud data
col = intake.open_esm_datastore(

# from the full `col` object, create a subset using facet search
cat =
    ],  # make sure that we only get models which have all of the above experiments

# convert the sub-catalog into a datatree object, by opening each dataset into an xarray.Dataset (without loading the data)
kwargs = dict(
    preprocess=combined_preprocessing,  # apply xMIP fixes to each dataset
    ),  # ensure all datasets use the same time index
        "token": "anon"
    },  # anonymous/public authentication to google cloud storage

cat.esmcat.aggregation_control.groupby_attrs = ["source_id", "experiment_id"]
dt = cat.to_datatree(**kwargs)
--> The keys in the returned dictionary of datasets are constructed as follows:
100.00% [1/1 00:02<00:00]

We now have a “datatree” containing the data we searched for. A datatree is a high-level container of xarray data, useful for organizing many related datasets together. You can think of a single DataTree object as being like a (nested) dictionary of xarray.Dataset objects. Each dataset in the tree is known as a “node” or “group”, and we can also have empty nodes. This DataTree object may seem complicated when we load only one dataset, but it will prove to be very useful in later tutorials where you will work with multiple models, experiments, and ensemble members

You can explore the nodes of the tree and its contents interactively in a similar way to how you can explore the contents of an xarray.Dataset:

<xarray.DatasetView> Size: 0B
Dimensions:  ()
Data variables:

Now that we have the model datasets organized within thie datatree (dt) we can plot the datasets. Let’s start by plotting a map of SST from the TaiESM1 CMIP6 model in July 2023.

Note that CMIP6 experiments were run several years ago, so the cut-off between past (observed forcing) and future (scenario-based/projected forcing) was at the start of 2015. This means that July 2023 is about 8 years into the CMIP6 future and so it is unlikely to look exactly like Earth’s current SST state.

# select just a single model (TaiESM1) and experiment (ssp585) to plot
sst_ssp585 = dt["TaiESM1"]["ssp585"].ds.tos

fig, (ax_present) = plt.subplots(
    ncols=1, nrows=1, figsize=[8, 4], subplot_kw={"projection": ccrs.Robinson()}

# select the model data for July 2023
sst_present = sst_ssp585.sel(time="2023-07").squeeze()

# plot the model data
ax_present.set_title("July 2023")
Text(0.5, 1.0, 'July 2023')

Coding Exercise 1.1#

Now that we can plot maps of CMIP6 data, let’s look at some projected future changes using this data!

In this coding exercise your goals are to:

  1. Create a map of the projected sea surface temperature in July 2100 under the SSP5-8.5 high-emissions scenario (we’ll discuss scenarios in the next mini-lecture) using data from the TaiESM1 CMIP6 model.

  2. Create a map showing how this sea surface temperature projection is different from the current (July 2023) sea surface temperature in this model

  3. Plot a similar map for this model that shows how January 2100 is different from January 2023

To get you started, we have provided code to load the required data set into a variable called sst_ssp585, and we have plotted the current (July 2023) sea surface temperature from this data set.

Note: differences between two snapshots of SST are not the same as the anomalies that you encountered earlier in the course, which were the difference relative to the average during a reference period.

# select just a single model and experiment
sst_ssp585 = dt["TaiESM1"]["ssp585"].ds.tos

fig, ([ax_present, ax_future], [ax_diff_july, ax_diff_jan]) = plt.subplots(
    ncols=2, nrows=2, figsize=[12, 6], subplot_kw={"projection": ccrs.Robinson()}

# plot a timestep for 2023
sst_present = sst_ssp585.sel(time="2023-07").squeeze()
ax_present.set_title("July 2023")

# repeat for 2100
# complete the following line to extract data for July 2100
sst_future = ...
_ = ...
ax_future.set_title("July 2100")

# now find the difference between July 2100 and July 2023
# complete the following line to extract the July difference
sst_difference_july = ...
_ = ...
ax_diff_july.set_title("2100 vs. 2023 Difference (July)")

# finally, find the difference between January of the two years used above
# complete the following line to extract the January difference
sst_difference_jan = ...
_ = ...
ax_diff_jan.set_title("2100 vs. 2023 Difference (January)")
Text(0.5, 1.0, '2100 vs. 2023 Difference (January)')

Click for solution

Example output:

Solution hint

Questions 1.1: Climate Connection#

  1. Comparing only the top two panels, how is the July SST projected to change in this particular model simulation? Do these changes agree with the map of July change that you plotted in the bottom left, and are these changes easier to see in this bottom map?

  2. In what ways are the July and January maps similar or dissimilar, and can you think of any physical explanations for these (dis)similarities?

  3. Why do you think the color bar axes vary? (i.e., the top panels say “Sea Surface Temperature [\(^oC\)]” while the bottom panels say “tos”)

Many of the changes seen in the maps are a result of a changing climate under this high-emissions scenarios. However, keep in mind that these are differences between two months that are almost 80 years apart, so some of the changes are due to weather/synoptic differences between these particular months.

Click for solution

Section 1.2: Horizontal Regridding#

Many CMIP6 models use distinct spatial grids, we call this the model’s native grid.

You are likely familiar with the regular latitude-longitude grid where we separate the planet into boxes that have a fixed latitude and longitude span like this image we saw in the tutorial:


Let’s look at the grid used for the ocean component of the TaiESM1 CMIP6 model:

# create a scatter plot with a symbol at the center of each ocean grid cell in TaiESM1
fig, ax = plt.subplots()
ax.scatter(x=sst_ssp585.lon,, s=0.1)
ax.set_title("Grid cell locations in TaiESM1")
Text(0.5, 1.0, 'Grid cell locations in TaiESM1')

Questions 1.2#

  1. How would this plot look for a regular latitude-longitude grid like the globe image shown above and in the slides? In what ways is the TaiESM1 grid different from this regular grid?

  2. Can you think of a reason the Northern and Southern hemisphere ocean grids differ?*

*Hint: from an oceanographic context, how are the North and South poles different from each other?

Click for solution

If you want to compare spatial maps from different models/observations, e.g. plot a map averaged over several models or the bias of this map relative to observations, you must first ensure the data from all the models (and observations) is on the same spatial grid. This is where regridding becomes essential!

Regridding is applied lazily, but it is still taking time to compute when it is applied. So if you want to compare for example the mean over time of several models it is often much quicker to compute the mean in time over the native grid and then regrid the result of that, instead of regridding each timestep and then calculating the mean!

# define a 'target' grid. This is simply a regular lon/lat grid that we will interpolate our data on
ds_target = xr.Dataset(
        "lat": (["lat"], np.arange(-90, 90, 1.0), {"units": "degrees_north"}),
        "lon": (["lon"], np.arange(0, 360, 1.0), {"units": "degrees_east"}),
)  # you can try to modify the parameters above to e.g. just regrid onto a region or make the resolution coarser etc
<xarray.Dataset> Size: 4kB
Dimensions:  (lat: 180, lon: 360)
  * lat      (lat) float64 1kB -90.0 -89.0 -88.0 -87.0 ... 86.0 87.0 88.0 89.0
  * lon      (lon) float64 3kB 0.0 1.0 2.0 3.0 4.0 ... 356.0 357.0 358.0 359.0
Data variables:
# define the regridder object (from our source dataarray to the target)
regridder = xe.Regridder(
    sst_ssp585, ds_target, "bilinear", periodic=True
)  # this takes some time to calculate a weight matrix for the regridding
xESMF Regridder 
Regridding algorithm:       bilinear 
Weight filename:   
Reuse pre-computed weights? False 
Input grid shape:           (384, 320) 
Output grid shape:          (180, 360) 
Periodic in longitude?      True
# now we can apply the regridder to our data
sst_ssp585_regridded = regridder(sst_ssp585)  # this is a lazy operation!
# so it does not slow us down significantly to apply it to the full data!
# we can work with this array just like before and the regridding will only be
# applied to the parts that we later load into memory or plot.
<xarray.DataArray (member_id: 1, dcpp_init_year: 1, time: 1032, lat: 180,
                   lon: 360)> Size: 267MB
dask.array<astype, shape=(1, 1, 1032, 180, 360), dtype=float32, chunksize=(1, 1, 87, 180, 320), chunktype=numpy.ndarray>
  * time            (time) object 8kB 2015-01-17 00:29:59.999993 ... 2100-12-...
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) float64 8B nan
  * lat             (lat) float64 1kB -90.0 -89.0 -88.0 -87.0 ... 87.0 88.0 89.0
  * lon             (lon) float64 3kB 0.0 1.0 2.0 3.0 ... 357.0 358.0 359.0
    regrid_method:  bilinear
# compare the shape to the original array
<xarray.DataArray 'tos' (member_id: 1, dcpp_init_year: 1, time: 1032, y: 384,
                         x: 320)> Size: 507MB
dask.array<broadcast_to, shape=(1, 1, 1032, 384, 320), dtype=float32, chunksize=(1, 1, 87, 384, 320), chunktype=numpy.ndarray>
    lat             (y, x) float64 983kB dask.array<chunksize=(384, 320), meta=np.ndarray>
    lon             (y, x) float64 983kB dask.array<chunksize=(384, 320), meta=np.ndarray>
  * time            (time) object 8kB 2015-01-17 00:29:59.999993 ... 2100-12-...
  * y               (y) int64 3kB 0 1 2 3 4 5 6 ... 377 378 379 380 381 382 383
  * x               (x) int64 3kB 0 1 2 3 4 5 6 ... 313 314 315 316 317 318 319
  * member_id       (member_id) object 8B 'r1i1p1f1'
  * dcpp_init_year  (dcpp_init_year) float64 8B nan
    cell_measures:  area: areacello
    cell_methods:   area: mean where sea time: mean
    comment:        Temperature of upper boundary of the liquid ocean, includ...
    history:        2021-04-16T01:30:38Z altered by CMOR: replaced missing va...
    long_name:      Sea Surface Temperature
    original_name:  TEMP
    standard_name:  sea_surface_temperature
    units:          °C

Section 1.3: Visually Comparing Data with Different Map Projections#

Let’s use the code from above to plot a map of the model data on its original (native) grid, and a map of the model data after it is regridded.

fig, ([ax_regridded, ax_native]) = plt.subplots(
    ncols=2, figsize=[12, 3], subplot_kw={"projection": ccrs.Robinson()}

# Native grid data
sst_future = sst_ssp585.sel(time="2100-07").squeeze()
ax_native.set_title("July 2100 Native Grid")

# Regridded data
sst_future_regridded = sst_ssp585_regridded.sel(time="2100-07").squeeze()
ax_regridded.set_title("July 2100 Regridded")
Text(0.5, 1.0, 'July 2100 Regridded')

Questions 1.3#

  1. Is this what you expected to see after regridding the data?

Click for solution


In Tutorial 1 you have:

  • Loaded and manipulated data from a CMIP6 model under a high-emissions future scenario experiment

  • Created maps of future projected changes in the Earth system using CMIP6 data

  • Converted/regridded CMIP6 model data onto a desired grid. This is a critical processing step that allows us to directly compare data from different models and/or observations


This tutorial uses data from the simulations conducted as part of the CMIP6 multi-model ensemble.

For examples on how to access and analyze data, please visit the Pangeo Cloud CMIP6 Gallery

For more information on what CMIP is and how to access the data, please see this page.