Training school and workshop on dust

CAMS global atmospheric composition forecasts - Dust Aerosol Optical Depth¶

About¶

This notebooks provides an introduction to the CAMS global atmospheric composition forecasts and shows you how the variable Dust Aerosol Optical Depth can be used to monitor and predict dust events.

CAMS produces global forecasts for atmospheric composition for the next five days twice a day. The forecasts consist of more than 50 chemical species (e.g. ozone, nitrogen dioxide, carbon dioxide) and seven different types of aerosol (desert dust, sea salt, organic matter, black carbon, sulphate, nitrate and ammonium aerosol).

The initial conditions of each forecast are obtained by combining a previous forecast with current satellite observations through a process called data assimilation. This best estimate of the state of the atmosphere at the initial forecast time step, called the analysis, provides a globally complete and consistent dataset allowing for estimates at locations where observation data coverage is low or for atmospheric pollutants for which no direct observations are available.

The notebook examines the Saharan Dust event which occured over Europe at the beginning of February 2021.

Basic Facts¶

Spatial resolution: 0.4° x 0.4°
Spatial coverage: Global
Temporal resolution: 1-hourly up to leadtime hour 120
Temporal coverage: since 2015 to present
Data format: GRIB or zipped NetCDF

How to access the data¶

CAMS global atmospheric composition forecasts are available for download via the Copernicus Atmosphere Data Store (ADS). You will need to create an ADS account here.

Data from the ADS can be downloaded in two ways:

manually via the ADS web interface
programmatically with a Python package called cdsapi (more information)

Module outline¶

Optional: Retrieve CAMS global atmospheric forecasts programmatically
1 - Load and browse dust aerosol optical depth (duaod) at 550nm of the CAMS global near-real-time forecast
2 - Retrieve the data variable dust AOD at 550nm as xarray.DataArray
3 - Visualize dust AOD at 550nm
4 - Create a geographical subset for Europe
5 - Animate dust AOD at 550nm changes over time

Load required libraries¶

In [1]:

import xarray as xr

import matplotlib.pyplot as plt
import matplotlib.colors
from matplotlib.cm import get_cmap
from matplotlib import animation
from matplotlib.axes import Axes

import cartopy.crs as ccrs
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import cartopy.feature as cfeature
from cartopy.mpl.geoaxes import GeoAxes
GeoAxes._pcolormesh_patched = Axes.pcolormesh

from IPython.display import HTML

import warnings
warnings.simplefilter(action = "ignore", category = RuntimeWarning)

Load helper functions¶

In [2]:

%run ../functions.ipynb

Optional: Retrieve CAMS global atmospheric composition forecasts programmatically¶

The CDS Application Program Interface (CDS API) is a Python library which allows you to access data from the ADS programmatically. In order to use the CDS API, follow the steps below:

Self-register at the ADS registration page (if you do not have an account yet)
Login to the ADS portal and go to the api-how-to page
Copy the CDS API key displayed in the black terminal window and replace the ###### of the KEY variable below with your individual CDS API key

Note: You find your CDS API key displayed in the black terminal box under the section Install the CDS API key. If you do not see a URL or key appear in the black terminal box, please refresh your browser tab.

In [3]:

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

The next step is then to request the data with a so called API request. Via the ADS web interface, you can select the data and at the end of the web interface, you can open the ADS request via Show API request. Copy the request displayed there in the cell below. Once you execute the cell, the download of the data starts automatically.

In [ ]:

import cdsapi

c = cdsapi.Client(url=URL, key=KEY)

c.retrieve(
    'cams-global-atmospheric-composition-forecasts',
    {
        'variable': 'dust_aerosol_optical_depth_550nm',
        'date': '2021-02-04/2021-02-04',
        'time': '12:00',
        'leadtime_hour': [
            '0', '12', '15',
            '18', '21', '24',
            '27', '3', '30',
            '33', '36', '39',
            '42', '45', '48',
            '51', '54', '57',
            '6', '60', '63',
            '66', '69', '72',
            '75', '78', '81',
            '84', '87', '9',
            '90',
        ],
        'type': 'forecast',
        'format': 'netcdf_zip',
    },
    './20210204_duaod.netcdf_zip')

