CAMS Regional: European Air Quality Forecast Optimised at Observation Sites data documentation

"""Demonstration code to download and plot air-quality point forecasts from
   the CAMS Atmosphere Data Store"""
 
import os
import sys
import yaml
import json
import cdsapi
import zipfile
import hashlib
import pandas as pd
from math import ceil, sqrt
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
 
 
def main(cdsapirc_file=None):
 
    # The data to download
    request = {'variable': ['nitrogen_dioxide', 'ozone',
                            'particulate_matter_2.5um',
                            'particulate_matter_10um'],
               'country': 'netherlands',
               'type': ['raw', 'mos_optimised'],
               'leadtime_hour': ['0-23', '24-47'],
               'year': ['2024'],
               'month': ['01'],
               'day': ['27', '28', '29', '30', '31'],
               'include_station_metadata': 'yes',
               'format': 'zip'}
 
    # The stations to plot. If None, plot them all.
    #stations = None
    stations = ['NL00014']
 
    # Plotting preferences
    style = {
        'type': {
            'raw': {'color': 'tab:blue'},
            'mos': {'color': 'tab:orange'}
        },
        'leadtime_day': {
            0: {'linestyle': 'solid'},
            1: {'linestyle': 'dashed'}
        }
    }
 
    # Download the data. Use a filename that depends on the request so we
    # don't have to re-download if the data already exists.
    data_file = data_filename(request)
    if not os.path.exists(data_file):
        get_data(request, data_file, cdsapirc_file)
 
    # Read the data
    station_data, data = read_data(data_file, stations)
 
    # Make a plot for each station
    for station in data.station_id.unique():
 
        # Extract station metadata for just this site. If there's more than one
        # entry we take the latest one that's valid within this time period
        sdata = station_data.loc[
            (station_data.id == station) &
            (station_data.date_start <= data.datetime.iloc[-1]) &
            (station_data.date_end >= data.datetime.iloc[0])
        ]
        assert len(sdata), 'No metadata for site?'
        sdata = sdata.iloc[-1, :]
 
        # Extract air quality data for just this site
        adata = data[data.station_id == station]
 
        if len(adata) > 0:
            station_plot(sdata, adata, style)
        else:
            print('No data for ' + station)
 
 
def data_filename(request):
    """Return a data filename containing a hash which depends on the request"""
    hash = hashlib.md5()
    hash.update(json.dumps(request, sort_keys=True).encode())
    return 'data_' + hash.hexdigest() + '.zip'
 
 
def get_data(request, data_file, cdsapirc_ifile):
    """Download requested data from the ADS"""
 
    # Read the login credentials if provided
    if cdsapirc_file:
        with open(cdsapirc_file, 'r') as f:
            credentials = yaml.safe_load(f)
        kwargs = {'url': credentials['url'],
                  'key': credentials['key'],
                  'verify': credentials['url'].startswith('https://ads.')}
    else:
        kwargs = {}
 
    client = cdsapi.Client(**kwargs)
    client.retrieve(
        'cams-europe-air-quality-forecasts-optimised-at-observation-sites',
        request,
        data_file)
 
 
