Create Maps from CAMS Data¶

This tutorial demonstrates how to visualise data from the Copernicus Atmosphere Monitoring Service (CAMS) in the form of two dimensional maps.

Run the tutorial via free cloud platforms:

Install and import packages¶

!pip install cdsapi

Requirement already satisfied: cdsapi in c:\users\cxcs\anaconda3\lib\site-packages (0.5.1)
Requirement already satisfied: requests>=2.5.0 in c:\users\cxcs\anaconda3\lib\site-packages (from cdsapi) (2.26.0)
Requirement already satisfied: tqdm in c:\users\cxcs\anaconda3\lib\site-packages (from cdsapi) (4.62.3)
Requirement already satisfied: charset-normalizer~=2.0.0 in c:\users\cxcs\anaconda3\lib\site-packages (from requests>=2.5.0->cdsapi) (2.0.4)
Requirement already satisfied: urllib3<1.27,>=1.21.1 in c:\users\cxcs\anaconda3\lib\site-packages (from requests>=2.5.0->cdsapi) (1.26.7)
Requirement already satisfied: certifi>=2017.4.17 in c:\users\cxcs\anaconda3\lib\site-packages (from requests>=2.5.0->cdsapi) (2021.10.8)
Requirement already satisfied: idna<4,>=2.5 in c:\users\cxcs\anaconda3\lib\site-packages (from requests>=2.5.0->cdsapi) (3.2)
Requirement already satisfied: colorama in c:\users\cxcs\anaconda3\lib\site-packages (from tqdm->cdsapi) (0.4.4)

# CDS API
import cdsapi

# Library to extract data
from zipfile import ZipFile

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

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

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

Access data¶

Copy your API key into the code cell below, replacing ####### with your key. (Remember, to access data from the ADS, you will need first to register/login https://ads.atmosphere.copernicus.eu and obtain an API key from https://ads.atmosphere.copernicus.eu/api-how-to.)

URL = 'https://ads.atmosphere.copernicus.eu/api/v2'

# Replace the hashtags with your key:
KEY = '########################################'

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

DATADIR = './'

For our first plotting example, we will use CAMS Global Atmospheric Composition Forecast data. The code below shows the subset characteristics that we will extract from this dataset for the purpose of this tutorial as an API request.

Note

Before running this code, ensure that you have accepted the terms and conditions. This is something you only need to do once for each CAMS dataset. You will find the option to do this by selecting the dataset in the ADS, then scrolling to the end of the Download data tab.

c = cdsapi.Client(url=URL, key=KEY)
c.retrieve(
    'cams-global-atmospheric-composition-forecasts',
    {
        'variable': [
            'dust_aerosol_optical_depth_550nm', 'organic_matter_aerosol_optical_depth_550nm', 'total_aerosol_optical_depth_550nm',
        ],
        'date': '2021-08-01/2021-08-01',
        'time': '12:00',
        'leadtime_hour': '0',
        'type': 'forecast',
        'format': 'netcdf_zip',
    },
    f'{DATADIR}/2021-08-01_AOD.zip')

2022-08-23 23:13:56,731 INFO Welcome to the CDS
2022-08-23 23:13:56,736 INFO Sending request to https://ads.atmosphere.copernicus.eu/api/v2/resources/cams-global-atmospheric-composition-forecasts
2022-08-23 23:13:56,813 INFO Request is queued
2022-08-23 23:13:57,880 INFO Request is running
2022-08-23 23:14:01,797 INFO Request is completed
2022-08-23 23:14:01,799 INFO Downloading https://download-0002-ads-clone.copernicus-climate.eu/cache-compute-0002/cache/data8/adaptor.mars_constrained.external-1661292784.7694805-11001-4-cbd2c200-e549-4301-b6c5-af2915739396.zip to ./DATA/2021-08-01_AOD.zip (2.3M)
2022-08-23 23:14:02,601 INFO Download rate 2.9M/s                                                                      