CAMS global atmospheric composition forecasts can be retrieved either in GRIB or in a zipped NetCDF. Above, we requested the data in a zipped NetCDF and for this reason, we have to unzip the file before we can open it. You can unzip zip archives in Python with the Python package zipfile and the function extractall(). You will see a new file called data.nc appearing in the same folder as this notebook. This is just for demonstration purposes.

In [5]:

import zipfile
with zipfile.ZipFile('../../../eodata/dust/part1/3_model/cams/global_forecast/20210204_duaod.netcdf_zip', 'r') as zip_ref:
    zip_ref.extractall('./')

1. Load and browse CAMS global atmospheric composition forecasts ¶

CAMS global near-real-time forecast data is available either in GRIB or netCDF. The data for the present example has been downloaded as netCDF.

You can use xarray's function xr.open_dataset() to open the netCDF file as xarray.Dataset.

In [3]:

file = xr.open_dataset('../../../eodata/dust/part1/3_model/cams/global_forecast/data.nc')
file

The data above has three dimensions (latitude, longitude, time) and one data variable:

duaod550: Dust Aerosol Optical Depth at 550nm

Let us inspect the coordinates of the file more in detail.

Below, you see that the data set consists of 31 time steps, ranging from 4 February 2021 12 UTC to 8 February 2021 06 UTC in a 3-hour timestep.

In [4]:

file.time

Out[4]:

<xarray.DataArray 'time' (time: 31)>
array(['2021-02-04T12:00:00.000000000', '2021-02-04T15:00:00.000000000',
       '2021-02-04T18:00:00.000000000', '2021-02-04T21:00:00.000000000',
       '2021-02-05T00:00:00.000000000', '2021-02-05T03:00:00.000000000',
       '2021-02-05T06:00:00.000000000', '2021-02-05T09:00:00.000000000',
       '2021-02-05T12:00:00.000000000', '2021-02-05T15:00:00.000000000',
       '2021-02-05T18:00:00.000000000', '2021-02-05T21:00:00.000000000',
       '2021-02-06T00:00:00.000000000', '2021-02-06T03:00:00.000000000',
       '2021-02-06T06:00:00.000000000', '2021-02-06T09:00:00.000000000',
       '2021-02-06T12:00:00.000000000', '2021-02-06T15:00:00.000000000',
       '2021-02-06T18:00:00.000000000', '2021-02-06T21:00:00.000000000',
       '2021-02-07T00:00:00.000000000', '2021-02-07T03:00:00.000000000',
       '2021-02-07T06:00:00.000000000', '2021-02-07T09:00:00.000000000',
       '2021-02-07T12:00:00.000000000', '2021-02-07T15:00:00.000000000',
       '2021-02-07T18:00:00.000000000', '2021-02-07T21:00:00.000000000',
       '2021-02-08T00:00:00.000000000', '2021-02-08T03:00:00.000000000',
       '2021-02-08T06:00:00.000000000'], dtype='datetime64[ns]')
Coordinates:
  * time     (time) datetime64[ns] 2021-02-04T12:00:00 ... 2021-02-08T06:00:00
Attributes:
    long_name:  time

xarray.DataArray

'time'

time: 31

2021-02-04T12:00:00 2021-02-04T15:00:00 ... 2021-02-08T06:00:00

array(['2021-02-04T12:00:00.000000000', '2021-02-04T15:00:00.000000000',
       '2021-02-04T18:00:00.000000000', '2021-02-04T21:00:00.000000000',
       '2021-02-05T00:00:00.000000000', '2021-02-05T03:00:00.000000000',
       '2021-02-05T06:00:00.000000000', '2021-02-05T09:00:00.000000000',
       '2021-02-05T12:00:00.000000000', '2021-02-05T15:00:00.000000000',
       '2021-02-05T18:00:00.000000000', '2021-02-05T21:00:00.000000000',
       '2021-02-06T00:00:00.000000000', '2021-02-06T03:00:00.000000000',
       '2021-02-06T06:00:00.000000000', '2021-02-06T09:00:00.000000000',
       '2021-02-06T12:00:00.000000000', '2021-02-06T15:00:00.000000000',
       '2021-02-06T18:00:00.000000000', '2021-02-06T21:00:00.000000000',
       '2021-02-07T00:00:00.000000000', '2021-02-07T03:00:00.000000000',
       '2021-02-07T06:00:00.000000000', '2021-02-07T09:00:00.000000000',
       '2021-02-07T12:00:00.000000000', '2021-02-07T15:00:00.000000000',
       '2021-02-07T18:00:00.000000000', '2021-02-07T21:00:00.000000000',
       '2021-02-08T00:00:00.000000000', '2021-02-08T03:00:00.000000000',
       '2021-02-08T06:00:00.000000000'], dtype='datetime64[ns]')

