Antarctic Ozone Hole Monitoring Tutorial

This tutorial provides guided examples of ozone monitoring using data from the Copernicus Atmosphere Monitoring Service (CAMS). It is divided into two parts:

1. View animation of Antarctic ozone hole

2. Calculate the size of the Antarctic ozone hole

It uses CAMS global reanalysis (EAC4) data freely available from the Atmosphere Data Store (ADS)

Run the tutorial via free cloud platforms:

Install and import packages

!pip install cdsapi

# CDS API
import cdsapi

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

# Libraries for plotting and animating data
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import colors
import matplotlib.animation as animation
import cartopy.crs as ccrs
from IPython.display import HTML

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

from collections import OrderedDict

Download and read and pre-process ozone data

Download 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-beta.atmosphere.copernicus.eu/ and obtain an API key from https://ads-beta.atmosphere.copernicus.eu/how-to-api.)

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

# 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 this tutorial, we will use CAMS Global Reanalysis (EAC4) data. The code below shows the subset characteristics that we will extract from this dataset 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.

dataset = "cams-global-reanalysis-eac4"
request = {
    'variable': ['total_column_ozone'],
    'date': ['2020-07-01/2021-01-31'],
    'time': ['00:00'],
    'data_format': 'netcdf',
    'area': [0, -180, -90, 180]
}

client = cdsapi.Client(url=URL, key=KEY)
client.retrieve(dataset, request).download(
    f'{DATADIR}/TCO3_202007-202101_SHem.nc')

2024-09-12 14:55:56,391 INFO Request ID is 783294a0-7119-4379-b0dc-f8d994b1b3c1
2024-09-12 14:55:56,536 INFO status has been updated to accepted
2024-09-12 14:55:58,187 INFO status has been updated to running
2024-09-12 14:56:28,825 INFO Creating download object as as_source with files:
['data_sfc.nc']
2024-09-12 14:56:28,827 INFO status has been updated to successful

'./TCO3_202007-202101_SHem.nc'

Import and inspect data

Here we read the data we downloaded into an xarray Dataset, and view its contents:

fn = f'{DATADIR}/TCO3_202007-202101_SHem.nc'
ds = xr.open_dataset(fn)
ds

<xarray.Dataset> Size: 50MB
Dimensions:     (valid_time: 215, latitude: 121, longitude: 480)
Coordinates:
  * valid_time  (valid_time) datetime64[ns] 2kB 2020-07-01 ... 2021-01-31
  * latitude    (latitude) float64 968B 0.0 -0.75 -1.5 ... -88.5 -89.25 -90.0
  * longitude   (longitude) float64 4kB -180.0 -179.2 -178.5 ... 178.5 179.2
Data variables:
    gtco3       (valid_time, latitude, longitude) float32 50MB ...
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:                 2024-09-12T14:56 GRIB to CDM+CF via cfgrib-0.9.1...

xarray.Dataset

• Dimensions:
  • valid_time: 215
  • latitude: 121
  • longitude: 480
