logo

Tutorial on July 2023 record-breaking global surface temperatures using climate data from C3S#

This notebook can be run on free online platforms, such as Binder, Kaggle and Colab, or they can be accessed from GitHub. The links to run this notebook in these environments are provided here, but please note they are not supported by ECMWF.

binder kaggle colab github

Example of work

Learning objectives 🎯#

In this tutorial we will access data from the Climate Data Store (CDS) of the Copernicus Climate Change Service (C3S), and analyse air and sea surface temperatures comparing July 2023 record breaking values with climatologies. We will:

  • Calculate a global surface temperature climatology

  • Compute and visualise anomalies with respect to the climatology

  • View time series and rank global surface temperature records

  • Analyse North Atlantic sea surface temperature trends

Prepare your environment#

Set up the CDSAPI and your credentials#

The code below will ensure that the cdsapi package is installed. If you have not setup your ~/.cdsapirc file with your credenials, you can replace None with your credentials that can be found on the how to api page (you will need to log in to see your credentials).

!pip install -q cdsapi
# If you have already setup your .cdsapirc file you can leave this as None
cdsapi_key = None
cdsapi_url = None

(Install and) Import libraries#

We will be working with data in NetCDF format. To best handle this data we will use libraries for working with multidimensional arrays, in particular Xarray. We will also need libraries for plotting and viewing data, in this case we will use Matplotlib and Cartopy.

import urllib3
urllib3.disable_warnings()

import os

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cdsapi

plt.style.use('bmh')

Specify data directory#

# Directory to store data
# Please ensure that data_dir is a location where you have write permissions
DATADIR = './data_dir/'
# Create this directory if it doesn't exist
os.makedirs(DATADIR, exist_ok=True)

Explore data#

To search for data, visit the CDS website: https://cds.climate.copernicus.eu. Here you can search for ERA5 data using the search bar. The data we need for this tutorial is the ERA5 monthly averaged data on single levels from 1940 to present. ERA5 is the 5th version of the ECMWF Reanalysis dataset. Reanalysis uses a state of the art forecast model and data assimilation system to create a consistent “map without gaps” of observed and modelled climate variables over the past decades.

Search for the data#

Having selected the correct dataset, we now need to specify what product type, variables, temporal and geographic coverage we are interested in. These can all be selected in the “Download data” tab. In this tab a form appears in which we will select the following parameters to download:

Parameters of data to download
  • Product type: Monthly averaged reanalysis

  • Variable: 2m temperature

  • Year: 1940 to present year

  • Month: july

  • Time: 00:00 (default)

  • Geographical area: Sub-region extraction and leave the default values, 'area': [90, -180, -90, 180], # maxlat, minlon, minlat, maxlon (this will be explained further below)

  • Format: NetCDF

At the end of the download forms, select “Show API request”. This will reveal a block of code, which you can simply copy and paste into a cell of your Jupyter Notebook (see cell below). Having copied the API request into the cell below, running this will retrieve and download the data you requested into your local directory.

Warning

Please remember to accept the terms and conditions of the dataset, at the bottom of the CDS download form!

Download the data#

With the API request copied into the cell below, running this cell will retrieve and download the data you requested into your local directory.

