Tutorial per lo scaricamento dei dati sulla qualità dell'aria (CAMS) da Copernicus¶

Questo tutorial permette di scaricare i dati sulla qualità dell'aria (CAMS) dal sito di Copernicus

In [ ]:
import cdsapi
import geopandas as gpd
import xarray as xr
import os
import yaml
import hashlib
import plotly
from shapely import geometry
from datetime import datetime, timedelta
from math import dist
from scipy.stats import zscore
import plotly.graph_objects as go

plotly.offline.init_notebook_mode()
In [2]:
def get_dates(start_month, n_month=1):
    start_date = datetime.strptime(start_month, '%Y-%m')
    end_date = start_date + timedelta(days=n_month*31)
    return start_date.strftime('%Y-%m-01'), end_date.strftime('%Y-%m-01')

def sanitize_string(string_in):
    delete_chars = '\\/:*?"<>|,'
    for char in delete_chars:
        string_in = string_in.replace(char,'')
    ugly_chars = " '-" 
    for char in ugly_chars:
        string_in = string_in.replace(char,'_')
    return string_in

Per ottenere le credenziali per utilizzare cdsapi, seguire le istruzioni a questo link: https://ads.atmosphere.copernicus.eu/api-how-to

Per selezionare tutti i parametri di interesse (produce un json in fondo): https://ads.atmosphere.copernicus.eu/cdsapp#!/dataset/cams-europe-air-quality-forecasts?tab=form

Noi abbiamo analizzato i dati della Sardegna partendo dai confini amministrativi forniti dall'ISTAT: https://www.istat.it/storage/cartografia/confini_amministrativi/non_generalizzati/2024/Limiti01012024.zip

Abbiamo anche selezionato il mese di gennaio 2024 con tutte le ore del giorno

In [3]:
regione_sel = 'Sardegna'
month = '2024-06'
hours = [f'{i:02}:00' for i in range(0,24)]
In [4]:
start_date, end_date = get_dates(month, n_month = 1)
month_out = month.replace("-","")
regione_sel_out = sanitize_string(regione_sel.lower())
cams_data_fname = f"{os.path.join(os.path.expanduser('~'), 'ILAB_DATA', 'COPERNICUS')}/NC/data_{regione_sel_out}_{month_out}.nc"
In [5]:
comuni = gpd.read_parquet(os.path.join(os.path.expanduser('~'), 'ILAB_DATA','ISTAT', 'DATA', 'comuni_pop.parquet'))
regione = comuni[comuni['cod_reg']==20].copy()
regione_bb = list(regione['geometry'].to_crs(4326).total_bounds)
bb = [regione_bb[3],regione_bb[0],regione_bb[1],regione_bb[2]]
In [6]:
with open ('../../.adsapirc', 'r') as f:
    credentials= yaml.safe_load(f)
c = cdsapi.Client(url=credentials['url'], key=credentials['key'])

Scarico i dati CAMS (raster)

In [ ]:
c.retrieve(
'cams-europe-air-quality-forecasts',
{
    'model': 'ensemble',
    'date': f'{start_date}/{end_date}',
    'format': 'netcdf',
    'variable': [
        'ammonia', 'carbon_monoxide', 'nitrogen_dioxide',
        'ozone', 'particulate_matter_10um', 'particulate_matter_2.5um',
        'pm10_wildfires', 'pm2.5_total_organic_matter', 'residential_elementary_carbon',
        'sulphur_dioxide', 'total_elementary_carbon',
    ],
    'level': '0',
    'type': 'analysis',
    'time': hours,
    'leadtime_hour': '0',
    'area': bb
},cams_data_fname)

Apro il raster e lo trasformo in un dataframe

In [9]:
ds = xr.open_dataset(cams_data_fname)
df = ds.to_dataframe().reset_index()
ds.close()

Rinomino le colonne delle variabili eliminando '_conc'

In [10]:
df.columns = [col.replace('_conc','') for col in df.columns]
In [ ]:
df.head()

Seleziono le colonne che contengono le variabili, sulle quali verranno calcolati i valori medi

In [12]:
variables = list(df.columns)
variables = [element for element in variables if element not in ['longitude','latitude','level','time']]
aggs_dict = {element: 'mean' for element in variables}

Ottengo un valore di tempo assoluto

In [13]:
timestart = datetime.strptime(start_date, '%Y-%m-%d')
df['time_abs'] = timestart + df['time']
df['date'] = df['time_abs'].dt.date

df['time_abs'] = df['time_abs'].astype('string')
df['time'] = df['time'].astype('string')
df['date'] =  df['date'].astype('string')
In [ ]:
df.head()

Ora trasformiamo le celle del file raster in vettoriale¶

Per calcolare il passo della griglia (le dimensioni della cella), seleziono le prime due celle relative al primo intervallo temporale