Coordinates: (1)

time

(time)

datetime64[ns]

2021-02-04T12:00:00 ... 2021-02-...

long_name :: time

array(['2021-02-04T12:00:00.000000000', '2021-02-04T15:00:00.000000000',
       '2021-02-04T18:00:00.000000000', '2021-02-04T21:00:00.000000000',
       '2021-02-05T00:00:00.000000000', '2021-02-05T03:00:00.000000000',
       '2021-02-05T06:00:00.000000000', '2021-02-05T09:00:00.000000000',
       '2021-02-05T12:00:00.000000000', '2021-02-05T15:00:00.000000000',
       '2021-02-05T18:00:00.000000000', '2021-02-05T21:00:00.000000000',
       '2021-02-06T00:00:00.000000000', '2021-02-06T03:00:00.000000000',
       '2021-02-06T06:00:00.000000000', '2021-02-06T09:00:00.000000000',
       '2021-02-06T12:00:00.000000000', '2021-02-06T15:00:00.000000000',
       '2021-02-06T18:00:00.000000000', '2021-02-06T21:00:00.000000000',
       '2021-02-07T00:00:00.000000000', '2021-02-07T03:00:00.000000000',
       '2021-02-07T06:00:00.000000000', '2021-02-07T09:00:00.000000000',
       '2021-02-07T12:00:00.000000000', '2021-02-07T15:00:00.000000000',
       '2021-02-07T18:00:00.000000000', '2021-02-07T21:00:00.000000000',
       '2021-02-08T00:00:00.000000000', '2021-02-08T03:00:00.000000000',
       '2021-02-08T06:00:00.000000000'], dtype='datetime64[ns]')

Attributes: (1)
long_name :
time

The latitude values have a 0.4 degrees resolution and have a global N-S coverage.

In [5]:

file.latitude

Out[5]:

<xarray.DataArray 'latitude' (latitude: 451)>
array([ 90. ,  89.6,  89.2, ..., -89.2, -89.6, -90. ], dtype=float32)
Coordinates:
  * latitude  (latitude) float32 90.0 89.6 89.2 88.8 ... -88.8 -89.2 -89.6 -90.0
Attributes:
    units:      degrees_north
    long_name:  latitude

The longitude values have a 0.4 degrees resolution as well, and are disseminated in a [0, 360] grid.

In [6]:

file.longitude

Out[6]:

<xarray.DataArray 'longitude' (longitude: 900)>
array([  0. ,   0.4,   0.8, ..., 358.8, 359.2, 359.6], dtype=float32)
Coordinates:
  * longitude  (longitude) float32 0.0 0.4 0.8 1.2 ... 358.4 358.8 359.2 359.6
Attributes:
    units:      degrees_east
    long_name:  longitude

Above, you see that the longitude variables are in the range of [0, 359.6]. Per default, ECMWF data are on a [0, 360] grid. Let us bring the longitude coordinates to a [-180, 180] grid. You can use the xarray function assign_coords() to assign new values to coordinates of a xarray data arra. The code below shifts your longitude coordinates from [0, 359.6] to [-180, 179.6].

In [7]:

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

2. Retrieve the variable Dust Aerosol Optical Depth at 550nm as data array ¶

Let us extract from the dataset above the data variable Dust Aerosol Optical Depth (AOD) at 550nm as xarray.DataArray with the name du_aod. You can load a data array from a xarray dataset by specifying the name of the variable (duaod550) in square brackets.

In [8]:

du_aod = file['duaod550']
du_aod

Out[8]:

<xarray.DataArray 'duaod550' (time: 31, latitude: 451, longitude: 900)>
[12582900 values with dtype=float32]
Coordinates:
  * longitude  (longitude) float32 -180.0 -179.6 -179.2 ... 178.8 179.2 179.6
  * latitude   (latitude) float32 90.0 89.6 89.2 88.8 ... -89.2 -89.6 -90.0
  * time       (time) datetime64[ns] 2021-02-04T12:00:00 ... 2021-02-08T06:00:00
