Open data - Pollution de l'air - Partie 1

Par Rémi Julien | December 14, 2018

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 1
Châteauroux 1
Dreux 1
Faverolles-en-Berry 1
Issoudun 1
Joué-lès-Tours 1
Lucé 1
Marigny-les-Usages 1
Montargis 1
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 0 0 367 11 0
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 364 364 0 0 364
Joué-lès-Tours 0 367 367 367 0 367
Lucé 0 360 360 0 360 360
Marigny-les-Usages 0 0 0 367 0 0
Montargis 0 367 367 367 367 367
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 11
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(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 de la ou des cartes par l'execution de la fonction 'daily_map'
daily_map(as.Date(max(data$date_debut)))
comments powered by Disqus