library(sf)
library(leaflet)
library(leaflet.extras) # nécessaire pour addHeatmap()
library(dplyr)
library(webshot2) # to generate pdf with the interactive mapsL’objectif est de produire une carte inspirée de celle de John Snow où les lieux de décès sont représentés par des cercles dont la taille dépend du nombre de décès et où les pompes à eau sont représentées par un symbole distinct. Un fond de carte actuel est utilisé comme support spatial pour visualiser les données historiques.
Les données proviennent du jeu de données géographiques mis à disposition sur le site Web de Robin Wilson. Nous les récupérons sous forme d’un fichier zip. L’URL est:
data_url = "https://rtwilson.com/downloads/SnowGIS_SHP.zip" Le script télécharge une copie locale du fichier uniquement si elle n’existe pas déjà, afin de rendre l’analyse reproductible sans retélécharger inutilement les données.
data_file = "SnowGIS_SHP.zip"
if (!file.exists(data_file)) {
download.file(data_url, data_file, method="auto")
}
# Dézipper dans un sous-dossier
unzip(data_file, exdir = "SnowGIS_SHP")Chaque .shp vient avec ses fichiers compagnons (.dbf, .prj, .shx) qui doivent rester dans le même dossier.
Le document dézippé contient un dossier *_MACOSX* et un dossier SnowGIS_SHP. Le premier, vu son nom, est certainement fait pour les utilisateurs de MAC. J’utilise donc, pour ma part, les documents shp du second dossier.
Avec le package sf, on lit les shapefile:
library(sf)
deaths <- st_read("SnowGIS_SHP/SnowGIS_SHP/Cholera_Deaths.shp")## Reading layer `Cholera_Deaths' from data source
## `C:\Users\foucteau\Documents\8. Encadrements & formations\MOOC-Recherche_reproductible\training_directory\mooc-rr\module3\exo3\SnowGIS_SHP\SnowGIS_SHP\Cholera_Deaths.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 250 features and 2 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 529160.3 ymin: 180857.9 xmax: 529655.9 ymax: 181306.2
## Projected CRS: OSGB36 / British National Gridpumps <- st_read("SnowGIS_SHP/SnowGIS_SHP/Pumps.shp")## Reading layer `Pumps' from data source
## `C:\Users\foucteau\Documents\8. Encadrements & formations\MOOC-Recherche_reproductible\training_directory\mooc-rr\module3\exo3\SnowGIS_SHP\SnowGIS_SHP\Pumps.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 8 features and 1 field
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 529183.7 ymin: 180660.5 xmax: 529748.9 ymax: 181193.7
## Projected CRS: OSGB36 / British National GridLes fichiers sont fournis en British National Grid (EPSG:27700). Ils
sont reprojetés en WGS84 (EPSG:4326), format attendu par
leaflet pour être superposés à un fond OpenStreetMap.
library(sf)
deaths_wgs84 <- st_transform(deaths, crs = 4326)
pumps_wgs84 <- st_transform(pumps, crs = 4326)ON utilise le package leaflet, qui charge les tuiles OpenStreetMap :
icone_pompe <- makeAwesomeIcon(
icon = "tint",
library = "fa",
markerColor = "blue",
iconColor = "white"
)
leaflet() %>%
addProviderTiles(providers$OpenStreetMap) %>% # OSM explicite, pas de watermark
addCircleMarkers(data = deaths_wgs84,
radius = ~Count,
color = "red",
fillOpacity = 0.7,
label = ~paste("Décès :", Count)) %>%
addAwesomeMarkers(data = pumps_wgs84,
icon = icone_pompe,
label = "Pompe")La carte montre une concentration importante des décès autour du secteur de Broad Street, à proximité d’une pompe à eau. Cette représentation ne démontre pas à elle seule une relation causale, mais elle permet de visualiser spatialement l’association qui a rendu la carte de John Snow célèbre.
Comme suggéré dans l’énoncé de l’exercice, il est possible de réprésenter les densitées de décès sur la carte sous forme d’une heatmap :
leaflet() %>%
addProviderTiles(providers$OpenStreetMap) %>%
addHeatmap(data = deaths_wgs84,
intensity = ~Count,
blur = 20,
max = 0.05,
radius = 15) %>%
addAwesomeMarkers(data = pumps_wgs84,
icon = icone_pompe,
label = "Pompe")Calcul de la distance entre chaque décès et chaque pompe : L’idée est de montrer statistiquement que la majorité des décès est plus proche de la pompe de Broad Street que de toute autre pompe.
D’abord, je donne un nom aux pompes selon la rue à laquelle elles appartiennent (ou le lieu d’intérêt le plus proche si le nom de la rue n’est pas visible):
noms_pompes <- c(
"1" = "Broad Street",
"2" = "Liberty",
"3" = "Ramillies Place",
"4" = "Dean Street",
"5" = "Rupert Street",
"6" = "Bridle Lane",
"7" = "Piccadilly Circus",
"8" = "Warwick Street"
)Je calcule ensuite les distances entre chaque décès et chaque pompe avec identification de la pompe la plus proche :
distances <- st_distance(deaths_wgs84, pumps_wgs84)
deaths_wgs84$pompe_proche <- apply(distances, 1, which.min)Je calcule des statistiques de décès par pompe (nombre total et pourcentage sur le total des décès)
stats_pompes <- deaths_wgs84 %>%
st_drop_geometry() %>%
group_by(pompe_proche) %>%
summarise(total_deces = sum(Count)) %>%
mutate(
pourcentage = round(total_deces / sum(total_deces) * 100, 1),
nom_pompe = noms_pompes[as.character(pompe_proche)]
)Je joins les statistiques avec les données des pompes:
pumps_wgs84$pompe_id <- as.integer(seq_len(nrow(pumps_wgs84)))
stats_pompes$pompe_proche <- as.integer(stats_pompes$pompe_proche)
pumps_wgs84 <- left_join(pumps_wgs84, stats_pompes, by = c("pompe_id" = "pompe_proche"))
# Remplacer les NA pour les pompes sans décès associés
pumps_wgs84$total_deces <- ifelse(is.na(pumps_wgs84$total_deces), 0, pumps_wgs84$total_deces)
pumps_wgs84$pourcentage <- ifelse(is.na(pumps_wgs84$pourcentage), 0, pumps_wgs84$pourcentage)
pumps_wgs84$nom_pompe <- ifelse(is.na(pumps_wgs84$nom_pompe),
noms_pompes[as.character(pumps_wgs84$pompe_id)],
pumps_wgs84$nom_pompe)Je défnis une palette de couleurs par pompe :
palette_pompes <- colorFactor(palette = "Set1", domain = deaths_wgs84$pompe_proche
)et génère la carte interactive :
leaflet() %>%
addProviderTiles(providers$OpenStreetMap) %>%
addCircleMarkers(data = deaths_wgs84,
radius = ~Count,
color = ~palette_pompes(pompe_proche),
fillOpacity = 0.7,
label = ~paste("Décès :", Count)) %>%
addAwesomeMarkers(data = pumps_wgs84,
icon = icone_pompe,
label = pumps_wgs84$label_carte) %>%
addLegend(position = "bottomright",
pal = palette_pompes,
values = deaths_wgs84$pompe_proche,
title = "Pompe la plus proche",
labFormat = labelFormat(
transform = function(x) noms_pompes[as.character(x)]
))Cette dernière représentation est la plus visuelle pour se rendre compte du lien entre la pompe de broad street et les morts du choléra en 1854 !
Tentative pour un rendu correct en pdf :
library(sf)
# Créer des triangles comme polygones sf autour de chaque pompe
create_triangle <- function(lon, lat, size = 0.0003) {
matrix(c(
lon, lat + size, # sommet haut
lon + size, lat - size/2, # sommet bas droite
lon - size, lat - size/2, # sommet bas gauche
lon, lat + size # fermer le polygone
), ncol = 2, byrow = TRUE)
}
# Carte
leaflet() %>%
addTiles() %>%
# Décès : cercles pleins colorés par pompe la plus proche,
# taille proportionnelle au nombre de victimes
addCircleMarkers(data = deaths_wgs84,
radius = ~Count,
color = ~palette_pompes(pompe_proche),
fillOpacity = 0.7,
stroke = FALSE,
label = ~paste("Décès :", Count)) %>%
# Pompes : cercles vides avec label permanent
addAwesomeMarkers(data = pumps_wgs84,
label = pumps_wgs84$label_carte,
labelOptions = labelOptions(
noHide = TRUE,
direction = "top",
textOnly = TRUE
)) %>%
# Légende
addLegend(position = "bottomright",
pal = palette_pompes,
values = deaths_wgs84$pompe_proche,
title = "Pompe la plus proche",
labFormat = labelFormat(
transform = function(x) noms_pompes[as.character(x)]
))