Attributes:
    units:      ~
    long_name:  Dust Aerosol Optical Depth at 550nm

Above, you see that the variable du_aod has two attributes, units and long_name. Let us define variables for those attributes. The variables can be used later for visualizing the data.

In [9]:

long_name = du_aod.long_name
units = du_aod.units

Let us do the same for the coordinates longitude and latitude.

In [10]:

latitude = du_aod.latitude
longitude = du_aod.longitude

3. Visualize Dust Aerosol Optical Depth at 550nm¶

The next step is to visualize the dataset. You can use the function visualize_pcolormesh, which makes use of matploblib's function pcolormesh and the Cartopy library.

With ?visualize_pcolormesh you can open the function's docstring to see what keyword arguments are needed to prepare your plot.

In [11]:

?visualize_pcolormesh

Signature:
visualize_pcolormesh(
    data_array,
    longitude,
    latitude,
    projection,
    color_scale,
    unit,
    long_name,
    vmin,
    vmax,
    set_global=True,
    lonmin=-180,
    lonmax=180,
    latmin=-90,
    latmax=90,
)
Docstring:
Visualizes a xarray.DataArray with matplotlib's pcolormesh function.

Parameters:
    data_array(xarray.DataArray): xarray.DataArray holding the data values
    longitude(xarray.DataArray): xarray.DataArray holding the longitude values
    latitude(xarray.DataArray): xarray.DataArray holding the latitude values
    projection(str): a projection provided by the cartopy library, e.g. ccrs.PlateCarree()
    color_scale(str): string taken from matplotlib's color ramp reference
    unit(str): the unit of the parameter, taken from the NetCDF file if possible
    long_name(str): long name of the parameter, taken from the NetCDF file if possible
    vmin(int): minimum number on visualisation legend
    vmax(int): maximum number on visualisation legend
    set_global(boolean): optional kwarg, default is True
    lonmin,lonmax,latmin,latmax(float): optional kwarg, set geographic extent is set_global kwarg is set to 
                                        False
File:      /tmp/ipykernel_783/2163691773.py
Type:      function

You can make use of the variables we defined above:

units
long_name
latitude
longitude

Additionally, you can specify the color scale and minimum (vmin) and maxium (vmax) data values.

In [12]:

time_index = 10
visualize_pcolormesh(data_array=du_aod[time_index,:,:],
                     longitude=longitude,
                     latitude=latitude,
                     projection=ccrs.PlateCarree(),
                     color_scale='hot_r',
                     unit=units,
                     long_name=long_name + ' ' + str(du_aod[time_index,:,:].time.data)[0:19],
                     vmin=0, 
                     vmax=1.5)

Out[12]:

(<Figure size 1440x720 with 2 Axes>,
 <GeoAxesSubplot:title={'center':'Dust Aerosol Optical Depth at 550nm 2021-02-05T18:00:00'}>)

4. Create a geographical subset for Europe ¶

The map above shows Dust Aerosol Optical Depth at 550nm globally. Let us create a geographical subset for Europe, in order to better analyse the Saharan dust event which impacts parts of central and southern Europe.

For geographical subsetting, you can use the function generate_geographical_subset. You can use ?generate_geographical_subset to open the docstring in order to see the function's keyword arguments.

In [13]:

?generate_geographical_subset

Signature:
generate_geographical_subset(
    xarray,
    latmin,
    latmax,
    lonmin,
    lonmax,
    reassign=False,
)
Docstring:
Generates a geographical subset of a xarray.DataArray and if kwarg reassign=True, shifts the longitude grid 
from a 0-360 to a -180 to 180 deg grid.

Parameters:
    xarray(xarray.DataArray): a xarray DataArray with latitude and longitude coordinates
    latmin, latmax, lonmin, lonmax(int): lat/lon boundaries of the geographical subset
    reassign(boolean): default is False
    
Returns:
    Geographical subset of a xarray.DataArray.
File:      /tmp/ipykernel_783/3979307327.py
Type:      function

Define the bounding box information for Europe

In [14]:

latmin = 28.
latmax = 71.
lonmin = -22.
lonmax = 43