• Coordinates: (3)
  • valid_time
    (valid_time)
    datetime64[ns]
    2020-07-01 ... 2021-01-31

    long_name :

    time

    standard_name :

    time

    array(['2020-07-01T00:00:00.000000000', '2020-07-02T00:00:00.000000000',
           '2020-07-03T00:00:00.000000000', ..., '2021-01-29T00:00:00.000000000',
           '2021-01-30T00:00:00.000000000', '2021-01-31T00:00:00.000000000'],
          dtype='datetime64[ns]')

  • latitude
    (latitude)
    float64
    0.0 -0.75 -1.5 ... -89.25 -90.0

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    latitude

    stored_direction :

    decreasing

    array([  0.  ,  -0.75,  -1.5 ,  -2.25,  -3.  ,  -3.75,  -4.5 ,  -5.25,  -6.  ,
            -6.75,  -7.5 ,  -8.25,  -9.  ,  -9.75, -10.5 , -11.25, -12.  , -12.75,
           -13.5 , -14.25, -15.  , -15.75, -16.5 , -17.25, -18.  , -18.75, -19.5 ,
           -20.25, -21.  , -21.75, -22.5 , -23.25, -24.  , -24.75, -25.5 , -26.25,
           -27.  , -27.75, -28.5 , -29.25, -30.  , -30.75, -31.5 , -32.25, -33.  ,
           -33.75, -34.5 , -35.25, -36.  , -36.75, -37.5 , -38.25, -39.  , -39.75,
           -40.5 , -41.25, -42.  , -42.75, -43.5 , -44.25, -45.  , -45.75, -46.5 ,
           -47.25, -48.  , -48.75, -49.5 , -50.25, -51.  , -51.75, -52.5 , -53.25,
           -54.  , -54.75, -55.5 , -56.25, -57.  , -57.75, -58.5 , -59.25, -60.  ,
           -60.75, -61.5 , -62.25, -63.  , -63.75, -64.5 , -65.25, -66.  , -66.75,
           -67.5 , -68.25, -69.  , -69.75, -70.5 , -71.25, -72.  , -72.75, -73.5 ,
           -74.25, -75.  , -75.75, -76.5 , -77.25, -78.  , -78.75, -79.5 , -80.25,
           -81.  , -81.75, -82.5 , -83.25, -84.  , -84.75, -85.5 , -86.25, -87.  ,
           -87.75, -88.5 , -89.25, -90.  ])

  • longitude
    (longitude)
    float64
    -180.0 -179.2 ... 178.5 179.2

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    longitude

    array([-180.  , -179.25, -178.5 , ...,  177.75,  178.5 ,  179.25])

• Data variables: (1)
  • gtco3
    (valid_time, latitude, longitude)
    float32
    ...

    GRIB_paramId :

    210206

    GRIB_dataType :

    an

    GRIB_numberOfPoints :

    58080

    GRIB_typeOfLevel :

    surface

    GRIB_stepUnits :

    1

    GRIB_stepType :

    instant

    GRIB_gridType :

    regular_ll

    GRIB_uvRelativeToGrid :

    0

    GRIB_NV :

    0

    GRIB_Nx :

    480

    GRIB_Ny :

    121

    GRIB_cfName :

    atmosphere_mass_content_of_ozone

    GRIB_cfVarName :

    gtco3

    GRIB_gridDefinitionDescription :

    Latitude/Longitude Grid

    GRIB_iDirectionIncrementInDegrees :

    0.75

    GRIB_iScansNegatively :

    0

    GRIB_jDirectionIncrementInDegrees :

    0.75

    GRIB_jPointsAreConsecutive :

    0

    GRIB_jScansPositively :

    0

    GRIB_latitudeOfFirstGridPointInDegrees :

    0.0

    GRIB_latitudeOfLastGridPointInDegrees :

    -90.0

    GRIB_longitudeOfFirstGridPointInDegrees :

    -180.0

    GRIB_longitudeOfLastGridPointInDegrees :

    179.25

    GRIB_missingValue :

    3.4028234663852886e+38

    GRIB_name :

    GEMS Total column ozone

    GRIB_shortName :

    gtco3

    GRIB_totalNumber :

    0

    GRIB_units :

    kg m**-2

    long_name :

    GEMS Total column ozone

    units :

    kg m**-2

    standard_name :

    atmosphere_mass_content_of_ozone

    GRIB_number :

    0

    GRIB_surface :

    0.0

    [12487200 values with dtype=float32]

• Indexes: (3)
  • valid_time
    PandasIndex

    PandasIndex(DatetimeIndex(['2020-07-01', '2020-07-02', '2020-07-03', '2020-07-04',
                   '2020-07-05', '2020-07-06', '2020-07-07', '2020-07-08',
                   '2020-07-09', '2020-07-10',
                   ...
                   '2021-01-22', '2021-01-23', '2021-01-24', '2021-01-25',
                   '2021-01-26', '2021-01-27', '2021-01-28', '2021-01-29',
                   '2021-01-30', '2021-01-31'],
                  dtype='datetime64[ns]', name='valid_time', length=215, freq=None))

  • latitude
    PandasIndex

    PandasIndex(Index([   0.0,  -0.75,   -1.5,  -2.25,   -3.0,  -3.75,   -4.5,  -5.25,   -6.0,
            -6.75,
           ...
           -83.25,  -84.0, -84.75,  -85.5, -86.25,  -87.0, -87.75,  -88.5, -89.25,
            -90.0],
          dtype='float64', name='latitude', length=121))

  • longitude
    PandasIndex

    PandasIndex(Index([ -180.0, -179.25,  -178.5, -177.75,  -177.0, -176.25,  -175.5, -174.75,
            -174.0, -173.25,
           ...
             172.5,  173.25,   174.0,  174.75,   175.5,  176.25,   177.0,  177.75,
             178.5,  179.25],
          dtype='float64', name='longitude', length=480))