Result(content_length=2442568,content_type=application/zip,location=https://download-0002-ads-clone.copernicus-climate.eu/cache-compute-0002/cache/data8/adaptor.mars_constrained.external-1661292784.7694805-11001-4-cbd2c200-e549-4301-b6c5-af2915739396.zip)

Plot data¶

There are several ways we can visualise this data in a plot. To create a static plot of one variable for one timestep, we will first create an xarray Data Array which includes only one variable. In this case, we will focus on Total AOD at 550 nm, or aod550 (this is how the variable is named in our xarray Dataset).

# Create Xarray Data Array
da = ds['aod550']
da

<xarray.DataArray 'aod550' (time: 1, latitude: 451, longitude: 900)>
[405900 values with dtype=float32]
Coordinates:
  * longitude  (longitude) float32 0.0 0.4 0.8 1.2 ... 358.4 358.8 359.2 359.6
  * latitude   (latitude) float32 90.0 89.6 89.2 88.8 ... -89.2 -89.6 -90.0
  * time       (time) datetime64[ns] 2021-08-01T12:00:00
Attributes:
    units:      ~
    long_name:  Total Aerosol Optical Depth at 550nm

Simple plot using xarray¶

A very simple way to visualise this data is to use xarray’s own plotting functionality. This can be done with the .plot() method:

da.plot()

<matplotlib.collections.QuadMesh at 0x12065458400>

This figure would be easier to interpret with some modifications. These may include adjusting the colourscale, or adding annotations such as coastlines. While xarray’s plotting functionality can be further customised, it is essentially a thin wrapper around the popular matplotlib library. The example below shows how to generate a more customised visualisation using both matplotlib (for visualisations) and cartopy (for geospatial data processing).

Customised map with matplotlib and cartopy¶

# create the figure panel and specify size
fig = plt.figure(figsize=(15, 10))

# create the map using the cartopy PlateCarree projection
ax = plt.subplot(1,1,1, projection=ccrs.PlateCarree())

# Add lat/lon grid
ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--')

# Set figure title
ax.set_title('AOD at 550nm, 2021-08-01T12:00', fontsize=12) 

# Plot the data
im = plt.pcolormesh(da.longitude, da.latitude, da[0,:,:], cmap='YlOrRd', vmin=0, vmax=2) 

# Add coastlines
ax.coastlines(color='black') 

# Specify the colourbar, including fraction of original axes to use for colorbar, 
# and fraction of original axes between colorbar and new image axes
cbar = plt.colorbar(im, fraction=0.025, pad=0.05) 

# Define the colourbar label
cbar.set_label('AOD at 550nm') 

# Save the figure
fig.savefig(f'{DATADIR}/CAMS_global_forecast_AOD_2021-08-01.png')

The colours used in the above figure are from one of many built-in colourmaps used in matplotlib. We can also create our own colourmap. The code below shows how to create the same colourmap as used in many CAMS products, such as the CAMS aerosol forecasts.

Create custom colourmap¶

# This code creates a colourmap similar to that used in CAMS products.

# Import the Matplotlib class 'ListedColormap' to generate a colormap object from a given list of colours
from matplotlib.colors import ListedColormap

# Array containing 15 colours ranging from white, grey, yellow and orange to red.
# Each row in the array contains three numbers for red, green and blue respectively.
# The numbers represent intensity, from 0 (black) to 256 (white).
matrix = np.array([[256, 256, 256],
                   [210, 214, 234],
                   [167, 174, 214],
                   [135, 145, 190],
                   [162, 167, 144],
                   [189, 188, 101],
                   [215, 209, 57],
                   [242, 230, 11],
                   [243, 197, 10],
                   [245, 164, 8],
                   [247, 131, 6],
                   [248, 98, 5],
                   [250, 65, 3],
                   [252, 32, 1],
                   [254, 0, 0]])
n = 17 # Multiplication number

cams = np.ones((253, 4)) # Initial empty colormap, to be filled by the colours in 'matrix'.

# This loop fills in the empty 'cams' colormap with each of the 15 colours in 'matrix'
# multiplied by 'n'. Each colour is divided by 256 to normalise from 0 (black) to 1 (white). 
for i in range(matrix.shape[0]):
    cams[(i*n):((i+1)*n),:] = np.array([matrix[i,0]/256, matrix[i,1]/256, matrix[i,2]/256, 1])

# The final colormap is given by 'camscmp', which uses the Matplotlib class 'ListedColormap(Colormap)'
# to generate a colormap object from the list of colours provided by 'cams'.
camscmp = ListedColormap(cams)

We now apply this custom colourmap to the same code above to generate the same figure, but with the CAMS colourmap. Note the only difference between this and the above visualisation code (apart from removing the comments) is to replace in the plt.pcolormesh() function, cmap='YlOrRd' (the standard yellow-orange-red matplotlib colourmap) with cmap=camscmp (our custom-made CAMS colourmap).

fig = plt.figure(figsize=(15, 10))
ax = plt.subplot(1,1,1, projection=ccrs.PlateCarree())
ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--') 
ax.set_title('AOD at 550nm, 2021-08-01T12:00', fontsize=12) 
im = plt.pcolormesh(da.longitude, da.latitude, da[0,:,:], cmap=camscmp, vmin=0, vmax=2) 
ax.coastlines(color='black') 
cbar = plt.colorbar(im, fraction=0.025, pad=0.05)
cbar.set_label('AOD at 550nm') 
fig.savefig(f'{DATADIR}/CAMS_global_forecast_AOD_2021-08-01_CAMS-cmap.png')

Create subplots¶

The examples above demonstrate how to create a single map. However, we may want to create multiple maps or plots in the same figure, for example to compare maps of different variables side by side. The example below shows one way to do this.

# Here we extract the names of the variables from the xarray Dataset

variables = list(ds.keys())
variables

['omaod550', 'duaod550', 'aod550']

# We would like to change the order in which variables are visualised, 
# beginning with Total AOD (last in the list, becomes first), 
# then breaking this down into Dust AOD (order remains the same), 
# and Organic Matter AOD (first in the list, becomes last).

variables[0], variables[2] = variables[2], variables[0]

Having extracted the names of variables, inserted them into the right order, we can now loop through them to produce a single figure with subplots containing maps of each variable side-by-side. Given that most AOD concentration in this data is in the Northern Hemisphere, we will also change the projection from PlateCarree to Orthographic, centred at 90 degrees north.

fig, axs = plt.subplots(1, 3, figsize = (20, 10), subplot_kw={'projection': ccrs.Orthographic(central_latitude=90)})

for i in range(3):
    da = ds[variables[i]]
    axs[i].gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--') # Add lat/lon grid
    axs[i].set_title(f'{variables[i]}', fontsize=12) # Set figure title
    im = axs[i].pcolormesh(da.longitude, da.latitude, da[0,:,:],
                            transform = ccrs.PlateCarree(), cmap=camscmp, vmin=0, vmax=2)
    axs[i].coastlines(color='black') # Add coastlines
    cbar = fig.colorbar(im, ax=axs[i], fraction=0.046, pad=0.04) # Specify the colourbar
    cbar.set_label(variables[i]) # Define the colourbar label

plt.show() # Display the figure
fig.savefig(f'{DATADIR}/AOD_NHem.png') # Save the figure

We can see that the high values of AOD seem to be mainly due to organic matter in North America and Siberia, which saw many wildfire activity in this period, while further south we can see a dust contribution.

CAMS Training

Create Maps from CAMS Data

Contents