Now, let us apply the function generate_geographical_subset to subset the du_aod xarray.DataArray. Let us call the new xarray.DataArray du_aod_subset.

In [15]:

du_aod_subset = generate_geographical_subset(xarray=du_aod, 
                                             latmin=latmin, 
                                             latmax=latmax, 
                                             lonmin=lonmin, 
                                             lonmax=lonmax)
du_aod_subset

Out[15]:

<xarray.DataArray 'duaod550' (time: 31, latitude: 107, longitude: 162)>
array([[[2.9331446e-04, 2.9331446e-04, 2.9331446e-04, ...,
         1.0535121e-03, 1.0535121e-03, 1.0839105e-03],
        [2.6291609e-04, 2.6291609e-04, 2.6291609e-04, ...,
         1.0839105e-03, 1.0839105e-03, 1.1143088e-03],
        [2.6291609e-04, 2.6291609e-04, 2.6291609e-04, ...,
         1.0839105e-03, 1.0839105e-03, 1.0839105e-03],
        ...,
        [4.4536591e-04, 4.4536591e-04, 5.3656101e-04, ...,
         1.2706566e-01, 1.3022816e-01, 1.6598833e-01],
        [5.0616264e-04, 5.0616264e-04, 5.3656101e-04, ...,
         1.2873816e-01, 1.2736976e-01, 1.5762603e-01],
        [5.9741735e-04, 5.9741735e-04, 6.2781572e-04, ...,
         1.2645751e-01, 1.2381196e-01, 1.4801705e-01]],

       [[2.9331446e-04, 2.9331446e-04, 2.9331446e-04, ...,
         9.0146065e-04, 9.3185902e-04, 9.9271536e-04],
        [3.2371283e-04, 2.9331446e-04, 2.9331446e-04, ...,
         1.0535121e-03, 1.1143088e-03, 1.2055635e-03],
        [3.2371283e-04, 2.9331446e-04, 2.9331446e-04, ...,
         1.1751652e-03, 1.2055635e-03, 1.2663603e-03],
...
        [1.8441081e-03, 1.8441081e-03, 1.8137097e-03, ...,
         2.2704244e-02, 3.1218588e-02, 4.4871926e-02],
        [1.6008615e-03, 1.6008615e-03, 1.5096664e-03, ...,
         2.5927544e-02, 3.4016132e-02, 4.6605229e-02],
        [1.3576150e-03, 1.3272166e-03, 1.2663603e-03, ...,
         3.2343686e-02, 3.9520085e-02, 4.9919724e-02]],

       [[7.7986717e-04, 8.1026554e-04, 8.1026554e-04, ...,
         8.4066391e-04, 8.4066391e-04, 8.4066391e-04],
        [7.7986717e-04, 7.7986717e-04, 7.7986717e-04, ...,
         7.7986717e-04, 7.7986717e-04, 7.7986717e-04],
        [7.4940920e-04, 7.7986717e-04, 7.7986717e-04, ...,
         7.4940920e-04, 7.4940920e-04, 7.4940920e-04],
        ...,
        [1.9353628e-03, 1.9657612e-03, 1.9657612e-03, ...,
         2.5866747e-02, 3.3499241e-02, 4.3564379e-02],
        [1.8745661e-03, 1.9049644e-03, 1.8441081e-03, ...,
         2.8846741e-02, 3.7178636e-02, 4.7821522e-02],
        [1.6921163e-03, 1.6616583e-03, 1.6008615e-03, ...,
         3.3316791e-02, 4.0675581e-02, 5.0041378e-02]]], dtype=float32)
Coordinates:
  * longitude  (longitude) float32 -21.6 -21.2 -20.8 -20.4 ... 42.0 42.4 42.8
  * latitude   (latitude) float32 70.8 70.4 70.0 69.6 ... 29.6 29.2 28.8 28.4
  * time       (time) datetime64[ns] 2021-02-04T12:00:00 ... 2021-02-08T06:00:00
Attributes:
    units:      ~
    long_name:  Dust Aerosol Optical Depth at 550nm

xarray.DataArray

'duaod550'

time: 31
latitude: 107
longitude: 162

0.0002933 0.0002933 0.0002933 0.0002933 ... 0.03332 0.04068 0.05004

