Qui troverai l'ultimo passo "dell'avventura" iniziata con questo articolo in cui riporto l'andamento del PM10 nel nord Italia nell'ultimo decennio. Il come ho proceduto all'analisi è riportato in quest'altro articolo.

Lo scopo dell'articolo che stai per leggere è visualizzare su mappa le aree di analisi, cosa che fino ad ora non avevo mostrato.

Tutto l'iter è stato bello lungo ed è per questo che l'ho dovuto suddividere in ben tre differenti articoli.

Prima di iniziare¶

Librerie¶

Come per tutti gli articoli di PyGIS-Blog inizio con l'elenco delle librerie usate. Una menzione speciale è per la funzione xarray.open_mfdataset perchè mi ha consentito di velocizzare tantissimo la lettura di tutti i file .nc grazie al nativo supporto a dask.delayed.

Devo approfondire Dask appena ho tempo e Dask Delayed sembra un buon punto di partenza pratico su cui lavorare per il mio solito training on the job :)

In [1]:
from pathlib import Path
import xarray as xr
import geopandas as gpd
import matplotlib.pyplot as plt
from matplotlib import patheffects

Fonti dati¶

Come già spiegato nei precedenti articoli di questa involontaria miniserie, i dati usati sono quelli distributi dal Copernicus Atmosphere Monitoring Service. In particolare ho usato dati con una copertura temporale che va dal 2013 a marzo 2024.

In [2]:
nc_data_path = Path("/home/max/Desktop/pianura_padana/processed/netcdf")

nc_files = list(nc_data_path.glob("*.nc"))

nc_files
Out[2]:
[PosixPath('/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2021-forecasts.nc'),
 PosixPath('/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2017-reanalyses.nc'),
 PosixPath('/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2020-reanalyses.nc'),
 PosixPath('/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2014-reanalyses.nc'),
 PosixPath('/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2013-reanalyses.nc'),
 PosixPath('/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2016-reanalyses.nc'),
 PosixPath('/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2023-forecasts.nc'),
 PosixPath('/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2019-reanalyses.nc'),
 PosixPath('/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2015-reanalyses.nc'),
 PosixPath('/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2018-reanalyses.nc'),
 PosixPath('/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2022-forecasts.nc'),
 PosixPath('/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2024-forecasts.nc')]
In [3]:
main_city_path = Path.cwd().joinpath('data', 'citta_significative.shp')
main_city = gpd.read_file(main_city_path)

target_zones_boundaries = Path.cwd().parent.parent.parent.parent.joinpath('sample_data').joinpath('ProvCM01012023_g').joinpath('ProvCM01012023_g_WGS84.shp')

target_zones = gpd.read_file(target_zones_boundaries)
target_zones = target_zones[target_zones['COD_REG'].isin(list(range(1, 12)))]
target_zones = target_zones.to_crs(4326).sort_values('DEN_UTS')

figsize_w = 10
figsize_h = figsize_w / 1
plot_dpi = 300

polygon_facecolor = "none"
polygon_edgecolor = "#76CFEF"
polygon_linewidth = 1.0

marker_edgecolor = "black"
marker_facecolor = "red"
marker_type = "*"
marker_s = 250
marker_label_facecolor = "white"
marker_label_alpha = 0.75
marker_font_color = "black"
marker_font_size = 8
marker_x_offset = 0.05
marker_y_offset = -0.05

path_effects_linewidth = 3
path_effects_foreground = "white"

raster_cmap = "autumn"

L'area geografica è sempre quella dei precedenti due articoli, la riporto comunque di seguito per completezza.

In [4]:
fig1, ax1 = plt.subplots(figsize=(figsize_w, figsize_h), dpi=plot_dpi)

target_zones.plot(
    ax=ax1,
    facecolor=polygon_facecolor, 
    edgecolor=polygon_edgecolor,
    linewidth=polygon_linewidth
)
for _index, _row in main_city.iterrows():
    coordinates = _row.geometry.xy
    ax1.scatter(
        *coordinates,
        s=marker_s,
        marker=marker_type,
        facecolor=marker_facecolor,
        edgecolor=marker_edgecolor
    )
    label = plt.text(
        x=coordinates[0][0] + marker_x_offset,
        y=coordinates[1][0] + marker_y_offset,
        s=_row.citta,
        fontdict=dict(color=marker_font_color, size=marker_font_size),
        bbox=dict(facecolor=marker_label_facecolor, alpha=marker_label_alpha),
        path_effects=[patheffects.withStroke(linewidth=path_effects_linewidth, foreground=path_effects_foreground)]
    )
    label.set_path_effects([patheffects.withStroke(linewidth=path_effects_linewidth, foreground=path_effects_foreground)])
plt.show()
No description has been provided for this image

Contenuti¶

  • Analisi sul singolo file
  • Analisi su tutti i file
  • Conclusione

Analisi sul singolo file¶

Lettura del .nc¶

Prima di iterare tutto il processo sull'intero arco temporale ho preferito effettuare delle verifiche su un solo anno, anche per prendere meglio confidenza con i dati.

In [5]:
single_dataset = xr.open_dataset(
    filename_or_obj='/home/max/Desktop/pianura_padana/processed/netcdf/particulate_matter_10um-2013-reanalyses.nc',
    engine="netcdf4",
    decode_coords="all",
)

single_dataset
Out[5]: