The CAMS service provides the Air Quality Forecasts for Europe. The data product is based on an ensemble of the outputs of nine Earth System models that provides the forecasts for many polluttants such as $NO_x$, Ozone, particulate matter (PM2.5 and PM10) and $SO_2$ with a spatial resolution of 10 km. In this test we want to retrieve the forecasts for $NO_2$ over the Italian peninsula at ground level every hours for the next 96 hours (lead times) and visualize the forecasts as a sequence of frames in an animation. The forecasts for the air quality are the result of computer simulations, performed by the ECMWF, using a technique called data assimilation. The technique uses new observations, in our case from the TROPOMI instrument onboard the Copernicus Sentinel-5p, to update the initial conditions of a numerical model of the atmosphere that will be used for the next simulation that will represent the most up to date forecast. More information about the origin of NO2 and its relevance on the air quality is available at the EUMETSAT website and from the Julien Chimot's blog.
The TROPOspheric Monitoring Instrument (TROPOMI) is the only payload of Sentinel-5p and is dedicated to the monitoring of trace gases and particulate in the atmosphere and in particular in its lower layer, the troposphere. It is a passive remote sensing instrument and depends on the sun light reflected from the Earth surface through the atmosphere. The trace gases monitored are: ozone $(O_3)$, nitrogen dioxide $(NO_2)$, sulfur dioxide $(SO_2)$, carbon monoxide $(CO)$, methane $(CH_4)$, formaldehyde $(HCHO)$. The instrument is composed of 4 spectrometers that collect the light coming from nadir in 4 spectral bands: UV (267-332 nm), UV-VIS (305-499 nm), NIR (661-786 nm), and SWIR (2300-2389 nm). The light arriving at the instrument is filtered, sent to the detector of each spectrometer, transformed into a voltage signal, amplified and finally digitized in order to be transferred to the ground station. The detector of the UV-VIS-NIR bands is a CCD. The detector of the SWIR band is CMOS device (see [3]). The total and tropospheric column of the trace gases is extracted from the detector signal using the Differential optical absorption spectroscopy (DOAS). The technique is based on the Beer-Lambert law for which the intensity of the light at a particular wavelenght received at the detector, depends on the amount of light that is absorbed in the path between the source and the detector, and on the density of the absorbing molecules. The Beer-Lambert law has the form
$$I(\lambda) = I_0(\lambda)exp(-\sum_{j=1}^n \sigma_j(\lambda)\rho_j)$$where $\lambda$ is the wavelenght of the detected light, $\sigma_j(\lambda)$ is the absorption cross section of the species j, and $\rho_j$ is its column density. Assuming $\sigma_j(\lambda)$ are known from laboratory measurements, taking n measures of $I(\lambda)$ we have a system of n equations in n unknown variables $\rho_j$ that can be solved (see [6]). After this very short and simplified description of the instrument and the inversion algorithm, we can proceed using the data.
import xarray as xr
import pandas as pd
import numpy as np
import cdsapi
import warnings
warnings.filterwarnings('ignore')
from platform import python_version
print("python version: %s"%python_version())
print("pandas version: %s"%pd.__version__)
print("xarray version: %s"%xr.__version__)
python version: 3.11.5 pandas version: 2.1.1 xarray version: 2023.12.0
The Copernicus Atmosphere Monitoring Service (CAMS) provides analysis and forecasts related to air quality, atmospheric composition, greenhouse gases, solar irradiance. The datasets released by the CAMS service is the result of assimilation processes in which observations from satellites and ground stations are used to update a numerical model. The CAMS is operated by the European Centre for Medium-Range Weather Forecasts (ECMWF) on behalf of the European Commission.
The CAMS provides its datasets as open data, available to all for free, through a web page and a web service. For the air quality forecasts a user can select, among other options
In this notebook we use the web service through its API. In order to use the API you need an API key that will be given to you after you have registered into the CADS. The steps to use the CAMS web service are
You can also pass your API key and the web service url to the client API instead of copying them into your .condarc file. For more information follow the how-to instructions.
We set some parameters before submitting our request, in particular the area of interest, the variable (NO2), the lead time hours.
path = 'data/air_quality_forecasts/20240220'
file_name = 'air_quality_forecasts_europe_20240220.nc'
date = '2024-02-20'
area_north = 47.12
area_south = 36.40
area_west = 6.57
area_east = 18.52
URL = "https://ads.atmosphere.copernicus.eu/api/v2"
KEY = "xxxxxxxxxxxxxxxxxxxxxxxxxxxx"
c = cdsapi.Client(url=URL, key=KEY) # uncomment with your API key if you do not use the .condarc file
#c = cdsapi.Client() # comment this line if you are sending your API key in the request
c.retrieve(
'cams-europe-air-quality-forecasts',
{
'model': 'ensemble',
'date': date + '/' + date,
'format': 'netcdf',
'variable': 'nitrogen_dioxide',
'level': '0',
'type': 'forecast',
'time': '00:00',
'leadtime_hour': [
'0', '1', '10',
'11', '12', '13',
'14', '15', '16',
'17', '18', '19',
'2', '20', '21',
'22', '23', '24',
'25', '26', '27',
'28', '29', '3',
'30', '31', '32',
'33', '34', '35',
'36', '37', '38',
'39', '4', '40',
'41', '42', '43',
'44', '45', '46',
'47', '48', '49',
'5', '50', '51',
'52', '53', '54',
'55', '56', '57',
'58', '59', '6',
'60', '61', '62',
'63', '64', '65',
'66', '67', '68',
'69', '7', '70',
'71', '72', '73',
'74', '75', '76',
'77', '78', '79',
'8', '80', '81',
'82', '83', '84',
'85', '86', '87',
'88', '89', '9',
'90', '91', '92',
'93', '94', '95',
'96',
],
'area': [
area_north, area_west, area_south,
area_east,
],
},
path + '/' + file_name)
2024-02-20 10:20:18,365 INFO Welcome to the CDS 2024-02-20 10:20:18,367 INFO Sending request to https://ads.atmosphere.copernicus.eu/api/v2/resources/cams-europe-air-quality-forecasts 2024-02-20 10:20:18,412 INFO Request is queued 2024-02-20 10:20:19,441 INFO Request is running 2024-02-20 10:20:39,349 INFO Request is completed 2024-02-20 10:20:39,350 INFO Downloading https://download-0004-ads-clone.copernicus-climate.eu/cache-compute-0004/cache/data8/adaptor.cams_regional_fc.retrieve-1708420834.390126-8128-2-63dd217d-48c9-4662-915e-d43a1ed57433.nc to data/air_quality_forecasts/20240220/air_quality_forecasts_europe_20240220.nc (4.7M) 2024-02-20 10:20:41,544 INFO Download rate 2.1M/s
Result(content_length=4942928,content_type=application/x-netcdf,location=https://download-0004-ads-clone.copernicus-climate.eu/cache-compute-0004/cache/data8/adaptor.cams_regional_fc.retrieve-1708420834.390126-8128-2-63dd217d-48c9-4662-915e-d43a1ed57433.nc)
If we have requested the forecasts starting from only one day, let's say the day in which we send the request, the webservice will send a netcdf file with the forecasts for all the lead times we have specified. If we want the forecasts starting from more days the webservice will send a zip file that contains one netCDF file for each start day. In our test we will request the forecasts starting from one single day, the same day in which we send the request.
#import os
#import zipfile
#with zipfile.ZipFile('data/air_quality_forecasts_europe_20201203.zip', 'r') as zip_ref:
# zip_ref.extractall(path)
#nc_files = os.listdir(path)
#nc_files
Let's open the file to get the forecasts for $NO_2$
def getforecast(path, file_name):
# The forecasts for NO2 are in the 'no2_conc' variable.
# The function puts the data in memory, closes the file
# and returns the NO2 forecast
ds = xr.open_dataset(path + '/' + file_name)
no2_forecast = ds['no2_conc']
ds.close()
return no2_forecast
no2_forecasts = getforecast(path, file_name)
no2_forecasts
<xarray.DataArray 'no2_conc' (time: 97, level: 1, latitude: 107, longitude: 119)> [1235101 values with dtype=float32] Coordinates: * longitude (longitude) float32 6.65 6.75 6.85 6.95 ... 18.25 18.35 18.45 * latitude (latitude) float32 47.05 46.95 46.85 46.75 ... 36.65 36.55 36.45 * level (level) float32 0.0 * time (time) timedelta64[ns] 00:00:00 01:00:00 ... 4 days 00:00:00 Attributes: species: Nitrogen Dioxide units: µg/m3 value: hourly values standard_name: mass_concentration_of_nitrogen_dioxide_in_air
As we can see from the description, the time interval of each forecast is provided as a delta in nanoseconds (ns), from the first lead time to the next. Since we have requested a forecast every hour the delta is 3600 * 10^9, or 3600000000000, nanoseconds. We will create a date time index to be used later on in the plots
delta_time = 3600000000000 # delta between one lead time and the next one
num_forecasts = len(no2_forecasts)
start_day = pd.to_datetime(date)
date_index = start_day + pd.to_timedelta(np.arange(num_forecasts), 'H')
date_index.size
start_day.strftime('%Y-%m-%d %H:%M:%S')
'2024-02-20 00:00:00'
We plot the forecast for one lead time hour to check that everything works fine and then we'll create an animation using all the forecasts.
import matplotlib
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy
print("matplotlib version: %s"%matplotlib.__version__)
print("cartopy version: %s"%cartopy.__version__)
matplotlib version: 3.8.0 cartopy version: 0.22.0
We create the figure that we will use to plot the data
def create_figure():
fig = plt.figure(figsize=(20,10))
ax = plt.axes(projection=ccrs.PlateCarree())
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.LAND, edgecolor='black')
ax.coastlines()
return fig, ax
Let's plot the first forecast.
_, ax = create_figure()
lead_time = 1
first_forecast = no2_forecasts.sel(time = str(delta_time * lead_time), level = 0.0)
first_forecast.plot.pcolormesh(ax=ax, x='longitude', y='latitude', add_colorbar=True, cmap='viridis')
plt.title('Air Quality Forecast - {0:s}'.format(date_index[lead_time - 1].strftime('%Y-%m-%d %H:%M:%S')))
#plt.savefig('images/forecast_example.png')
Text(0.5, 1.0, 'Air Quality Forecast - 2024-02-20 00:00:00')
We can notice the NO2 traces in the Mediterranean sea, in particular in the Strait of Sicily and over the major Italian ports. The traces are due to combustion of fossil fuels of the ships. Now we create the animation using all the forecasts at ground level. The animation will contain a frame for each forecast (or lead time)
import matplotlib.animation as animation
from IPython.display import HTML, display
forecasts = no2_forecasts.sel(level = 0.0)
def draw(frame, add_colorbar):
forecast_frame = forecasts.sel(time = str(delta_time * frame))
plot = forecast_frame.plot.pcolormesh(ax=ax, x='longitude', y='latitude', add_colorbar=add_colorbar, cmap='viridis')
title = 'Air Quality Forecast - {0:s}'.format(date_index[frame - 1].strftime('%Y-%m-%d %H:%M:%S'))
ax.set_title(title)
return plot
def init():
return draw(1, add_colorbar=True)
def animate(frame):
plt.close()
return draw(frame + 1, add_colorbar=False)
fig, ax = create_figure()
ani = animation.FuncAnimation(fig, animate, frames=96, interval=500, blit=False, init_func=init, repeat=False)
HTML(ani.to_jshtml())