array([[[2.9331446e-04, 2.9331446e-04, 2.9331446e-04, ...,
         1.0535121e-03, 1.0535121e-03, 1.0839105e-03],
        [2.6291609e-04, 2.6291609e-04, 2.6291609e-04, ...,
         1.0839105e-03, 1.0839105e-03, 1.1143088e-03],
        [2.6291609e-04, 2.6291609e-04, 2.6291609e-04, ...,
         1.0839105e-03, 1.0839105e-03, 1.0839105e-03],
        ...,
        [4.4536591e-04, 4.4536591e-04, 5.3656101e-04, ...,
         1.2706566e-01, 1.3022816e-01, 1.6598833e-01],
        [5.0616264e-04, 5.0616264e-04, 5.3656101e-04, ...,
         1.2873816e-01, 1.2736976e-01, 1.5762603e-01],
        [5.9741735e-04, 5.9741735e-04, 6.2781572e-04, ...,
         1.2645751e-01, 1.2381196e-01, 1.4801705e-01]],

       [[2.9331446e-04, 2.9331446e-04, 2.9331446e-04, ...,
         9.0146065e-04, 9.3185902e-04, 9.9271536e-04],
        [3.2371283e-04, 2.9331446e-04, 2.9331446e-04, ...,
         1.0535121e-03, 1.1143088e-03, 1.2055635e-03],
        [3.2371283e-04, 2.9331446e-04, 2.9331446e-04, ...,
         1.1751652e-03, 1.2055635e-03, 1.2663603e-03],
...
        [1.8441081e-03, 1.8441081e-03, 1.8137097e-03, ...,
         2.2704244e-02, 3.1218588e-02, 4.4871926e-02],
        [1.6008615e-03, 1.6008615e-03, 1.5096664e-03, ...,
         2.5927544e-02, 3.4016132e-02, 4.6605229e-02],
        [1.3576150e-03, 1.3272166e-03, 1.2663603e-03, ...,
         3.2343686e-02, 3.9520085e-02, 4.9919724e-02]],

       [[7.7986717e-04, 8.1026554e-04, 8.1026554e-04, ...,
         8.4066391e-04, 8.4066391e-04, 8.4066391e-04],
        [7.7986717e-04, 7.7986717e-04, 7.7986717e-04, ...,
         7.7986717e-04, 7.7986717e-04, 7.7986717e-04],
        [7.4940920e-04, 7.7986717e-04, 7.7986717e-04, ...,
         7.4940920e-04, 7.4940920e-04, 7.4940920e-04],
        ...,
        [1.9353628e-03, 1.9657612e-03, 1.9657612e-03, ...,
         2.5866747e-02, 3.3499241e-02, 4.3564379e-02],
        [1.8745661e-03, 1.9049644e-03, 1.8441081e-03, ...,
         2.8846741e-02, 3.7178636e-02, 4.7821522e-02],
        [1.6921163e-03, 1.6616583e-03, 1.6008615e-03, ...,
         3.3316791e-02, 4.0675581e-02, 5.0041378e-02]]], dtype=float32)

Coordinates: (3)

longitude

(longitude)

float32

-21.6 -21.2 -20.8 ... 42.4 42.8

