Open data - Pollution de l'air - Partie 1

Lig’Air - Open data

Lig’Air est une Association Agréée pour la Surveillance de la Qualité de l’Air et fait partie de la Fédération ATMO France.
Elle a pour rôles la surveillance de la qualité de l’air sur les 6 départements de la région Centre-Val de Loire (Cher, Eure-et-Loir, Indre, Indre-et-Loire, Loir-et-Cher et Loiret), l’information et la diffusion de ses résultats.
Lig’Air met à disposition du public ses données sur son portail Open Data.

Le jeu de données :
Issu des données ouvertes par Lig’Air, le jeu de données analysé ici concerne la pollution de l’air en région Centre - Val de Loire sur une année complète de Novembre 2017 à Novembre 2018.
Le jeu de données comporte les concentrations moyennes journalières issues du réseau fixe des mesures européennes des principaux polluants réglementés dans l’air sur la région Centre-Val de Loire: monoxyde d’azote NO et dioxyde d’azote NO2, particules en suspension PM10, particules en suspension PM2.5, ozone O3, monoxyde de carbone CO (plus d’infos ici).

La dataviz :
Dans cette première partie nous allons cartographier les niveaux d’émissions en polluant jour après jour.
La dataviz est réalisée grâce au logiciel R accompagné des librairies leaflet et rasterpour la cartographie.

Cartographie des émissions d’ozone en région Centre-Val de Loire

Perspectives :
* Etudier les dimensions temporelles et spatiales quant aux niveaux d’emissions de polluants (voir partie 2).
* Etudier les autres polluants du jeu de données : dioxyde d’azote, particules en suspension, etc.
* Evaluer un modèle de dispersion qui tienne compte des conditions météorologiques : vent, température, pression atmosphérique, etc. (modèle de dispersion inversement proportionnel à la distance appliqué dans cet exemple).

Le code :
Toujours dans un esprit Open Data, le code source (langage R) est fourni ci-après pour réutilisation et amélioration.

Code : chargement des librairies et des données
Les données sont téléchargées directement depuis l’API proposée par Lig’Air (voir la requête nommée query).

# Chargement des librairies
library(httr)
library(jsonlite)
library(tidyverse)
library(knitr)
library(raster)
library(gstat)
library(sp)
library(leaflet)


# Chargement du jeu de donnees depuis le portail Open Data de Lig'Air
query <- "https://services1.arcgis.com/HzzPcgRsxxyIZdlU/arcgis/rest/services/mes_centre_val_de_loire_journalier_poll_princ_1/FeatureServer/0/query?where=1%3D1&outFields=*&outSR=4326&f=json&resultType=standard" # resultType au format standard pour étendre la limite autorisée du transfert de données 
get_data <- httr::GET(query)
get_data_text <- httr::content(get_data, "text")
get_data_json <- jsonlite::fromJSON(get_data_text, flatten = TRUE)
data <- as.data.frame(get_data_json$features)


# Simplification des noms de colonne pour usage ulterieur
colnames(data) <- sub("attributes.", "", colnames(data))
colnames(data) <- sub("geometry.", "", colnames(data))

Code : exploration du jeu de données
Quelles sont les villes et les stations équipées de capteurs ?
Quels sont les paramètres mesurés ? Combien de mesures disponibles dans le jeu de données ?

# Calcul du nombre de station par commune
knitr::kable(apply(table(data$nom_com, data$nom_station)>0, 1, sum), caption = "Nombre de station par commune")

Table 1: Nombre de station par commune

	x
Blois	1
Bourges	2
Chanceaux-sur-Choisille	1
Chartres	2
Châteauroux	1
Dreux	1
Faverolles-en-Berry	1
Issoudun	1
Joué-lès-Tours	1
Lucé	1
Marigny-les-Usages	1
Montargis	2
Montierchaume	1
Orléans	2
Oysonville	1
Saint-Jean-de-Braye	1
Saint-Rémy-sur-Avre	1
Tours	2
Verneuil	1

# Calcul du nombre d'observations par commune et type de polluant
knitr::kable(table(data$nom_com, data$nom_polluant), caption = "Nombre d'observations par commune et polluant")

Table 1: Nombre d’observations par commune et polluant

	CO	NO	NO2	O3	PM10	PM2.5