def read_data(data_file, stations):
    """Read the downloaded zip file and return the station metadata and
       concentration data as Pandas DataFrames"""
 
    data = {}
 
    # Loop over zip file contents
    with zipfile.ZipFile(data_file) as zip:
        for name in sorted(zip.namelist()):
            with zip.open(name) as f:
 
                if name.startswith('station_list'):
 
                    # Read the station metadata file
                    date_fmt = '%Y-%m-%d'
                    station_data = pd.read_csv(f, sep=';',
                                               keep_default_na=False)
 
                    # Remove metadata for stations we're not interested in
                    if stations:
                        station_data = station_data[
                            station_data.id.isin(stations)]
 
                    # Set missing end dates to a date far in the future
                    no_end = (station_data.date_end == '')
                    station_data.loc[no_end, 'date_end'] = '2099-01-01'
 
                    # Parse start and end dates into datetime objects
                    station_data.date_start = pd.to_datetime(
                        station_data.date_start,
                        format=date_fmt)
                    station_data.date_end = pd.to_datetime(
                        station_data.date_end,
                        format=date_fmt)
 
                else:
 
                    # Read the data file
                    df = pd.read_csv(f, sep=';', parse_dates=['datetime'],
                                     infer_datetime_format=True)
 
                    # Remove data for stations we're not interested in
                    if stations:
                        df = df[df.station_id.isin(stations)]
 
                    # Get the name of the column containing concentration so we
                    # can group data by raw/mos type
                    data_col = [c for c in df.columns if c.startswith('conc_')]
                    assert len(data_col) == 1
                    data_col = data_col[0]
                    if data_col not in data:
                        data[data_col] = []
 
                    data[data_col].append(df)
 
    # Merge the data frames
    merged_data = None
    for data_col in list(data.keys()):
 
        # Concatenate all times/countries/species
        x = pd.concat(data[data_col])
 
        # Merge raw and mos into combined records
        if merged_data is None:
            merged_data = x
        else:
            merged_data = merged_data.merge(x, how='outer', validate='1:1',
                                            on=[c for c in x.columns
                                                if c != data_col])
 
    return station_data, merged_data
 
 
def station_plot(station, data, style):
    """Make a series of plots for a station - one for each species"""
 
    # Species to plot
    allspecies = data.species.unique()
 
    # Create the figure
    nplotsx = ceil(sqrt(len(allspecies)))
    nplotsy = ceil(len(allspecies) / nplotsx)
    fig = plt.figure(figsize=(nplotsx*8, nplotsy*5))
    fig.subplots_adjust(hspace=0.35)
 
    # Plot each species
    for iplot, species in enumerate(allspecies):
        ax = plt.subplot(nplotsy, nplotsx, iplot + 1)
 
        # Extract data for just this species
        sdata = data[data.species == species]
 
        plot_species(station, sdata, style, ax)
 
    plt.show()
 
 
def plot_species(station, data, style, ax):
    """Make a plot for one species at one station"""
 
    allspecies = data.species.unique()
    assert len(allspecies) == 1
 
    # Plotting both raw and mos-corrected or just one?
    types = [t for t in ['raw', 'mos']
             if f'conc_{t}_micrograms_per_m3' in data.columns]
 
    # Plot different lines for each post-processing type and lead time day
    for type in types:
        lead_time_day = data.lead_time_hour // 24
        for day in lead_time_day.unique():
            dt = data[lead_time_day == day]
 
            ax.plot(dt.datetime, dt[f'conc_{type}_micrograms_per_m3'],
                    label=f'{type} forecast D+{day}',
                    **prefs(type, day, style))
 
    # Nicer date labels
    ax.xaxis.set_major_formatter(mdates.DateFormatter('%H\n%a %d'))
 
    ax.set_ylabel('$\mu$g / m$^3$')
 
    date_range = 'from ' + ' to '.join(data.datetime.iat[i].strftime('%Y-%m-%d')
                                       for i in [0, -1])
    ax.set_title(
        ('{species} at {station} (lat={lat}, lon={lon}, altitude={alt}m)\n'
         '{dates}').format(species=allspecies[0],
                           station=station.id,
                           lat=station.lat,
                           lon=station.lon,
                           alt=station.alt,
                           dates=date_range))
 
    ax.legend()
 
 
def prefs(type, leadtime, style):
    """Return pyplot.plot() keyword arguments for given type and leadtime"""
 
    return {**style.get('type', {}).get(type, {}),
            **style.get('leadtime_day', {}).get(leadtime, {})}
 
 
if __name__ == '__main__':
 
    # The ADS credentials can be passed as an argument if they're not stored in
    # the default location
    cdsapirc_file = sys.argv[1] if len(sys.argv) > 1 else None
 
    main(cdsapirc_file=cdsapirc_file)