CAMS Global Reanalysis Practical#

Run the tutorial via free cloud platforms: binder kaggle colab

Learning objectives#

  • In this practical exercise you will learn how to download CAMS Global Reanalysis (EAC4) data programmatically using the Application Programming Interface (API) of the Atmosphere Data Store (ADS).

  • You will also learn how to read the data into a Python object and explore its characteristics, including data dimensions, units, etc.

  • Having explored the data, you will calculate the June climatology for AOD at 550nm for the period 2003 to 2024.

  • Finally, you will explore the famous Godzilla dust storm of June 2020, and analyse this event in the context of the 2003 to 2024 climatology.

This practical session is based on a real event reported in the following article https://atmosphere.copernicus.eu/tracking-massive-dust-cloud-africa-america

Initial setup#

Before we begin we must prepare our environment. This includes installing the Application Programming Interface (API) of the Atmosphere Data Store (ADS), intalling any other packages not already installed, setting up our ADS API credentials and importing the various Python libraries that we will need.

# Ensure that the cdsapi package is installed
!pip install -q cdsapi

# If you are running this notebook in Colab, uncomment the line below and run this cell.
#!pip install cartopy

Add your ADS API credentials#

To set up your ADS API credentials, please login/register on the ADS, then you will see your unique API key here.

You can add this API key to your current session by replacing ######### in the code cell below with your API key.

import os
os.environ['CDSAPI_URL'] = 'https://ads.atmosphere.copernicus.eu/api'
os.environ['CDSAPI_KEY'] = '###########################################'

Import libraries#

# CDS API
import cdsapi

# Libraries for working with multidimensional arrays
import numpy as np
import xarray as xr

# Libraries for plotting and visualising data
import matplotlib.pyplot as plt
import cartopy
import cartopy.crs as ccrs

# Disable warnings for data download via API
import urllib3 
urllib3.disable_warnings()

Here we specify a data directory in which we will download our data and all output files that we will generate:

DATADIR = '.'

Explore and download data#

Visit the download form for the CAMS global reanalysis data https://ads.atmosphere.copernicus.eu/datasets/cams-global-reanalysis-eac4?tab=download. View the parameters in the API script in the following cell and select the corresponding options.

At the end of the download form, select “Show API request”. This will reveal a block of code, which should be identical to the code cell below.

Please remember to accept the terms and conditions at the bottom of the download form.

Download 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.

dataset = "cams-global-reanalysis-eac4-monthly"
request = {
    'variable': ['total_aerosol_optical_depth_550nm', 'dust_aerosol_optical_depth_550nm', 'organic_matter_aerosol_optical_depth_550nm'],
    'year': ['2003', '2004', '2005', '2006', '2007', '2008', '2009', '2010', '2011', '2012', '2013', '2014', '2015', '2016', '2017', '2018', '2019', '2020', '2021', '2022', '2023', '2024'],
    'month': ['01','02','03','04','05','06','07','08','09','10','11','12'],
    'product_type': ['monthly_mean'],
    'data_format': 'netcdf'
}

client = cdsapi.Client()
client.retrieve(dataset, request).download(f'{DATADIR}/aod-550nm-global-reanalysis.nc')

Inspect data#

# Path to the downloaded file
netcdf_file = f'{DATADIR}/aod-550nm-global-reanalysis.nc'

# Create Xarray Dataset
ds = xr.open_dataset(netcdf_file)

# view the dataset
ds

month = 6
mstr = 'june'

ds_sel = ds.isel(valid_time=(ds.valid_time.dt.month == month))

# create xarray data array object (single variable)
da = ds_sel['aod550']
da

Calculate climatology for June 2003 to 2024#

climatology = da.mean(dim="valid_time")

climatology

Calculate anomaly for June 2020#

june2020 = da.sel(valid_time = '2020-06-01T00:00:00.000000000') - climatology

Plot global maps to view the June climatology and anomaly for 2020#

Define a function to plot a map#

As we will create several maps, creating a function to do this reduces the amount of boilerplate code, i.e. repeated copies of the same block of code.

def create_figure(): # Define a function
    """This function creates a map"""
    fig = plt.figure(figsize=(20, 10)) # Create a figure
    ax = plt.axes(projection=ccrs.PlateCarree()) # Create axes and define the map projection
    ax.coastlines() # Add coastlines
    gl = ax.gridlines(draw_labels=True, linestyle='--') # Add gridlines
    return fig, ax # Define what the function returns (the figure and axes objects)

Plot a map of the June climatology for dust AOD for the period 2003 to 2024#

_, ax = create_figure() # Call the function to create a map figure
climatology.plot.pcolormesh(ax=ax, x='longitude', y='latitude', add_colorbar=True, cmap='YlOrRd')
plt.title('AOD at 550nm, June climatology for period 2003 to 2024', fontsize=12)
plt.savefig('aod_june_2003-2024_climatology.png')

Plot a map of the dust AOD anomaly for June 2020#

_, ax = create_figure() # Call the function to create a map figure
june2020.plot.pcolormesh(ax=ax, x='longitude', y='latitude', add_colorbar=True, cmap='RdYlBu_r')
plt.title('AOD at 550nm anomaly for June 2020', fontsize=12)
plt.savefig('aod_june2020_anomaly.png')

The anomaly map puts events into context. In this case we can see the Godzilla dust storm produced higher than average dust in the atmosphere in June 2020, with respect to the 2003 to 2023 climatology. For more information about this event, visit https://atmosphere.copernicus.eu/tracking-massive-dust-cloud-africa-america

ragged_point_lat = 13.17
ragged_point_lon = -59.43

ragged_point_da = ds['aod550'].sel(latitude = ragged_point_lat, longitude = ragged_point_lon, method='nearest')
ragged_point_om = ds['omaod550'].sel(latitude = ragged_point_lat, longitude = ragged_point_lon, method='nearest')
ragged_point_du = ds['duaod550'].sel(latitude = ragged_point_lat, longitude = ragged_point_lon, method='nearest')

fig = plt.figure(figsize=(10, 5))
ragged_point_da.plot(marker='o')
ragged_point_om.plot(marker='o', color='magenta')
ragged_point_du.plot(marker='o', color='orange')
plt.suptitle("AOD at 550nm, Ragged Point")
plt.grid(True)

payerne_lat = 46.81
payerne_lon = 6.94
payerne_da = ds['aod550'].sel(latitude = payerne_lat, longitude = payerne_lon, method='nearest')
payerne_om = ds['omaod550'].sel(latitude = payerne_lat, longitude = payerne_lon, method='nearest')
payerne_du = ds['duaod550'].sel(latitude = payerne_lat, longitude = payerne_lon, method='nearest')

fig = plt.figure(figsize=(10, 5))
payerne_da.plot(marker='o')
payerne_om.plot(marker='o', color='magenta')
payerne_du.plot(marker='o', color='orange')
plt.suptitle("AOD at 550nm, Payerne")
plt.grid(True)

Suggested further steps#

Carry-out the same steps, but for organic matter AOD and with a focus on Africa (add 'area': [40, -20, -40, 55,] in the data retreive script, or select North 40, East 50, South -40, West -20 in the Restricted area field of the data download form.)

Were the wildfire emissions in June 2020 in Africa greater or less than the 2003 to 2024 climatological average?

How do 2025 monthly mean AOD values compare with the climatology?