In [15]:
lons = df[df['time']==f'0 days {hours[0]}:00'].sort_values(by=['latitude','longitude']).head(2)['longitude'].values
lats = df[df['time']==f'0 days {hours[0]}:00'].sort_values(by=['latitude','longitude']).head(2)['latitude'].values
distance = dist([lons[1]], [lons[0]])

Definisco una chiave univoca per ogni cella della griglia, data dal passo e dai valori di latitudine e longitudine

In [16]:
def get_id_cell(row, distance_in):
    grid_distance = str(distance_in)[:3].replace('.','_')
    id_cell = f'{grid_distance}__{row["latitude"]}_{row["longitude"]}'
    result = hashlib.md5(id_cell.encode()).hexdigest()
    return result

df['id_cell'] = df.apply(get_id_cell, distance_in=distance, axis='columns')

Seleziono un dataframe relativo al primo valore temporale immesso, che userò come base per le geometrie (non ha senso rifare le geometrie per ogni ora del giorno)

In [17]:
df_geom = df[df['time']==f'0 days {hours[0]}:00'].copy()
df_geom = df_geom[['longitude','latitude','id_cell']].copy()

Creo il campo geometria per ogni cella e lo trasformo in un geodataframe. Così otteniamo le geometrie di tutte le celle della Sardegna

In [18]:
def make_cell(row, distance_in):
    xmin = row['longitude'] - (distance_in/2)
    ymin = row['latitude'] - (distance_in/2)
    xmax = row['longitude'] + (distance_in/2)
    ymax = row['latitude'] + (distance_in/2)
    geom = geometry.Polygon(((xmin,ymin), (xmin,ymax), (xmax,ymax), (xmax,ymin)))
    return geom

df_geom['geometry'] = df_geom.apply(make_cell, distance_in=distance, axis='columns')

gdf = gpd.GeoDataFrame(df_geom, geometry=df_geom['geometry'])
gdf = gdf[['id_cell', 'geometry']].set_crs(4326)

Uniamo l'info sulla geometria ai dati CAMS precedentemente scaricati

In [19]:
gdf_out = gdf.merge(df, on='id_cell')
In [20]:
gdf_out.to_parquet(os.path.join(os.path.expanduser('~'), 'ILAB_DATA', 'COPERNICUS','DATA', f'cams_grid_{regione_sel_out}_{month_out}.parquet'))
In [ ]:
gdf_out.head()
In [24]:
# # cella a Roma
# fig = go.Figure()
# for var in variables:
#     temp = gdf_out[(gdf_out['id_cell']=='53d7e08a8fae6e9b00f5f516945b6a37')].copy()# & (gdf_out['date']=='2024-01-01')
#     temp['zscore'] = temp[[var]].apply(zscore)
#     fig.add_trace(go.Scatter(x=temp['time_abs'], y=temp['zscore'], name=var))

# fig.show()

Visualizziamo i dati¶

Vediamo l'andamento mensile di tutte le variabili scaricate per una specifica cella sul comune di Cagliari

In [ ]:
fig = go.Figure()
for var in variables:
    temp = gdf_out[(gdf_out['id_cell']=='54cfb985fad9555c04e4b4135a471259')].copy()# & (gdf_out['date']=='2024-01-01')
    temp['zscore'] = temp[[var]].apply(zscore)
    fig.add_trace(go.Scatter(x=temp['time_abs'], y=temp['zscore'], name=var))

fig.show()

Possiamo anche calcolare i valori medi comunali¶

Eseguo l'intersezione tra i poligoni dei comuni e le celle della griglia (facciamo una sola volta l'intersezione sulle geometrie e poi riaggiungiamo tutte le info dei comuni e di CAMS che sono ripetute per giorno e per ora)

In [25]:
intersection = gdf.overlay(regione[['geometry','pro_com_t','comune']].to_crs(4326), how="intersection")
intersection = intersection[['pro_com_t','id_cell']].copy()
m = intersection.merge(regione[['geometry','pro_com_t','comune']], on = 'pro_com_t')
m = m.merge(df, on='id_cell')

Eseguo il groupby per ogni comune per ottenere i valori mediati delle grandezze

In [26]:
g = m.groupby(['pro_com_t','comune'], as_index=False).agg(aggs_dict)
regione_stats = regione.merge(g)
regione_stats = gpd.GeoDataFrame(regione_stats, geometry=regione_stats['geometry'])
In [ ]:
regione_stats.to_parquet(os.path.join(os.path.expanduser('~'), 'ILAB_DATA', 'DATA', f'cams_comuni_{regione_sel_out}_{month_out}.parquet'))
In [ ]:
g.shape, regione_stats.shape
In [ ]:
regione_stats[regione_stats['comune']=='Cagliari'].head()
In [ ]: