Catégories
Non classé

COVID-19 decease animation map

Coronavirus decease in France
# Animation carto décès COVID 19 France
# avec lissage

# packages -----------------------------------------------------------------
library(tidyverse)
library(httr)
library(fs)
library(sf)
library(readxl)
library(janitor)
library(glue)
library(tmap)
library(grid)
library(classInt)
library(magick)
# + btb, raster, fasterize, plyr

# 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 <- path_expand("~/Downloads/ADMIN-EXPRESS_2-2__SHP__FRA_2020-02-24/ADMIN-EXPRESS/1_DONNEES_LIVRAISON_2020-02-24")

# config ------------------------------------------------------------------

rayon <- 100000 # distance de lissage (m)
pixel <- 10000 # résolution grille (m)

force_download <- FALSE # 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 (!dir_exists("donnees")) dir_create("donnees")
if (!dir_exists("resultats")) dir_create("resultats")
if (!dir_exists("resultats/animation")) dir_create("resultats/animation")

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)

pop <- read_xls(fichier_pop, skip = 2) %>%
  clean_names()


# prétraitement -----------------------------------------------------------

# contour métropole pour grille de référence
fichier_fr <- "donnees/fr.rds"

if (!file_exists(fichier_fr)) {
  fr <- dep %>%
    st_union() %>%
    st_sf() %>%
    write_rds(fichier_fr)
} else {
  fr <- read_rds(fichier_fr)
}

# jointures des données
creer_df <- function(territoire, date = NULL) {
  territoire %>%
    left_join(pop, by = c("insee_dep" = "x1")) %>%
    left_join(
      covid %>%
        filter(jour == if_else(is.null(date), max(jour), date),
               sexe == 0) %>%
               rename(deces = dc,
                      reanim = rea,
                      hospit = hosp),
      by = c("insee_dep" = "dep")) %>%
    st_point_on_surface()
}

covid_geo_pop <- creer_df(dep)


# lissage -----------------------------------------------------------------
# génération de la dernière grille mortalité
# et création des grilles pour 100000 habitants

# décès métropole 
d <- covid_geo_pop %>%
  lissage("deces", rayon, pixel, fr)


# population métropole et DOM
p <- covid_geo_pop %>%
  lissage("x2020_p", rayon, pixel, fr)

# grilles pour 100000 hab
d100k <- d * 100000 / p


# classification à réutiliser pour les autres cartes
set.seed(1234)
classes <- classIntervals(raster::values(d100k), n = 6, style = "kmeans", dataPrecision = 0)$brks


# animation ---------------------------------------------------------------