• Attributes: (6)

  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 :

  2024-09-12T14:56 GRIB to CDM+CF via cfgrib-0.9.14.0/ecCodes-2.36.0 with {"source": "data.grib", "filter_by_keys": {"typeOfLevel": "surface"}, "encode_cf": ["parameter", "time", "geography", "vertical"]}

We can see that the data has three coordinate dimensions (longitude, latitude and time) and one variable, total column ozone. By inspecting the coordinates we can see the data is of the southern hemisphere from 1st July 2020 to 31 January 2021 at 00:00 UTC each day. This includes the entire period in which the Antarctic ozone hole appears.

Note

This is only a subset of the CAMS Global Reanalysis data available on the ADS, which includes global data at 3 hourly resolution (in addition to monthly averages), from 2003 to the present, and at 60 model levels (vertical layers) in the atmosphere.

To facilitate further processing, we convert the Xarray Dataset into an Xarray Data Array containing the single variable of total column ozone.

tco3 = ds['gtco3']

Unit conversion

We can see from the attributes of our data are represented as mass concentration of ozone, in units of kg m**-2. We would like to convert this into Dobson Units, which is a standard for ozone measurements. The Dobson unit is defined as the thickness (in units of 10 μm) of a layer of pure gas (in this case O3) which would be formed by the total column amount at standard conditions for temperature and pressure.

convertToDU = 1 / 2.1415e-5

tco3 = tco3 * convertToDU

View Antarctic ozone hole

Let us now visualise our data in form of maps and animations.

We will first define the colour scale we would like to represent our data.

Define colour scale

Extract range of values

tco3max = tco3.max()
tco3min = tco3.min()
tco3range = (tco3max - tco3min)/20.

Define colourmap

cmap = plt.cm.jet
norm = colors.BoundaryNorm(np.arange(tco3min, tco3max, tco3range), cmap.N)

Plot map

We will now plot our data. Let us begin with a map for the first time step, 1 July 2020.

fig = plt.figure(figsize=(5, 5)) 
ax = plt.subplot(1,1,1, projection=ccrs.Orthographic(central_latitude=-90)) 
ax.coastlines(color='black') # Add coastlines
ax.gridlines(draw_labels=True, linewidth=1, color='gray', alpha=0.5, linestyle='--') 
ax.set_title(f'Total Column Ozone, {str(tco3.valid_time[0].values)[:-19]}', fontsize=12) 
im = plt.pcolormesh(tco3.longitude.values, tco3.latitude.values, tco3[0,:,:],
                cmap=cmap, norm=norm, transform=ccrs.PlateCarree())
cbar = plt.colorbar(im,fraction=0.046, pad=0.04) 
cbar.set_label('Total Column Ozone (DU)') 

/opt/conda/lib/python3.10/site-packages/cartopy/io/__init__.py:241: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/110m_physical/ne_110m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

Create animation

Let us now view the entire time series by means of an animation that shows all time steps. We see from the dataset description above, that the time dimension in our data includes 215 entries. These correspond to 215 time steps, i.e. 215 days, between 1 July 2020 and 31 January 2021. The number of frames of our animation is therefore 215.

frames = 215

Now we will create an animation. The last cell, which creates the animation in Javascript, may take a few minutes to run.

def animate(i):
    array = tco3[i,:,:].values
    im.set_array(array.flatten())
    ax.set_title(f'Total Column Ozone, {str(tco3.valid_time[i].values)[:-19]}', fontsize=12)

ani = animation.FuncAnimation(fig, animate, frames, interval=100)
ani.save(f'{DATADIR}/TCO3_202007-202101_SHem.gif')
HTML(ani.to_jshtml())

Once Loop Reflect

ACTIVITY:
In what date range does the ozone hole appear greatest? In the next section we will quantify the size of the ozone hole at each time step.

Calculate size of ozone hole