array([-21.599976, -21.200012, -20.799988, -20.400024, -20.      , -19.599976,
       -19.200012, -18.799988, -18.400024, -18.      , -17.599976, -17.200012,
       -16.799988, -16.400024, -16.      , -15.599976, -15.200012, -14.799988,
       -14.400024, -14.      , -13.599976, -13.200012, -12.799988, -12.400024,
       -12.      , -11.599976, -11.200012, -10.799988, -10.400024, -10.      ,
        -9.599976,  -9.200012,  -8.799988,  -8.400024,  -8.      ,  -7.599976,
        -7.200012,  -6.799988,  -6.400024,  -6.      ,  -5.599976,  -5.200012,
        -4.799988,  -4.400024,  -4.      ,  -3.599976,  -3.200012,  -2.799988,
        -2.400024,  -2.      ,  -1.599976,  -1.200012,  -0.799988,  -0.400024,
         0.      ,   0.399994,   0.800003,   1.199997,   1.600006,   2.      ,
         2.399994,   2.800003,   3.199997,   3.600006,   4.      ,   4.399994,
         4.800003,   5.199997,   5.600006,   6.      ,   6.399994,   6.800003,
         7.199997,   7.600006,   8.      ,   8.399994,   8.800003,   9.199997,
         9.600006,  10.      ,  10.399994,  10.800003,  11.199997,  11.600006,
        12.      ,  12.399994,  12.800003,  13.199997,  13.600006,  14.      ,
        14.399994,  14.800003,  15.199997,  15.600006,  16.      ,  16.399994,
        16.800003,  17.199997,  17.600006,  18.      ,  18.399994,  18.800003,
        19.199997,  19.600006,  20.      ,  20.399994,  20.800003,  21.199997,
        21.600006,  22.      ,  22.399994,  22.800003,  23.199997,  23.600006,
        24.      ,  24.399994,  24.800003,  25.199997,  25.600006,  26.      ,
        26.399994,  26.800003,  27.199997,  27.600006,  28.      ,  28.399994,
        28.800003,  29.199997,  29.600006,  30.      ,  30.399994,  30.800003,
        31.199997,  31.600006,  32.      ,  32.399994,  32.800003,  33.199997,
        33.600006,  34.      ,  34.399994,  34.800003,  35.199997,  35.600006,
        36.      ,  36.399994,  36.800003,  37.199997,  37.600006,  38.      ,
        38.399994,  38.800003,  39.199997,  39.600006,  40.      ,  40.399994,
        40.800003,  41.199997,  41.600006,  42.      ,  42.399994,  42.800003],
      dtype=float32)

latitude

(latitude)

float32

70.8 70.4 70.0 ... 29.2 28.8 28.4

units :: degrees_north
long_name :: latitude

array([70.8, 70.4, 70. , 69.6, 69.2, 68.8, 68.4, 68. , 67.6, 67.2, 66.8, 66.4,
       66. , 65.6, 65.2, 64.8, 64.4, 64. , 63.6, 63.2, 62.8, 62.4, 62. , 61.6,
       61.2, 60.8, 60.4, 60. , 59.6, 59.2, 58.8, 58.4, 58. , 57.6, 57.2, 56.8,
       56.4, 56. , 55.6, 55.2, 54.8, 54.4, 54. , 53.6, 53.2, 52.8, 52.4, 52. ,
       51.6, 51.2, 50.8, 50.4, 50. , 49.6, 49.2, 48.8, 48.4, 48. , 47.6, 47.2,
       46.8, 46.4, 46. , 45.6, 45.2, 44.8, 44.4, 44. , 43.6, 43.2, 42.8, 42.4,
       42. , 41.6, 41.2, 40.8, 40.4, 40. , 39.6, 39.2, 38.8, 38.4, 38. , 37.6,
       37.2, 36.8, 36.4, 36. , 35.6, 35.2, 34.8, 34.4, 34. , 33.6, 33.2, 32.8,
       32.4, 32. , 31.6, 31.2, 30.8, 30.4, 30. , 29.6, 29.2, 28.8, 28.4],
      dtype=float32)

time

(time)

datetime64[ns]

2021-02-04T12:00:00 ... 2021-02-...

long_name :: time

array(['2021-02-04T12:00:00.000000000', '2021-02-04T15:00:00.000000000',
       '2021-02-04T18:00:00.000000000', '2021-02-04T21:00:00.000000000',
       '2021-02-05T00:00:00.000000000', '2021-02-05T03:00:00.000000000',
       '2021-02-05T06:00:00.000000000', '2021-02-05T09:00:00.000000000',
       '2021-02-05T12:00:00.000000000', '2021-02-05T15:00:00.000000000',
       '2021-02-05T18:00:00.000000000', '2021-02-05T21:00:00.000000000',
       '2021-02-06T00:00:00.000000000', '2021-02-06T03:00:00.000000000',
       '2021-02-06T06:00:00.000000000', '2021-02-06T09:00:00.000000000',
       '2021-02-06T12:00:00.000000000', '2021-02-06T15:00:00.000000000',
       '2021-02-06T18:00:00.000000000', '2021-02-06T21:00:00.000000000',
       '2021-02-07T00:00:00.000000000', '2021-02-07T03:00:00.000000000',
       '2021-02-07T06:00:00.000000000', '2021-02-07T09:00:00.000000000',
       '2021-02-07T12:00:00.000000000', '2021-02-07T15:00:00.000000000',
       '2021-02-07T18:00:00.000000000', '2021-02-07T21:00:00.000000000',
       '2021-02-08T00:00:00.000000000', '2021-02-08T03:00:00.000000000',
       '2021-02-08T06:00:00.000000000'], dtype='datetime64[ns]')