image_animation <- function(date) {
  m <- creer_df(dep, date) %>%
    lissage("deces", rayon, pixel, fr) %>%
    magrittr::divide_by(p) %>%
    magrittr::multiply_by(100000) %>%
    tm_shape() +
    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_layout(title = glue("COVID-19 - France métropolitaine - cumul au {date}"),
              legend.position = c("left", "bottom"),
              frame = FALSE) +
    tm_credits(glue("http://r.iresmi.net/
                  lissage noyau bisquare {rayon / 1000} km sur grille {pixel / 1000} km
                  classif. kmeans
                  projection LAEA Europe
                  données départementales Santé publique France,
                  INSEE RP 2020, IGN Adminexpress 2020"),
               size = .5,
               position = c(.5, .025))
  
  tmap_save(m, glue("resultats/animation/covid_fr_{date}.png"),
            width = 800, height = 800, scale = .4,)
}

unique(covid$jour) %>%
  walk(image_animation)

animation <- glue("resultats/deces_covid19_fr_{max(covid$jour)}.gif")

dir_ls("resultats/animation") %>%
  map(image_read) %>%
  image_join() %>%
  #image_scale("500x500") %>%
  image_morph(frames = 1) %>%
  image_animate(fps = 2, optimize = TRUE) %>%
  image_write(animation)


13 réponses sur « COVID-19 decease animation map »

Bonjour, j’ai installé R 3.6.3 (je connais pas R), quand je lance le code avec `Rscript`, j’obtiens :

« `
Linking to ImageMagick 6.9.10.23
Enabled features: fontconfig, freetype, fftw, lcms, pango, webp, x11
Disabled features: cairo, ghostscript, rsvg
Using ‘,’ as decimal and ‘.’ as grouping mark. Use read_delim() for more control.
Parsed with column specification:
cols(
dep = col_character(),
sexe = col_double(),
jour = col_date(format = «  »),
hosp = col_double(),
rea = col_double(),
rad = col_double(),
dc = col_double()
)
New names:
* «  -> …1
* «  -> …2
Warning message:
In st_point_on_surface.sf(.) :
st_point_on_surface assumes attributes are constant over geometries of x
reprojecting data…
reprojecting study zone…
generating reference grid…
Error: Columns 1 and 2 must be named.
Backtrace:

1. └─covid_geo_pop %>% lissage(« deces », rayon, pixel, fr)
2. ├─base::withVisible(eval(quote(`_fseq`(`_lhs`)), env, env))
3. └─base::eval(quote(`_fseq`(`_lhs`)), env, env)
4. └─base::eval(quote(`_fseq`(`_lhs`)), env, env)
5. └─`_fseq`(`_lhs`)
6. └─magrittr::freduce(value, `_function_list`)
7. ├─base::withVisible(function_list[[k]](value))
8. └─function_list[[k]](value)
9. └─global::lissage(., « deces », rayon, pixel, fr)
10. └─`%>%`(…)
11. ├─base::withVisible(eval(quote(`_fseq`(`_lhs`)), env, env))
12. └─base::eval(quote(`_fseq`(`_lhs`)), env, env)
13. └─base::eval(quote(`_fseq`(`_lhs`)), env, env)

xecution halted
« `

Après avoir changer la ligne 78 en :
`as_tibble(.name_repair = ~c(« X », « Y »)) %>%`, j’obtiens :

« `
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
✔ ggplot2 3.3.0 ✔ purrr 0.3.4
✔ tibble 3.0.1 ✔ dplyr 0.8.5
✔ tidyr 1.1.0 ✔ stringr 1.4.0
✔ readr 1.3.1 ✔ forcats 0.5.0
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
Linking to GEOS 3.8.1, GDAL 3.0.4, PROJ 7.0.1

Attaching package: ‘janitor’

The following objects are masked from ‘package:stats’:

chisq.test, fisher.test

Attaching package: ‘glue’

The following object is masked from ‘package:dplyr’:

collapse

Linking to ImageMagick 6.9.10.23
Enabled features: fontconfig, freetype, fftw, lcms, pango, webp, x11
Disabled features: cairo, ghostscript, rsvg
Using ‘,’ as decimal and ‘.’ as grouping mark. Use read_delim() for more control.
Parsed with column specification:
cols(
dep = col_character(),
sexe = col_double(),
jour = col_date(format = «  »),
hosp = col_double(),
rea = col_double(),
rad = col_double(),
dc = col_double()
)
New names:
* «  -> …1
* «  -> …2
Warning message:
In st_point_on_surface.sf(.) :
st_point_on_surface assumes attributes are constant over geometries of x
reprojecting data…
reprojecting study zone…
generating reference grid…
computing kernel…
Note: Using an external vector in selections is ambiguous.
ℹ Use `all_of(field)` instead of `field` to silence this message.
ℹ See .
This message is displayed once per session.
Error: Assigned data `(1:nrow(dtCentroidesUniques)) – 1L` must be compatible with existing data.
✖ Existing data has 0 rows.
✖ Assigned data has 2 rows.
ℹ Only vectors of size 1 are recycled.
Backtrace:

1. └─covid_geo_pop %>% lissage(« deces », rayon, pixel, fr)
2. ├─base::withVisible(eval(quote(`_fseq`(`_lhs`)), env, env))
3. └─base::eval(quote(`_fseq`(`_lhs`)), env, env)
4. └─base::eval(quote(`_fseq`(`_lhs`)), env, env)
5. └─`_fseq`(`_lhs`)
6. └─magrittr::freduce(value, `_function_list`)
7. ├─base::withVisible(function_list[[k]](value))
8. └─function_list[[k]](value)
9. └─global::lissage(., « deces », rayon, pixel, fr)
10. └─`%>%`(…)
11. ├─base::withVisible(eval(quote(`_fseq`(`
Execution halted
« `

Un problème de version ?

Merci pour l’aide

oui, je crois

> require(« tidyverse »)
Loading required package: tidyverse
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
✔ ggplot2 3.3.0 ✔ purrr 0.3.4
✔ tibble 3.0.1 ✔ dplyr 0.8.5
✔ tidyr 1.1.0 ✔ stringr 1.4.0
✔ readr 1.3.1 ✔ forcats 0.5.0
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
> require(« httr »)
Loading required package: httr
> require(« fs »)
Loading required package: fs
> require(« sf »)
Loading required package: sf
Linking to GEOS 3.8.1, GDAL 3.0.4, PROJ 7.0.1
> require(« readxl »)
Loading required package: readxl
> require(« janitor »)
Loading required package: janitor

Attaching package: ‘janitor’

The following objects are masked from ‘package:stats’:

chisq.test, fisher.test

> require(« glue »)
Loading required package: glue

Attaching package: ‘glue’

The following object is masked from ‘package:dplyr’:

collapse

> require(« tmap »)
Loading required package: tmap
> require(« grid »)
Loading required package: grid
> require(« classInt »)
Loading required package: classInt
> require(« magick »)
Loading required package: magick
Linking to ImageMagick 6.9.10.23
Enabled features: fontconfig, freetype, fftw, lcms, pango, webp, x11
Disabled features: cairo, ghostscript, rsvg
> require(« btb »)
Loading required package: btb
> require(« raster »)
Loading required package: raster
Loading required package: sp

Attaching package: ‘raster’

The following object is masked from ‘package:glue’:

trim

The following object is masked from ‘package:janitor’:

crosstab

The following object is masked from ‘package:dplyr’:

select

The following object is masked from ‘package:tidyr’:

extract

> require(« fasterize »)
Loading required package: fasterize

Attaching package: ‘fasterize’

The following object is masked from ‘package:graphics’:

plot

> require(« plyr »)
Loading required package: plyr
——————————————————————————
You have loaded plyr after dplyr – this is likely to cause problems.
If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
library(plyr); library(dplyr)
——————————————————————————

Attaching package: ‘plyr’

The following objects are masked from ‘package:dplyr’:

arrange, count, desc, failwith, id, mutate, rename, summarise,
summarize

The following object is masked from ‘package:purrr’:

compact

>

En gros, qu’est-ce que fais le code pour obtenir le lissage à une date donnée ? On colorie d’abord les départements avec leur couleur puis on applique une fonction de lissage ?

On crée une grille pour la France, puis on génère les « centroïdes » des départements contenant le champ numérique à lisser et on utilise btb pour faire le lissage des points sur la grille.

Laisser un commentaire