See also the animated map.
From the official data by Santé Publique France, we spatially smooth the decease (produced by SPF at the département scale) and normalize by a similarly smoothed population grid. For that we use the {btb} package.
# Carto décès COVID 19 France
# avec lissage
# sources -----------------------------------------------------------------
# https://www.data.gouv.fr/fr/datasets/donnees-hospitalieres-relatives-a-lepidemie-de-covid-19/
fichier_covid <- "donnees/covid.csv"
url_donnees_covid <- "https://www.data.gouv.fr/fr/datasets/r/63352e38-d353-4b54-bfd1-f1b3ee1cabd7"
# https://www.insee.fr/fr/statistiques/2012713#tableau-TCRD_004_tab1_departements
fichier_pop <- "donnees/pop.xls"
url_donnees_pop <- "https://www.insee.fr/fr/statistiques/fichier/2012713/TCRD_004.xls"
# Adminexpress : à télécharger manuellement
# https://geoservices.ign.fr/documentation/diffusion/telechargement-donnees-libres.html#admin-express
#aex <- "donnees/1_DONNEES_LIVRAISON_2019-03-14/"
aex <- path_expand("~/Downloads/ADMIN-EXPRESS_2-2__SHP__FRA_2020-02-24/ADMIN-EXPRESS/1_DONNEES_LIVRAISON_2020-02-24")
# config ------------------------------------------------------------------
library(tidyverse)
library(httr)
library(fs)
library(sf)
library(readxl)
library(janitor)
library(glue)
library(tmap)
library(grid)
library(classInt)
# + btb, raster, fasterize, plyr
rayon <- 100000 # distance de lissage (m)
pixel <- 10000 # résolution grille (m)
force_download <- TRUE # retélécharger même si le fichier existe et a été téléchargé aujourd'hui ?
#' Kernel weighted smoothing with arbitrary bounding area
#'
#' @param df sf object (points)
#' @param field weight field in the df
#' @param bandwidth kernel bandwidth (map units)
#' @param resolution output grid resolution (map units)
#' @param zone sf study zone (polygon)
#' @param out_crs EPSG (should be an equal-area projection)
#'
#' @return a raster object
#' @import btb, raster, fasterize, dplyr, plyr, sf
lissage <- function(df, field, bandwidth, resolution, zone, out_crs = 3035) {
if (st_crs(zone)$epsg != out_crs) {
message("reprojecting data...")
zone <- st_transform(zone, out_crs)
}
if (st_crs(df)$epsg != out_crs) {
message("reprojecting study zone...")
df <- st_transform(df, out_crs)
}
zone_bbox <- st_bbox(zone)
# grid generation
message("generating reference grid...")
zone_xy <- zone %>%
dplyr::select(geometry) %>%
st_make_grid(
cellsize = resolution,
offset = c(
plyr::round_any(zone_bbox[1] - bandwidth, resolution, f = floor),
plyr::round_any(zone_bbox[2] - bandwidth, resolution, f = floor)
),
what = "centers"
) %>%
st_sf() %>%
st_join(zone, join = st_intersects, left = FALSE) %>%
st_coordinates() %>%
as_tibble() %>%
dplyr::select(x = X, y = Y)
# kernel
message("computing kernel...")
kernel <- df %>%
cbind(., st_coordinates(.)) %>%
st_set_geometry(NULL) %>%
dplyr::select(x = X, y = Y, field) %>%
btb::kernelSmoothing(
dfObservations = .,
sEPSG = out_crs,
iCellSize = resolution,
iBandwidth = bandwidth,
vQuantiles = NULL,
dfCentroids = zone_xy
)
# rasterization
message("\nrasterizing...")
raster::raster(
xmn = plyr::round_any(zone_bbox[1] - bandwidth, resolution, f = floor),
ymn = plyr::round_any(zone_bbox[2] - bandwidth, resolution, f = floor),
xmx = plyr::round_any(zone_bbox[3] + bandwidth, resolution, f = ceiling),
ymx = plyr::round_any(zone_bbox[4] + bandwidth, resolution, f = ceiling),
resolution = resolution
) %>%
fasterize::fasterize(kernel, ., field = field)
}
# téléchargement--------------------------------------------------------------
if (!file_exists(fichier_covid) |
file_info(fichier_covid)$modification_time < Sys.Date() |
force_download) {
GET(url_donnees_covid,
progress(),
write_disk(fichier_covid, overwrite = TRUE))
}
if (!file_exists(fichier_pop)) {
GET(url_donnees_pop,
progress(),
write_disk(fichier_pop))
}
# données -----------------------------------------------------------------
covid <- read_csv2(fichier_covid)
# adminexpress prétéléchargé
dep <- read_sf(path(aex, "ADE_2-2_SHP_LAMB93_FR/DEPARTEMENT.shp")) %>%
clean_names() %>%
st_set_crs(2154)
dep_971 <- read_sf(path(aex, "ADE_2-2_SHP_RGAF09UTM20_D971/DEPARTEMENT.shp")) %>%
clean_names() %>%
st_set_crs(5490)
dep_972 <- read_sf(path(aex, "ADE_2-2_SHP_RGAF09UTM20_D972/DEPARTEMENT.shp")) %>%
clean_names() %>%
st_set_crs(5490)
dep_973 <- read_sf(path(aex, "ADE_2-2_SHP_UTM22RGFG95_D973/DEPARTEMENT.shp")) %>%
clean_names() %>%
st_set_crs(2972)
dep_974 <- read_sf(path(aex, "ADE_2-2_SHP_RGR92UTM40S_D974/DEPARTEMENT.shp")) %>%
clean_names() %>%
st_set_crs(2975)
dep_976 <- read_sf(path(aex, "ADE_2-2_SHP_RGM04UTM38S_D976/DEPARTEMENT.shp")) %>%
clean_names() %>%
st_set_crs(4471)
pop <- read_xls(fichier_pop, skip = 2) %>%
clean_names()
# prétraitement -----------------------------------------------------------
# contour métropole
fr <- dep %>%
st_union() %>%
st_sf()
# jointures des données
creer_df <- function(territoire) {
territoire %>%
left_join(pop, by = c("insee_dep" = "x1")) %>%
left_join(
covid %>%
filter(jour == max(jour),
sexe == 0) %>%
group_by(dep) %>%
summarise(deces = sum(dc, na.rm = TRUE),
reanim = sum(rea, na.rm = TRUE),
hospit = sum(hosp, na.rm = TRUE)),
by = c("insee_dep" = "dep")) %>%
st_point_on_surface()
}
covid_geo_pop <- creer_df(dep)
covid_geo_pop_971 <- creer_df(dep_971)
covid_geo_pop_972 <- creer_df(dep_972)
covid_geo_pop_973 <- creer_df(dep_973)
covid_geo_pop_974 <- creer_df(dep_974)
covid_geo_pop_976 <- creer_df(dep_976)
# lissage -----------------------------------------------------------------
# génération des grilles mortalité, hospitalisation et réanimation et population
# et création des grilles pour 100000 habitants
# décès métropole et DOM
d <- covid_geo_pop %>%
lissage("deces", rayon, pixel, fr)
d_971 <- covid_geo_pop_971 %>%
lissage("deces", rayon, pixel, dep_971, 5490)
d_972 <- covid_geo_pop_972 %>%
lissage("deces", rayon, pixel, dep_972, 5490)
d_973 <- covid_geo_pop_973 %>%
lissage("deces", rayon, pixel, dep_973, 2972)
d_974 <- covid_geo_pop_974 %>%
lissage("deces", rayon, pixel, dep_974, 2975)
d_976 <- covid_geo_pop_976 %>%
lissage("deces", rayon, pixel, dep_976, 4471)
# population métropole et DOM
p <- covid_geo_pop %>%
lissage("x2020_p", rayon, pixel, fr)
p_971 <- covid_geo_pop_971 %>%
lissage("x2020_p", rayon, pixel, dep_971, 5490)
p_972 <- covid_geo_pop_972 %>%
lissage("x2020_p", rayon, pixel, dep_972, 5490)
p_973 <- covid_geo_pop_973 %>%
lissage("x2020_p", rayon, pixel, dep_973, 2972)
p_974 <- covid_geo_pop_974 %>%
lissage("x2020_p", rayon, pixel, dep_974, 2975)
p_976 <- covid_geo_pop_976 %>%
lissage("x2020_p", rayon, pixel, dep_976, 4471)
# grilles pour 100000 hab
d100k <- d * 100000 / p
d100k_971 <- d_971 * 100000 / p_971
d100k_972 <- d_972 * 100000 / p_972
d100k_973 <- d_973 * 100000 / p_973
d100k_974 <- d_974 * 100000 / p_974
d100k_976 <- d_976 * 100000 / p_976
# réanimation et hospitalisation métropole uniquement
r <- covid_geo_pop %>%
lissage("reanim", rayon, pixel, fr)
r100k <- r * 100000 / p
h <- covid_geo_pop %>%
lissage("hospit", rayon, pixel, fr)
h100k <- h * 100000 / p
# carto -------------------------------------------------------------------
# décès métropole et DOM
# classification à réutiliser pour les 6 cartes
set.seed(1234)
classes <- classIntervals(raster::values(d100k), n = 5, style = "kmeans", dataPrecision = 0)$brks
# métro et DOM
(carte_d <- tm_layout(title = paste("COVID-19 - France - cumul au", max(covid$jour)),
legend.position = c("left", "bottom"),
frame = FALSE) +
tm_shape(d100k) +
tm_raster(title = glue("décès à l'hôpital
pour 100 000 hab."),
style = "fixed",
breaks = classes,
palette = "viridis",
legend.format = list(text.separator = "à moins de",
digits = 0),
legend.reverse = TRUE) +
tm_shape(dep) +
tm_borders() +
tm_credits(glue("http://r.iresmi.net/
lissage noyau bisquare {rayon / 1000} km sur grille {pixel / 1000} km
classif. kmeans
projections LAEA Europe (métropole) et locales (DOM)
données départementales Santé publique France,
INSEE RP 2020, IGN Adminexpress 2020"),
size = .5,
position = c(.5, .025))
)
tm_971 <- tm_shape(d100k_971, ext = 0.7) +
tm_raster(style = "fixed",
breaks = classes,
palette = "viridis",
legend.show = FALSE) +
tm_shape(dep_971) +
tm_borders() +
tm_layout(frame = FALSE,
bg.color = NA)
tm_972 <- tm_shape(d100k_972, ext = 0.7) +
tm_raster(style = "fixed",
breaks = classes,
palette = "viridis",
legend.show = FALSE) +
tm_shape(dep_972) +
tm_borders() +
tm_layout(frame = FALSE,
bg.color = NA)
tm_973 <- tm_shape(d100k_973) +
tm_raster(style = "fixed",
breaks = classes,
palette = "viridis",
legend.show = FALSE) +
tm_shape(dep_973) +
tm_borders() +
tm_layout(frame = FALSE,
bg.color = NA)
tm_974 <- tm_shape(d100k_974, ext = 0.75) +
tm_raster(style = "fixed",
breaks = classes,
palette = "viridis",
legend.show = FALSE) +
tm_shape(dep_974) +
tm_borders()+
tm_layout(frame = FALSE,
bg.color = NA)
tm_976 <- tm_shape(d100k_976, ext = 0.6) +
tm_raster(style = "fixed",
breaks = classes,
palette = "viridis",
legend.show = FALSE) +
tm_shape(dep_976) +
tm_borders()+
tm_layout(frame = FALSE,
bg.color = NA)
# assemblage
fichier_carto <- glue("resultats/covid_fr_{max(covid$jour)}.png")
tmap_save(carte_d, fichier_carto, width = 900, height = 900, scale = .4,
insets_tm = list(tm_971, tm_972, tm_973, tm_974, tm_976),
insets_vp = list(viewport(x = .1, y = .65, width = .15, height = .15),
viewport(x = .1, y = .58, width = .15, height = .15),
viewport(x = .15, y = .4, width = .35, height = .45),
viewport(x = .9, y = .4, width = .15, height = .15),
viewport(x = .9, y = .5, width = .15, height = .15)))
2 réponses sur « Coronavirus : spatially smoothed decease in France »
[…] by data_admin [This article was first published on r.iresmi.net, and kindly contributed to R-bloggers]. (You can report issue about the content on this page […]
[…] will again use the same COVID-19 dataset. We manually add the overseas départements once we have found the right grid (by trying […]