Attributes: (2)
units :
~
long_name :
Dust Aerosol Optical Depth at 550nm

Let us visualize the subsetted xarray.DataArray again. This time, you set the set_global kwarg to False and you specify the longitude and latitude bounds specified above.

Additionally, in order to have the time information as part of the title, we add the string of the datetime information to the long_name variable: long_name + ' ' + str(du_aod_subset[##,:,:].time.data)[0:19].

In [16]:

time_index = 14
visualize_pcolormesh(data_array=du_aod_subset[time_index,:,:],
                     longitude=du_aod_subset.longitude,
                     latitude=du_aod_subset.latitude,
                     projection=ccrs.PlateCarree(),
                     color_scale='hot_r',
                     unit=units,
                     long_name=long_name + ' ' + str(du_aod_subset[time_index,:,:].time.data)[0:19],
                     vmin=0, 
                     vmax=1,
                     lonmin=lonmin,
                     lonmax=lonmax,
                     latmin=latmin,
                     latmax=latmax,
                     set_global=False)

Out[16]:

(<Figure size 1440x720 with 2 Axes>,
 <GeoAxesSubplot:title={'center':'Dust Aerosol Optical Depth at 550nm 2021-02-06T06:00:00'}>)

5. Animate changes of Dust Aerosol Optical Depth at 550nm over time ¶

In the last step, you can animate the Dust Aerosol Optical Depth at 550nm in order to see how the trace gas develops over a period of 5 days, from 4th to 8th February 2021.

You can do animations with matplotlib's function animation. Jupyter's function HTML can then be used to display HTML and video content.

The animation function consists of 4 parts:

Setting the initial state:

Here, you define the general plot your animation shall use to initialise the animation. You can also define the number of frames (time steps) your animation shall have.

Functions to animate:

An animation consists of three functions: draw(), init() and animate(). draw() is the function where individual frames are passed on and the figure is returned as image. In this example, the function redraws the plot for each time step. init() returns the figure you defined for the initial state. animate() returns the draw() function and animates the function over the given number of frames (time steps).

Create a animate.FuncAnimation object:

The functions defined before are now combined to build an animate.FuncAnimation object.

Play the animation as video:

As a final step, you can integrate the animation into the notebook with the HTML class. You take the generate animation object and convert it to a HTML5 video with the to_html5_video function

In [17]:

# Setting the initial state:
# 1. Define figure for initial plot
fig, ax = visualize_pcolormesh(data_array=du_aod_subset[0,:,:],
                               longitude=du_aod_subset.longitude, 
                               latitude=du_aod_subset.latitude,
                               projection=ccrs.PlateCarree(), 
                               color_scale='hot_r', 
                               unit='-',
                               long_name=long_name + ' '+ str(du_aod_subset.time[0].data)[0:19],
                               vmin=0,
                               vmax=1, 
                               lonmin=lonmin, 
                               lonmax=lonmax, 
                               latmin=latmin, 
                               latmax=latmax,
                               set_global=False)

frames = 30

def draw(i):
    img = plt.pcolormesh(du_aod_subset.longitude, 
                         du_aod_subset.latitude, 
                         du_aod_subset[i,:,:], 
                         cmap='hot_r', 
                         transform=ccrs.PlateCarree(),
                         vmin=0,
                         vmax=1,
                         shading='auto')
    
    ax.set_title(long_name + ' '+ str(du_aod_subset.time[i].data)[0:19], fontsize=20, pad=20.0)
    return img

def init():
    return fig

def animate(i):
    return draw(i)

ani = animation.FuncAnimation(fig, animate, frames, interval=800, blit=False,
                              init_func=init, repeat=True)

HTML(ani.to_html5_video())
plt.close(fig)

Play the animation video as HTML5 video¶

In [18]:

HTML(ani.to_html5_video())

Out[18]:

References¶

Copernicus Service information 2021
Generated using Copernicus Atmosphere Monitoring Service Information 2021

Logo EU Copernicus

This project is licensed under GNU General Public License v3.0 only and is developed under a Copernicus contract.

EUMETSAT Training | Contact the training team