In this section we will calculate the actual size of the ozone hole, and plot its evolution over our time series. We will then apply the exact same method on monthly averaged data over the past decades to compare ozone hole extents at each yearly occurence since 2003.

Define constants

In order to calculate the size of the ozone hole, we need to define some constants. These include the following:

Ozone hole limit, i.e. minimum ozone, in Dobson units, beneath which we consider there to be a hole in the ozone layer.

OZONE_HOLE_LIMIT = 220.

Radius of the Earth, in order to calculate the size of each geographic grid cell of our data.

Rearth = 6371009.

Calculate area of each geographic grid cell of our data

Here we calculate the size of each grid cell of our data, which varies as a function of latitude. The formula we apply is the following:

area = |(sin(top latitude) - sin(bottom latitude))| * |Δlongitude| * Earth radius squared

First we define a function to apply this formula, taking as arguments the min and max of the latitudes and longitudes represented by each grid cell:

def geo_area_array(lats1,lons1,lats2,lons2):
    area= np.abs(np.sin(np.deg2rad(lats1)) - np.sin(np.deg2rad(lats2))) * np.abs(np.deg2rad(lons1-lons2))
    area = area * Rearth * Rearth 
    return area

Then we apply this function to our data, defining first the longitude and latitude delta at each grid cell.

DeltaLat = 0.75 # Latitudinal spatial resolution, in degrees
DeltaLon = 0.75 # Longitudinal spatial resolution, in degrees

tco3['areas'] = geo_area_array(tco3['latitude']-DeltaLat/2, tco3['longitude']-DeltaLat/2, 
                             tco3['latitude']+DeltaLat/2, tco3['longitude']+DeltaLat/2)

Mask data belonging to Antarctic ozone hole

We now calculate the grid cells that meet the criteria of belonging to the Antarctic ozone hole: i.e. they have less than the minimum threshold of ozone, and they are below -60 degrees latitude.

Here we define a mask with the threshold conditions for ozone and latitude:

mask = tco3.where((tco3 < OZONE_HOLE_LIMIT) & (tco3["latitude"] < -60.))

All data points that do not meet the conditions are set to NaN (Not a Number) in the resulting array. Dividing this array by itself gives us a mask with each valid data point set to 1. Multiplying this by the corresponding areas leaves us with an array of grid cell areas that meet the conditions of belonging to the Antarctic ozone hole.

area = (mask / mask) * mask['areas']

Calculate ozone hole area

If we now sum these grid cell areas, we have the total ozone hole extent. This is applied across each time step. We multiply by the number 1e-12 to convert from meters squared to million kilometers squared.

ozone_hole = area.sum(dim=["latitude", "longitude"], skipna=True) * 1e-12

Let us update the Data Array attributes:

ozone_hole.attrs['long_name'] = 'Ozone hole area'
ozone_hole.attrs['units'] = 'million km^2'

Plot ozone hole area at each time step

Finally we can plot the evolution of the ozone hole, in million km squared, throughout the time series.

ozone_hole.plot()

[<matplotlib.lines.Line2D at 0x7ef2d87e9570>]

We can also calculate the maximum extent of the ozone hole:

index = ozone_hole.argmax() # Here we find the index of the maximum value
print('Max extent of O3 hole:', ozone_hole[index].values, 'million km2')

Max extent of O3 hole: 23.757288326941282 million km2

… and find out when this was reached:

print('This was reached on:', ozone_hole[index].valid_time.values)

This was reached on: 2020-09-28T00:00:00.000000000

ACTIVITY:
Now repeat the steps above with the data below (see API request in cell below). See what you find out about the evolution of the ozone hole since 2003!

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

client.retrieve(dataset, request).download(
    f'{DATADIR}/TCO3_monthly_2003-2023.nc')

2024-09-12 15:13:26,368 INFO Request ID is c15a068f-9137-464c-8555-f092b0731fd0
2024-09-12 15:13:26,516 INFO status has been updated to accepted
2024-09-12 15:13:28,169 INFO status has been updated to running
2024-09-12 15:13:47,051 INFO Creating download object as as_source with files:
['data_allhours_sfc.nc']
2024-09-12 15:13:58,592 INFO status has been updated to successful

'./TCO3_monthly_2003-2023.nc'