Blois	0	367	367	367	367	0
Bourges	0	734	734	367	734	0
Chanceaux-sur-Choisille	0	0	0	367	0	0
Chartres	0	34	34	367	63	34
Châteauroux	0	365	365	365	365	0
Dreux	0	367	367	367	367	0
Faverolles-en-Berry	0	0	0	366	0	0
Issoudun	0	367	367	0	0	367
Joué-lès-Tours	0	367	367	367	0	367
Lucé	0	308	308	0	308	308
Marigny-les-Usages	0	0	0	367	0	0
Montargis	0	396	396	367	367	396
Montierchaume	0	0	0	367	0	0
Orléans	0	731	731	364	731	0
Oysonville	0	0	0	367	0	0
Saint-Jean-de-Braye	0	367	367	0	0	367
Saint-Rémy-sur-Avre	0	367	367	0	367	0
Tours	367	734	734	0	734	63
Verneuil	0	0	0	367	0	367
Code : cartographie
Dans cette section nous no	us foc	alison	s sur	le pol	luant O	zone (O3). La cartographie est réalisée pour un jour donné par l’appel de la fonction daily_map, développée ici.
A partir des concentration	s rele	vées d	ans le	s stat	ions, n	ous allons appliquer un modèle (simple) de dispersion inversement proportionnelle à la distance (idw), afin de couvrir l’ensemble territoire de la région Centre-Val de Loire par interpolation.

# Initialisation des parametres utiles
## Chargement de la carte de france et extraction de la region Centre - Val de Loire (rCVL)
fr <- raster::getData("GADM", country="france", level=1)
rCVL <- fr[match(toupper("centre"), toupper(fr$NAME_1)),]

## Analyse d'un des polluants du jeu de donnees : choix de O3 car régulièrement mesuré dans la majorité des stations
polluant <- "O3"

## Definition des echelles et unité de representation du polluant
scale <- as.numeric(unlist(data %>%
                             dplyr::filter(nom_polluant == polluant) %>%
                             dplyr::select(valeur)))

breaks <- pretty(scale, n = 6)
unit <- as.character(unlist(data %>%
                            dplyr::filter(nom_polluant == polluant) %>%
                            dplyr::select(unite)))[1]


# Creation de la fonction 'daily_map' : cartographie pour un jour donne
daily_map <- function(date) {
  
  ## Initialisation du raster
  frame <- data %>%
    dplyr::arrange(date_debut) %>%
    dplyr::filter(statut_valid == 1, # utilisation des données exploitables (statut 0 signifie non exploitable)
                  nom_polluant == polluant,
                  date_debut == date) %>%
    dplyr::select(valeur, x, y)

    frame.xy <- frame[c("x", "y")]
    coordinates(frame.xy ) <- ~x+y
    
    ## Initialisation de la grille (au format rectangle selon les coordonnees de la region CVL)
    grd <- expand.grid(x = seq(from = xmin(rCVL), to = xmax(rCVL), by = 0.01), 
                       y = seq(from = ymin(rCVL), to = ymax(rCVL), by = 0.01))
    coordinates(grd) <- ~x+y
    gridded(grd) <- TRUE
    
    ## Interpolation des emissions de polluant sur toute le surface du raster (maillage)
    interpolation <- gstat::idw(formula = frame$valeur ~ 1, locations = frame.xy, newdata = grd)  # application d'un modele de dispersion simple 'idw' (inverse distance weighting)
    raster <- raster::raster(interpolation, "var1.pred")
    
    ## Adaptation du raster (format rectangle) en un polygone adaptes aux coordonnees exactes de la region CVL
    raster.rCVL <- raster::mask(raster, rCVL)
    crs(raster.rCVL) <- sp::CRS("+init=epsg:4326")

    ## Initialisation des couleurs et du fond de carte
    pal <- leaflet::colorNumeric("RdYlBu", breaks, na.color = "transparent", reverse = FALSE)
    map <- leaflet::leaflet(data) %>%
      leaflet::fitBounds(lng1 = xmin(rCVL), lat1 = ymin(rCVL), lng2 = xmax(rCVL), lat2 = ymax(rCVL)) %>%
      leaflet::addTiles()
    
    ## Integration des donnees de maillage sur le fond de carte
    map <- leaflet::addRasterImage(map, x=raster.rCVL, opacity = 0.5) %>%
      leaflet::addLegend(pal = pal, values = breaks, title = paste0("Concentration ", polluant, " [", unit, "]"),
                labFormat = labelFormat(transform = function(x) sort(x, decreasing = TRUE))) %>%
      leaflet::addControl(paste0("Date : ", date), position = "topleft")
    
    ## Integration des contours sur le fond de carte
    if (sd(frame$valeur) != 0) {
      contour.rCVL <- raster::rasterToContour(raster.rCVL)
      map <- leaflet::addPolylines(map=map, data=contour.rCVL, fillOpacity=0.5, color = "grey95", fillColor = "transparent", opacity=10, weight=0.8, dashArray="4 1")
    }
    
    ## Sauvegarde
    htmlwidgets::saveWidget(map, paste0("map_", date, ".html"), selfcontained = FALSE)

}


# Creation des cartes par l'execution de la fonction 'daily_map'
selection.date <- as.Date(as.character(data %>%
    dplyr::filter(statut_valid == 1, # utilisation des données exploitables (statut 0 signifie non exploitable)
                  nom_polluant == polluant) %>%
    dplyr::summarise(max(date_debut))))