c = cdsapi.Client()
c.retrieve(
    'reanalysis-era5-single-levels-monthly-means',
    {
        'data_format': 'netcdf',
        'product_type': 'monthly_averaged_reanalysis',
        'variable': '2m_temperature',
        'year': [
            '1940', '1941', '1942',
            '1943', '1944', '1945',
            '1946', '1947', '1948',
            '1949', '1950', '1951',
            '1952', '1953', '1954',
            '1955', '1956', '1957',
            '1958', '1959', '1960',
            '1961', '1962', '1963',
            '1964', '1965', '1966',
            '1967', '1968', '1969',
            '1970', '1971', '1972',
            '1973', '1974', '1975',
            '1976', '1977', '1978',
            '1979', '1980', '1981',
            '1982', '1983', '1984',
            '1985', '1986', '1987',
            '1988', '1989', '1990',
            '1991', '1992', '1993',
            '1994', '1995', '1996',
            '1997', '1998', '1999',
            '2000', '2001', '2002',
            '2003', '2004', '2005',
            '2006', '2007', '2008',
            '2009', '2010', '2011',
            '2012', '2013', '2014',
            '2015', '2016', '2017',
            '2018', '2019', '2020',
            '2021', '2022', '2023',
        ],
        'month': '07',
        'time': '00:00',
        'area': [90, -180, -90, 180], # maxlat, minlon, minlat, maxlon
    },
    f'{DATADIR}t2m_global_july_1940-2023.nc'
)
2025-07-15 10:14:14,343 INFO [2024-09-26T00:00:00] Watch our [Forum](https://forum.ecmwf.int/) for Announcements, news and other discussed topics.
2025-07-15 10:14:14,592 INFO Request ID is d03b586b-7198-451c-a56b-25da8fbd041e
2025-07-15 10:14:14,695 INFO status has been updated to accepted
2025-07-15 10:14:23,417 INFO status has been updated to running
2025-07-15 10:14:28,561 INFO status has been updated to accepted
2025-07-15 10:14:36,230 INFO status has been updated to running
2025-07-15 10:14:47,836 INFO status has been updated to successful
'./data_dir/t2m_global_july_1940-2023.nc'

Inspect data#

Now that we have downloaded the data, we can inspect it. We have requested the data in NetCDF format. This is a commonly used format for array-oriented scientific data. To read and process this data we will make use of the Xarray library. Xarray is an open source project and Python package that makes working with labelled multi-dimensional arrays simple and efficient. We will read the data from our NetCDF file into an xarray.Dataset.

ds = xr.open_dataset(f'{DATADIR}t2m_global_july_1940-2023.nc')

Now we can query our newly created Xarray dataset… Let’s have a look at the ds.

ds
<xarray.Dataset> Size: 349MB
Dimensions:     (valid_time: 84, latitude: 721, longitude: 1440)
Coordinates:
    number      int64 8B ...
  * valid_time  (valid_time) datetime64[ns] 672B 1940-07-01 ... 2023-07-01
  * latitude    (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude   (longitude) float64 12kB -180.0 -179.8 -179.5 ... 179.5 179.8
    expver      (valid_time) <U4 1kB ...
Data variables:
    t2m         (valid_time, latitude, longitude) float32 349MB ...
Attributes:
    GRIB_centre:             ecmf
    GRIB_centreDescription:  European Centre for Medium-Range Weather Forecasts
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             European Centre for Medium-Range Weather Forecasts
    history:                 2025-07-15T08:14 GRIB to CDM+CF via cfgrib-0.9.1...

We see that the dataset has one variable called t2m, which stands for “2 metre temperature”, and three coordinates of longitude, latitude and time.

Before we have commented to choose the sub-region extraction option but selecting the whole region in the request to the CDSAPI, 'area': [90, -180, -90, 180], # maxlat, minlon, minlat, maxlon. If we had chosen the whole available region option then the longitudes in the Xarray Dataset would be on a [0, 360] grid instead of [-180, 180] and we should correct this performing some extra steps.

If your longitudes are on a [0, 360] grid you could do the following to bring the longitude coordinates to a [-180, 180] grid:

ds_180 = ds.assign_coords(
    longitude=(((ds.longitude + 180) % 360) - 180)
).sortby('longitude')

There is also an expver coordinate. More on this later.

Select the icons to the right of the table above to expand the attributes of the coordinates and data variables. What are the units of the temperature data?

While an Xarray dataset may contain multiple variables, an Xarray data array holds a single multi-dimensional variable and its coordinates. To make the processing of the t2m data easier, we convert it into an Xarray data array. We will call it da_tmp (a temporary data array) because we will transform the data in some ways.

da_tmp = ds['t2m']

Let’s view this data:

da_tmp
<xarray.DataArray 't2m' (valid_time: 84, latitude: 721, longitude: 1440)> Size: 349MB
[87212160 values with dtype=float32]
Coordinates:
    number      int64 8B ...
  * valid_time  (valid_time) datetime64[ns] 672B 1940-07-01 ... 2023-07-01
  * latitude    (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude   (longitude) float64 12kB -180.0 -179.8 -179.5 ... 179.5 179.8
    expver      (valid_time) <U4 1kB ...
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            avgua
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                2 metre temperature
    units:                                    K
    standard_name:                            unknown
    GRIB_surface:                             0.0

From the result of the cell above you can see that now we have a xarray.DataArray.

Change temperature units from Kelvin to Celsius#

Notice that the ERA5 temperature data are in units of Kelvin, the base unit for temperature in the International System of Units (SI). If you want to convert the values from Kelvin to degrees Celsius, you have to subtract 273.15.

# Subtract 273.15 to convert to Celsius
da_celsius = da_tmp - 273.15
# Update the attributes of the DataArray to add Celsius units
da_celsius = da_celsius.assign_attrs(da_tmp.attrs)
da_celsius.attrs['units'] = '°C'

Data to be used#

The da_celsius data array will be used in the rest of the surface temperature exercise. Let’s check what we have:

da_celsius
<xarray.DataArray 't2m' (valid_time: 84, latitude: 721, longitude: 1440)> Size: 349MB
array([[[  0.85317993,   0.85317993,   0.85317993, ...,   0.85317993,
           0.85317993,   0.85317993],
        [  0.8590393 ,   0.8590393 ,   0.8590393 , ...,   0.8590393 ,
           0.8590393 ,   0.8590393 ],
        [  0.85317993,   0.85317993,   0.85317993, ...,   0.85317993,
           0.85317993,   0.85317993],
        ...,
        [-55.322586  , -55.320633  , -55.320633  , ..., -55.326492  ,
         -55.326492  , -55.32454   ],
        [-55.234695  , -55.234695  , -55.232742  , ..., -55.23665   ,
         -55.23665   , -55.234695  ],
        [-54.98079   , -54.98079   , -54.98079   , ..., -54.98079   ,
         -54.98079   , -54.98079   ]],

       [[  0.898468  ,   0.898468  ,   0.898468  , ...,   0.898468  ,
           0.898468  ,   0.898468  ],
        [  0.89456177,   0.89456177,   0.89456177, ...,   0.89456177,
           0.89456177,   0.89456177],
        [  0.8808899 ,   0.8808899 ,   0.8808899 , ...,   0.882843  ,
           0.8808899 ,   0.8808899 ],
...
         -56.495087  , -56.49118   ],
        [-56.72751   , -56.725555  , -56.723602  , ..., -56.733368  ,
         -56.731415  , -56.72946   ],
        [-56.506805  , -56.506805  , -56.506805  , ..., -56.506805  ,
         -56.506805  , -56.506805  ]],

       [[  0.7371521 ,   0.7371521 ,   0.7371521 , ...,   0.7371521 ,
           0.7371521 ,   0.7371521 ],
        [  0.7312927 ,   0.7312927 ,   0.7312927 , ...,   0.7312927 ,
           0.7312927 ,   0.7312927 ],
        [  0.7215271 ,   0.7215271 ,   0.7215271 , ...,   0.7215271 ,
           0.7215271 ,   0.7215271 ],
        ...,
        [-53.130035  , -53.12613   , -53.122223  , ..., -53.1398    ,
         -53.135895  , -53.13199   ],
        [-53.23941   , -53.237457  , -53.235504  , ..., -53.243317  ,
         -53.243317  , -53.241364  ],
        [-53.294098  , -53.294098  , -53.294098  , ..., -53.294098  ,
         -53.294098  , -53.294098  ]]],
      shape=(84, 721, 1440), dtype=float32)
Coordinates:
    number      int64 8B ...
  * valid_time  (valid_time) datetime64[ns] 672B 1940-07-01 ... 2023-07-01
  * latitude    (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude   (longitude) float64 12kB -180.0 -179.8 -179.5 ... 179.5 179.8
    expver      (valid_time) <U4 1kB ...
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            avgua
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                2 metre temperature
    units:                                    °C
    standard_name:                            unknown
    GRIB_surface:                             0.0

Now we can see the updated values in Celsius and the units attribute updated accordingly.

Plotting one timestep#

Just to check what we have so far, let’s plot a map of 2m temperature for the first (July 1940) and the last (July 2023) timesteps. We will plot these maps using the convenience method plot available for xarray.DataArray. This allows the creation of simple plots using one line of code. Also, with the xarray method sel(), you can select a data array based on coordinate labels.

da_celsius.sel(valid_time='1940').plot()
<matplotlib.collections.QuadMesh at 0x148623050>
../../_images/0c66955def6477e282e846abea5f4458edbcae98d5c5c27bd7a710143dc38b53.png
da_celsius.sel(valid_time='2023').plot()
<matplotlib.collections.QuadMesh at 0x15b026900>
../../_images/780b88152cb0e368d0e7050bc86954ef383cd41f0acb1afe5e321f18df046caa.png

Calculate a surface temperature climatology: reference period 1991-2020#

Standard reference periods and climatologies#

Anthropogenic activities and natural variations from years to decades shape the Earth’s climate. In order to evaluate anomalous conditions for a specific month or year, the World Meteorological Organization (WMO) defines standard reference periods used to create climatologies, also known as climate normals. Climatologies can be considered as the typical climate for the period they are based on.

Until 2020, the most current and widely used standard reference period was the 30-year range of 1981-2010. With the start of 2021, the WMO recommended updating the climate normal reference period to the range 1991-2020.

First, let us calculate the temperature climatology for July during the reference period 1991-2020. For this, we will select the reference period and then average along the time dimension. The resulting object contains the average July mean surface air temperature on each grid point.

t2m_ref_per = da_celsius.sel(valid_time=slice('1991-01-01', '2020-12-31')).mean(dim='valid_time')

If we have a look at this data object we will see now we have only two coordinates, latitude and longitude.

t2m_ref_per
<xarray.DataArray 't2m' (latitude: 721, longitude: 1440)> Size: 4MB
array([[  0.7482117 ,   0.7482117 ,   0.7482117 , ...,   0.7482117 ,
          0.7482117 ,   0.7482117 ],
       [  0.7403341 ,   0.7403992 ,   0.74046427, ...,   0.7397481 ,
          0.74000853,   0.74020386],
       [  0.72770387,   0.72783405,   0.7280294 , ...,   0.726532  ,
          0.72692263,   0.7273783 ],
       ...,
       [-55.417484  , -55.414036  , -55.410713  , ..., -55.427708  ,
        -55.42432   , -55.42087   ],
       [-55.408566  , -55.406937  , -55.405376  , ..., -55.413837  ,
        -55.412014  , -55.41039   ],
       [-55.04561   , -55.04561   , -55.04561   , ..., -55.04561   ,
        -55.04561   , -55.04561   ]], shape=(721, 1440), dtype=float32)
Coordinates:
    number     int64 8B ...
  * latitude   (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude  (longitude) float64 12kB -180.0 -179.8 -179.5 ... 179.5 179.8

We can also make a quick plot to have an exploratory view of this new xarray.DataArray:

t2m_ref_per.plot()
<matplotlib.collections.QuadMesh at 0x16a0a8680>
../../_images/c8003d85fa5cc438eb3058a1de64e5800a3f6c60df3a6a4d971b0476d47fe20b.png

Visualise surface temperature anomalies#

The next step is now to calculate the anomaly for July 2023 with respect to the climatology (1991-2020). The term anomaly refers to the deviation of a value from the long-term average. Positive or negative anomalies indicate that the average temperatures of a particular month were respectively warmer or cooler than the reference value for the same month.

Let us calculate the temperature anomaly for the year 2023. In a first step, we select the average near-surface temperature values for the year 2023 from the xarray.DataArray object da_celsius. As commented before, with the xarray method sel(), you can select a data array based on coordinate labels. The coordinate label of interest is year='2023'. As we are not doing, for instance, an aggregation operation like the mean over a dimension we will have the time dimension in the resulting data array. To remove the time dimension we can use the squeeze method.

t2m_july2023 = da_celsius.sel(valid_time='2023')
t2m_july2023
<xarray.DataArray 't2m' (valid_time: 1, latitude: 721, longitude: 1440)> Size: 4MB
array([[[  0.7371521,   0.7371521,   0.7371521, ...,   0.7371521,
           0.7371521,   0.7371521],
        [  0.7312927,   0.7312927,   0.7312927, ...,   0.7312927,
           0.7312927,   0.7312927],
        [  0.7215271,   0.7215271,   0.7215271, ...,   0.7215271,
           0.7215271,   0.7215271],
        ...,
        [-53.130035 , -53.12613  , -53.122223 , ..., -53.1398   ,
         -53.135895 , -53.13199  ],
        [-53.23941  , -53.237457 , -53.235504 , ..., -53.243317 ,
         -53.243317 , -53.241364 ],
        [-53.294098 , -53.294098 , -53.294098 , ..., -53.294098 ,
         -53.294098 , -53.294098 ]]], shape=(1, 721, 1440), dtype=float32)
Coordinates:
    number      int64 8B ...
  * valid_time  (valid_time) datetime64[ns] 8B 2023-07-01
  * latitude    (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude   (longitude) float64 12kB -180.0 -179.8 -179.5 ... 179.5 179.8
    expver      (valid_time) <U4 16B ...
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            avgua
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                2 metre temperature
    units:                                    °C
    standard_name:                            unknown
    GRIB_surface:                             0.0
t2m_july2023 = t2m_july2023.squeeze('valid_time')
t2m_july2023
<xarray.DataArray 't2m' (latitude: 721, longitude: 1440)> Size: 4MB
array([[  0.7371521,   0.7371521,   0.7371521, ...,   0.7371521,
          0.7371521,   0.7371521],
       [  0.7312927,   0.7312927,   0.7312927, ...,   0.7312927,
          0.7312927,   0.7312927],
       [  0.7215271,   0.7215271,   0.7215271, ...,   0.7215271,
          0.7215271,   0.7215271],
       ...,
       [-53.130035 , -53.12613  , -53.122223 , ..., -53.1398   ,
        -53.135895 , -53.13199  ],
       [-53.23941  , -53.237457 , -53.235504 , ..., -53.243317 ,
        -53.243317 , -53.241364 ],
       [-53.294098 , -53.294098 , -53.294098 , ..., -53.294098 ,
        -53.294098 , -53.294098 ]], shape=(721, 1440), dtype=float32)
Coordinates:
    number      int64 8B ...
    valid_time  datetime64[ns] 8B 2023-07-01
  * latitude    (latitude) float64 6kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0
  * longitude   (longitude) float64 12kB -180.0 -179.8 -179.5 ... 179.5 179.8
    expver      <U4 16B ...
Attributes: (12/32)
    GRIB_paramId:                             167
    GRIB_dataType:                            an
    GRIB_numberOfPoints:                      1038240
    GRIB_typeOfLevel:                         surface
    GRIB_stepUnits:                           1
    GRIB_stepType:                            avgua
    ...                                       ...
    GRIB_totalNumber:                         0
    GRIB_units:                               K
    long_name:                                2 metre temperature
    units:                                    °C
    standard_name:                            unknown
    GRIB_surface:                             0.0

The anomaly will be the difference between t2m_july2023 and t2m_ref_per. A positive value means July 2023 is above the expected mean:

anom = t2m_july2023 - t2m_ref_per

The previous operation results in the anomaly on each longitude and latitude location stored in the anom data array. We can plot this in a map to check where the anomaly was positive (July 2023 warmer than the climatology) or negative (July 2023 colder than the climatology). This time we will create the plot using the matplotlib and cartopy libraries.

The code below is briefly commented in the code cell. The key part is plt.pcolormesh, which is used to plot the 2D data array on each (longitude, latitude) grid cell. The colours are defined using the RdBr_r colormap in a range between -12 and 12 (ºC). This is a diverging colormap that goes from red to blue but in reverse order, this is why it has the _r suffix. So, in this specific case, -12 degrees is shown in dark blue, 12 in dark red and 0, equal to the climate normal, in white.

# create the figure panel and the map using the Cartopy PlateCarree projection
fig, ax = plt.subplots(1, 1, figsize = (16, 8), subplot_kw={'projection': ccrs.PlateCarree()})

# Plot the data
im = plt.pcolormesh(
    anom.longitude, 
    anom.latitude, 
    anom, 
    cmap='RdBu_r',
    vmin=-12, 
    vmax=12
) 

# Set the figure title, add lat/lon grid and coastlines
ax.set_title('Temperature anomaly for July 2023 (with respect to 1991-2020 July mean)', fontsize=16)
ax.gridlines(
    draw_labels=True, 
    linewidth=1, 
    color='gray', 
    alpha=0.5, 
    linestyle='--'
) 
ax.coastlines(color='black')

# Specify the colorbar and set a label for the colorbar
cbar = plt.colorbar(im, fraction=0.05, pad=0.04)
cbar.set_label('Temperature anomaly') 

# Show or save the figure, uncomment the line/s in case of need
#fig.show() # not needed in a notebook inline context
#fig.savefig('near_sfc_t2m_anomaly_july2023.png') # not needed in a notebook inline context
../../_images/93a63a43bc6af9702ba022a351ffbfe4625166ea753e7056cafe1a07244574e7.png