vaccinationFigs.Rmd

---
title: "Vaccination Figures"
author: "FD"
output: 
  html_document: 
      code_folding: hide
      toc: TRUE
      toc_float: TRUE
editor_options:
  chunk_output_type: console
---

Source <https://www.data.gouv.fr/en/datasets/donnees-relatives-aux-personnes-vaccinees-contre-la-covid-19-1/>

# Load data

```{r}
URL <- "https://www.data.gouv.fr/en/datasets/r/54dd5f8d-1e2e-4ccb-8fb8-eac68245befd"
dataFile <- paste0("data/FranceVacc.csv") # name file with today's date
download.file(URL, dataFile) # download file from repo
dat <- read.csv(dataFile, sep = ";", stringsAsFactors = FALSE)
head(dat)

plot(dat$n_cum_complet/dat$couv_complet)
```

Age classes in the dataset

```{r}
unique(dat$clage_vacsi)
```

Compute population sizes

```{r}
# Select later dates to avoid issues with low numbers
ii <- which(dat$jour > "2021-05-20")
# Population sizes
nn <- 100*(dat$n_cum_complet/dat$couv_complet)[ii]
# Remove Inf and NaN
nn[nn>10^8] <- NA
nn[is.nan(nn)] <- NA
# Check values
plot(nn)
# Compute population size by averaging computed values
popsizes <- aggregate(nn, by = list(age = dat[ii, "clage_vacsi"]), FUN = mean, na.rm = TRUE)
popsizes
```

# Plot

## Generic

```{r}
# Define colors 

col1D <- "#af8dc3"  # 1 dose
colFin <- "#7fbf7b" # Parcours fini
col0 <- "#d9d9d9"   # non vacciné
```

```{r}
# Source: https://gist.github.com/Jfortin1/72ef064469d1703c6b30
lighten <- function(color, factor = 1.4){
    col <- col2rgb(color)
    col <- col*factor
    col <- rgb(t(as.matrix(apply(col, 1, function(x) if (x > 255) 255 else x))), maxColorValue=255)
    col
}
```

Age class codes

> 0 : Tous âges
> 9 : 0-9
> 17 : 10-17
> 24 : 18-24
> 29 : 25-29
> 39 : 30-39
> 49 : 40-49
> 59 : 50-59
> 69 : 60-69
> 74 : 70-74
> 79 : 75-79
> 80 : 80 et +



```{r}
# I am not using this after all, but keeping it to avoid having to type it again if I ever need to use it
dic <- c("Tous âges", "0-9", "10-17", "18-24", "25-29", "30-39", "40-49", "50-59", "60-69", "70-74", "75-79", "80 et +")
names(dic) <- c(0, 9, 17, 24, 29, 39, 49, 59, 69, 74, 79, 80)
```


Age Categories

```{r}
# Define age categories as list: 
# - each list element is a big age categories, and 
# - it contains as a vector the corresponding age classes
cats <- list(c(24, 29, 39, 49), c(59), c(64, 69, 74), c(79, 80))
catNames <- c("18-49", "50-59", "60-74", "75+")

# Compute population sizes for each big age class
pop <- rep(0, length(cats)) # initialize
for (i in seq_along(pop)){
  pop[i] <- sum(popsizes[is.element(popsizes$age, unlist(cats[i])), "x"]) # Sum population sizes of the age cats making up the big age class
}
```


```{r}

# Function to compute the key values for each age class, 
# for the chosen subset of data
computeVals <- function(subdat, cats, pop){
  # subdat : subset of the data for the chosen date
  # cats : list of age classes, with ages categories as vectors
  # pop : vector of population sizes for the corresponding age classes, same order as elements of cats
  
  # Initialize output
  newdat <- as.data.frame(cbind(seq_along(cats)))
  newdat$couv_dose1 <- 0
  newdat$couv_complet <- 0
  
  for(i in newdat[, 1]){
    # Select subset of values for age class i
    ssdat <- subdat[is.element(subdat$clage_vacsi, unlist(cats[i])), ]
    # Compute proportions and numbers
    newdat[i, "n_cum_dose1"] <- sum(ssdat$n_cum_dose1)
    newdat[i, "n_cum_complet"] <- sum(ssdat$n_cum_complet)
    newdat[i, "couv_dose1"] <- sum(ssdat$n_cum_dose1)/pop[i]
    newdat[i, "couv_complet"] <- sum(ssdat$n_cum_complet)/pop[i]
  }
  # En cours is at least one minus finished
  newdat$n_cum_encours <- newdat$n_cum_dose1 - newdat$n_cum_complet
  newdat$couv_encours <- newdat$couv_dose1 - newdat$couv_complet

  # Output
  newdat
}
```

## Bars

```{r}
# Lines of the last days
iend <- which(dat$jour >= max(dat$jour, na.rm = TRUE))

subdat <- dat[iend, ]
subdat
```


```{r}
# Indices of the values to be plotted
ival <- iend


# Size parameters
recL <- 100 # Length of the rectangles
popFactor <- 10^(-6) # Scaling factor for population size
deltaZ <- 7 # Space between rectangles



# Function to plot the figure
plotProps <- function(lejour){
  # lejour: "YYYY-MM-DD" day to be plotted
  
  # Indices in the table corresponding to the chosen day
  ival <- which(dat$jour == lejour)
  
  # Subset of the data for the chosen day
  subdat <- dat[ival, ]
  
  # Compute key values
  newdat <- computeVals(subdat, cats, pop)
  
#--- 
  # Plotting

  # (I wrote x but it is actually y)  
  # xmin values of the rectangles
  xmin <- c(0, cumsum(pop)[-length(pop)])*popFactor + (seq_along(pop)-1)*deltaZ
  # xmax values of the rectangles
  xmax <- c(cumsum(pop))*popFactor + (seq_along(pop)-1)*deltaZ
  # Mid values
  xmid <- apply(cbind(xmin, xmax), 1, mean)
  
  # Initialize plot
  par(mar = c(3, 7, 4.5, 0.5))
  plot(c(rep(0, length(pop)), rep(recL, length(pop))), c(xmin, xmax), 
       type = "n", axes = FALSE, 
       xlab = "", ylab = "")
  
  # Plot base rectangles, color unvaccinated
  for (icat in seq_along(cats)){
    rect(0, xmin[icat], recL, xmax[icat], col = col0, border = gray(0, 0), lwd = 0)
  }
  
  # Add legend of the rectangles
  par(xpd = TRUE)
  adjx <- 1.
  xtext <- -2
  # Age categories
  text(rep(xtext, length(cats)), xmid, labels = paste0(catNames, " ans\n"), adj = c(adjx, 0.5))
  # Population size
  text(rep(xtext, length(cats)), xmid, labels = paste0("(", round(pop/10^6, 1), " millions hab)"), adj = c(adjx, 1), cex = 0.9)
  
  # Plot at least one dose rectangles
  for (icat in seq_along(cats)){
    rect(0, xmin[icat], newdat[icat, "couv_dose1"]*recL, xmax[icat], col = col1D, border = gray(0, 0), lwd = 0)
  }
  
  # Plot vaccination complete rectangles
  for (icat in seq_along(cats)){
    rect(0, xmin[icat], newdat[icat, "couv_complet"]*recL, xmax[icat], col = colFin, border = gray(0, 0), lwd = 0)
  }
  
  # Add legend
  par(xpd = TRUE)
  legend(x = "topleft", col = c(colFin, col1D, col0), pch = 15, legend = c("Vaccination finie", "En cours (1 dose)", "Non vacciné"), inset = c(0, -0.1), ncol = 3, cex = 0.9, bty = "n")
  
  # Axis
  par(mgp = c(1, 0.3, 0), tck = -0.02)
  axis(1, at = seq(0, 100, by = 10))
  
  mtext("@flodebarre | Données https://www.data.gouv.fr/en/datasets/r/54dd5f8d-1e2e-4ccb-8fb8-eac68245befd", side = 1, line = 1.8, cex = 0.65)
  
  # Add title: plotted day
  mtext(subdat[1, "jour"], side = 3, line = 2.5, cex = 1.3)
  
}

# Plot the figure for the last day
plotProps(max(dat$jour, na.rm = TRUE))
```

```{r}
# Vector of all vaccination days
vacDays <- sort(unique(dat$jour))

# Figure for each day
for(dd in vacDays){
  png(filename = paste0("pics/vacc_", dd, ".png"), width = 1100, height = 700, pointsize = 20)

  plotProps(dd)
  dev.off()
}
    
# Convert into gif
system("convert -quality 100% -delay 0.1 -loop 0 pics/vacc*.png pics/animationVaccination_bars.gif")
```

## Donuts

```{r}
drawArc <- function(r1, r2, a1, a2, npts = 100, deltax = 0, deltay = 0, ...){
  # r1 radius of the big circle
  # r2 radius of the inner circle
  # a1 angle start point, as proportion [0, 1]
  # a2 angle end point
  # npts nb of points in the arc
  # deltax x shift
  # deltay y shift
  
  dpi <- + pi/2
  polygon(x = deltax + c(r2 * cos(-2 * pi * a1 + dpi), r1 * cos(-2 * pi * seq(a1, a2, length.out = npts) + dpi), r2 * cos(-2 * pi * seq(a2, a1, length.out = npts) + dpi)), 
          y = deltay + c(r2 * sin(-2 * pi * a1 + dpi), r1 * sin(-2 * pi * seq(a1, a2, length.out = npts) + dpi), r2 * sin(-2 * pi * seq(a2, a1, length.out = npts) + dpi)), ...)
}

drawDonut <- function(vals, cols = c(col1D, colFin, col0), r1, fracRadius = 0.6, deltax = 0, deltay = 0){
  # vals : fraction values
  # cols : corresponding colors
  
  # Check categories add
  if(round(sum(vals), 3) != 1){
    stop("values need to sum to 1")
  }
  
  cumvals <- cumsum(vals)

  drawArc(r1, r2 = r1*fracRadius, a1 = 0, a2 = vals[1], deltax = deltax, deltay = deltay, 
          col = cols[1], border = gray(0,0))
  for(i in 2:length(vals)){
    drawArc(r1, r2 = r1*fracRadius, a1 = cumvals[i-1], a2 = cumvals[i], deltax = deltax, deltay = deltay, 
            col = cols[i], border = gray(0,0))
  }
}


plot(-150:150, asp = 1, type = "n")
drawDonut(c(0.3, 0.2, 0.5), r1 = 67, deltax = 0, deltay =0 )

drawArc(100, 50, 0., 0.3, col = "red", border = "blue", deltax = 100)
```

```{r}
# Function to plot the figure
plotDonuts <- function(lejour, donutSize = TRUE, innerSize = 0.5, cols = c(col1D, colFin, col0)){
  # lejour: "YYYY-MM-DD" day to be plotted
  # donutSize: bool whether size of the category is through size of outer cicle
  # innerSize: if donutSize = TRUE, relative size of the inner circle
  
  # Indices in the table corresponding to the chosen day
  ival <- which(dat$jour == lejour)
  
  # Subset of the data for the chosen day
  subdat <- dat[ival, ]
  
  # Compute key values
  newdat <- computeVals(subdat, cats, pop)

#--- 
  # Plotting
  
  maxpop <- max(pop)
  dx <- 100
  rayons <- sqrt(pop/pi)
  scaledRayons <- rayons / max(rayons)*dx
  scaledRayons
  
  par(mar = c(1, 0.5, 2.5, 0.5))
  xx <- cumsum(2*scaledRayons)-scaledRayons + (seq_along(cats) - 1)*0.1*dx

  plot(xx, rep(0, length(cats)), 
       xlim = c(min(xx) - dx, max(xx) + dx), 
       ylim = c(-dx, dx), 
       type = "n", xlab = "", ylab = "", asp = 1, axes = FALSE)
  
  for(i in seq_along(cats)){
    v <- c(newdat[i, "couv_complet"], newdat[i, "couv_encours"], 1 - (newdat[i, "couv_encours"] + newdat[i, "couv_complet"]))
    drawDonut(v, cols = c(colFin, col1D, col0), scaledRayons[i], innerSize, deltax = xx[i])
  }
  
  text(xx, y = 0, labels = catNames, adj = c(0.5, 0))
  text(xx, y = 0, labels = paste0(round(pop/10^6, 1), " M"), adj = c(0.5, 1.5), cex = 0.7)
  

  # Add legend
  par(xpd = TRUE)
  legend(x = "topleft", col = c(colFin, col1D, col0), pch = 15, legend = c("Vaccination finie", "En cours (1 dose)", "Non vacciné"), inset = c(0, -0.075), ncol = 3, cex = 0.9, bty = "n")
  
  
  mtext("Visualisation @nicolasberrod, animation @flodebarre | Données https://www.data.gouv.fr/en/datasets/r/54dd5f8d-1e2e-4ccb-8fb8-eac68245befd\n Code sur https://github.com/flodebarre/covid_vaccination", side = 1, line = 0, cex = 0.65)
  
  # Add title: plotted day
  mtext(subdat[1, "jour"], side = 3, line = 1, cex = 1.3)
}

plotDonuts("2021-06-02")
```



```{r}
# Vector of all vaccination days
vacDays <- sort(unique(dat$jour))

# Figure for each day
for(dd in vacDays){
  png(filename = paste0("pics/vacc-donut_", dd, ".png"), width = 900, height = 400, pointsize = 20)

  plotDonuts(dd)
  dev.off()
}
    
# Convert into gif
system("convert -quality 100% -delay 0.1 -loop 0 pics/vacc-donut_*.png pics/animationVaccination_donuts.gif")
```

## Time dynamics

```{r}
# Function to compute the key values for each age class, 
# for the chosen subset of data
computeTimeVals <- function(subdat, popAge){
  # subdat : data
  # popAge : population size of the age class

  days <- sort(unique(subdat$jour))
  
  # Initialize output
  newdat <- as.data.frame(cbind(seq_along(days)))
  newdat$jour <- days
  newdat$couv_dose1 <- 0
  newdat$couv_complet <- 0
  
  for(i in newdat[, 1]){
    # Select subset of values for day i
    ssdat <- subdat[is.element(subdat$jour, days[i]), ]
    ssdat
    # Compute proportions and numbers
    newdat[i, "n_cum_dose1"] <- sum(ssdat$n_cum_dose1)
    newdat[i, "n_cum_complet"] <- sum(ssdat$n_cum_complet)
    newdat[i, "couv_dose1"] <- sum(ssdat$n_cum_dose1)/popAge
    newdat[i, "couv_complet"] <- sum(ssdat$n_cum_complet)/popAge
  }
  # En cours is at least one minus finished
  newdat$n_cum_encours <- newdat$n_cum_dose1 - newdat$n_cum_complet
  newdat$couv_encours <- newdat$couv_dose1 - newdat$couv_complet

  # Output
  newdat
}
```


```{r}
# Function to plot the time dynamics for an age class
plotDyn <- function(icat, geog = "France"){
  newdat <- computeTimeVals(dat[is.element(dat$clage_vacsi, unlist(cats[icat])),], pop[icat])
  
  par(mar = c(4, 4, 3, 2))
  par(las = 1)
  par(xpd = FALSE)
  par(mgp = c(2, 0.5, 0), tck = -0.02)
  
  plot(as.Date(newdat$jour), newdat$couv_dose1, 
       ylim = c(0, 1), xlim = range(as.Date(dat$jour)), 
       type = "n", xaxs = "i", yaxs = "i", frame.plot = FALSE, xlab = "", ylab = "Proportion de la catégorie d'âge", axes = FALSE)
  
  opacityFill <- 0.5
  # Plot at least one dose
  polygon(as.Date(c(newdat$jour, rev(newdat$jour), newdat[1, "jour"])), 
          c(newdat$couv_dose1, rev(0* newdat$couv_complet), newdat[1, "couv_dose1"]), border = gray(0, 0), col = adjustcolor(col1D, opacityFill))
  
  # Plot vaccination completed
  polygon(as.Date(c(newdat$jour, rev(newdat$jour), newdat[1, "jour"])), 
          c(newdat$couv_complet, rep(0, nrow(newdat)), newdat[1, "couv_dose1"]), border = gray(0, 0), col = adjustcolor(colFin, opacityFill))
  # Grid lines
  for(i in seq(0, 1, by = 0.1)){
    abline(h = i, col = gray(0.85), lwd = 0.8)
  }
  lwdLines <- 2
  
  legend(x = "topleft", box.lwd = 0, lty = c(1, 1), col = c(col1D, colFin), legend = c("Au moins une dose", "Vaccination complète"), lwd = lwdLines)
  
  
  # Lines 
  lines(as.Date(newdat$jour), newdat$couv_dose1, col = col1D, lwd = lwdLines)
  lines(as.Date(newdat$jour), newdat$couv_complet, col = colFin, lwd = lwdLines)
  axis(2)
  axis(4)
  months <- c("2021-01-01", "2021-02-01", "2021-03-01", "2021-04-01", "2021-05-01", "2021-06-01")
  axis(1, at = as.Date(months), labels = substr(months, 0, 7))
  
    mtext("@flodebarre | Données https://www.data.gouv.fr/en/datasets/r/54dd5f8d-1e2e-4ccb-8fb8-eac68245befd\n Code: https://github.com/flodebarre/covid_vaccination", side = 1, line = 3, cex = 0.65)
    
  title(main = paste0(catNames[icat], " ans, ", geog))

}

plotDyn(1)
for(i in seq_along(cats)){
  
  png(filename = paste0("pics/vaccDyn_", i, ".png"), width = 800, height = 500, pointsize = 20)

 plotDyn(i) 
 
 dev.off()
}
```


## Squares


Age Categories

```{r}
# Define age categories as list: 
# - each list element is a big age categories, and 
# - it contains as a vector the corresponding age classes
cats <- list(c(24, 29), c(39), c(49), c(59), c(64, 69), c(74, 79), c(80))
catNames <- c("18-29", "30-39", "40-49", "50-59", "60-69", "70-79", "80+")

# Compute population sizes for each big age class
pop <- rep(0, length(cats)) # initialize
for (i in seq_along(pop)){
  pop[i] <- round(sum(popsizes[is.element(popsizes$age, unlist(cats[i])), "x"])) # Sum population sizes of the age cats making up the big age class
}
```


```{r}
nCats <- length(cats)

# Plot parameters
nbRow <- 4 # Nb of rows for each age class
scaleSquare <- 5*10^4 # Number of people represented by a square
nbSquares <- round(pop/scaleSquare) # Number of squares per category
nbCols <- ceiling(nbSquares / nbRow) # Corresponding number of columns

cote <- 0.9 # Space for square
deltaSquare <- 0.1 # Space between squares
# NB in the current version, only cote+deltaSquare is used
deltaCat <- 0.65 # Space between categories


plotSquares <- function(lejour){
  # lejour: day to be plotted
  
  # Initialize plot
  plot(0, type = "n", xlim = c(0, max(nbCols)+1), ylim = c(0,  nCats*(deltaCat + nbRow*(cote + deltaSquare))), axes = FALSE, xlab = "", ylab = "")
  
  # Subset of the data for the given day
  subdat <- dat[which(dat$jour == lejour), ]
  
  # Compute key quantities for each category
  newdat <- computeVals(subdat, cats, pop)
  
  # Add information about the corresp. number of squares
  newdat$nbSquares_dose1 <- round(newdat$n_cum_dose1/scaleSquare)
  newdat$nbSquares_complet <- round(newdat$n_cum_complet/scaleSquare)

  for(icat in seq_len(nCats)){
    # x positions
    xpos <- rep(1:max(nbCols), each = nbRow)
    # y positions
    ypos <- rep((icat - 1) * (nbRow * (cote + deltaSquare) + deltaCat) + ((1:nbRow)-1)*(cote + deltaSquare), max(nbCols))

    thepch <- 15 # Plot squares
    
    # Points of unvaccinated, until total size of the category
    points(xpos[1:nbSquares[icat]], ypos[1:nbSquares[icat]], pch = thepch, col = col0)
    
    # Condition necessary because indices start at 1, so otherwise would plot 1 square when still at 0
    if(newdat[icat, "nbSquares_dose1"] >= 1){
      # Points of at least one dose
      points(xpos[1:newdat[icat, "nbSquares_dose1"]], ypos[1:newdat[icat, "nbSquares_dose1"]], pch = thepch, col = col1D)
    }
  
    if(newdat[icat, "nbSquares_complet"] >= 1){
      # Points of fully vaccinated
      points(xpos[1:newdat[icat, "nbSquares_complet"]], ypos[1:newdat[icat, "nbSquares_complet"]], pch = thepch, col = colFin)
    }
    
    # Add legend
    par(xpd = TRUE)
    # Age categories
    text(0, mean(c(ypos[2], ypos[3])), adj = c(1, 0), labels = paste(catNames[icat], "ans"))
    # Corresponding population size
    text(0, mean(c(ypos[2], ypos[3])), adj = c(1, 1.5), labels = paste(round(pop[icat]/10^6, 1), "M"), cex = 0.75)
    par(xpd = FALSE)
  }
  
  mtext("@flodebarre \nDonnées https://www.data.gouv.fr/en/datasets/r/54dd5f8d-1e2e-4ccb-8fb8-eac68245befd\nCode https://github.com/flodebarre/covid_vaccination", side = 1, line = 2, cex = 0.65, adj = 0)
  
  title(main = lejour)
  
  par(xpd = TRUE)
  legend(x = "topleft", col = c(colFin, col1D, col0, 1), pch = c(rep(thepch, 3), 0), legend = c("Vaccination finie", "Vaccination en cours", "Non vacciné", paste(toString(scaleSquare), "personnes")), inset = c(-0.1, -0.05), bty = "n", cex = 0.8, pt.cex = 1, ncol = 4, text.width = c(5, 11, 12, 11))

}

for(i in dev.list())dev.off()
plotSquares(max(dat$jour))
plotSquares("2021-06-01")
```

```{r}
# Vector of all vaccination days
vacDays <- sort(unique(dat$jour))

# Figure for each day
for(dd in vacDays){
  png(filename = paste0("pics/vacc-squares_", dd, ".png"), width = 950, height = 750, pointsize = 22)
  plotSquares(dd)
  dev.off()
}
    
# Convert into gif
system("convert -quality 100% -delay 0.1 -loop 0 pics/vacc-squares_*.png pics/animationVaccination_squares.gif")
```

## Squares and time dynamics

```{r}
# Plot parameters
nbRow <- 8 # Nb of rows for each age class
scaleSquare <- 5*10^4 # Number of people represented by a square
nbSquares <- round(pop/scaleSquare) # Number of squares per category
nbCols <- ceiling(nbSquares / nbRow) # Corresponding number of columns

cote <- 0.9 # Space for square
deltaSquare <- 0.1 # Space between squares
# NB in the current version, only cote+deltaSquare is used
deltaCat <- 1.5 # Space between categories

distDyn <- 2
widthDyn <- 25

plotSquaresDyn <- function(lejour){
  # lejour: day to be plotted

  # Initialize plot
  plot(0, type = "n", xlim = c(0, max(nbCols)+1 + distDyn + widthDyn), ylim = c(0,  nCats*(deltaCat + nbRow*(cote + deltaSquare))), axes = FALSE, xlab = "", ylab = "")
  
  # Subset of the data for the given day
  subdat <- dat[which(dat$jour == lejour), ]
  
  # Compute key quantities for each category
  newdat <- computeVals(subdat, cats, pop)
  
  # Add information about the corresp. number of squares
  newdat$nbSquares_dose1 <- round(newdat$n_cum_dose1/scaleSquare)
  newdat$nbSquares_complet <- round(newdat$n_cum_complet/scaleSquare)

  for(icat in seq_len(nCats)){
    # x positions
    xpos <- rep(1:max(nbCols), each = nbRow)
    # y positions
    ypos <- rep((icat - 1) * (nbRow * (cote + deltaSquare) + deltaCat) + ((1:nbRow)-1)*(cote + deltaSquare), max(nbCols))
    
    yTmin <- min(ypos[1:max(nbCols)])
    yTmax <- max(ypos[1:max(nbCols)])
    xTmin <- max(nbCols) + distDyn
    xTmax <- xTmin + widthDyn

    # Frame of time plot
    rect(xTmin, yTmin, xTmax, yTmax, border = gray(0.5))

    par(xpd = TRUE)    
    thecex <- 0.6
    if(icat == 1){
      text(mean(c(xTmin, xTmax)), y = yTmin, adj = c(0.5, 1.5), labels = "temps", cex = thecex)
      
    }
    text(rep(xTmax, 2), c(yTmin,yTmax), labels = c(" 0%", " 100%"), adj = c(0, 0.5), cex = thecex)

    par(xpd = FALSE)
    thepch <- 15 # Plot squares
    
    # Points of unvaccinated, until total size of the category
    points(xpos[1:nbSquares[icat]], ypos[1:nbSquares[icat]], pch = thepch, col = col0)
    
    # Condition necessary because indices start at 1, so otherwise would plot 1 square when still at 0
    if(newdat[icat, "nbSquares_dose1"] >= 1){
      # Points of at least one dose
      points(xpos[1:newdat[icat, "nbSquares_dose1"]], ypos[1:newdat[icat, "nbSquares_dose1"]], pch = thepch, col = col1D)
    }
  
    if(newdat[icat, "nbSquares_complet"] >= 1){
      # Points of fully vaccinated
      points(xpos[1:newdat[icat, "nbSquares_complet"]], ypos[1:newdat[icat, "nbSquares_complet"]], pch = thepch, col = colFin)
    }
    
    # Add legend
    par(xpd = TRUE)
    # Age categories
    text(0, mean(ypos), adj = c(1, 0), labels = paste(catNames[icat], "ans"))
    # Corresponding population size
    text(0, mean(ypos), adj = c(1, 1.5), labels = paste(round(pop[icat]/10^6, 1), "M"), cex = 0.75)
    par(xpd = FALSE)
    
    
  # Add time dynamics
  # Dataset
  newdatT <- computeTimeVals(dat[is.element(dat$clage_vacsi, unlist(cats[icat])),], pop[icat])
  
  vX <- as.numeric(as.Date(newdatT$jour))
  xT.all <- xTmin + (xTmax - xTmin) * (vX - min(vX))/(max(vX) - min(vX))
  
  yT_complet <- yTmin + (yTmax - yTmin) * newdatT$couv_complet
  yT_1D <- yTmin + (yTmax - yTmin) * newdatT$couv_dose1

  
  i.jour <- which(sort(unique(dat$jour)) == lejour)
  opacityFill <- 0.5
  lwdT <- 2
  
  # Fill areas
  polygon(c(xT.all[1:i.jour], rev(xT.all[1:i.jour]), xT.all[1]),
          c(rep(yTmax, i.jour), rev(yT_1D[1:i.jour]), 1), border = gray(0, 0), col = adjustcolor(col0, opacityFill))

  polygon(c(xT.all[1:i.jour], rev(xT.all[1:i.jour]), xT.all[1]),
          c(yT_1D[1:i.jour], rev(yT_complet[1:i.jour]), 1), border = gray(0, 0), col = adjustcolor(col1D, opacityFill))
  
  polygon(c(xT.all[1:i.jour], rev(xT.all[1:i.jour]), xT.all[1]),
          c(yT_complet[1:i.jour], rep(yTmin, i.jour), 1), border = gray(0, 0), col = adjustcolor(colFin, opacityFill))
  
  lines(xT.all[1:i.jour], yT_1D[1:i.jour], col = col1D, lwd = lwdT)
  lines(xT.all[1:i.jour], yT_complet[1:i.jour], col = colFin, lwd = lwdT)
  
  # Plot at least one dose
  polygon(as.Date(c(newdatT$jour, rev(newdatT$jour), newdatT[1, "jour"])), 
          c(newdatT$couv_dose1, rev(0* newdatT$couv_complet), newdatT[1, "couv_dose1"]), border = gray(0, 0), col = adjustcolor(col1D, opacityFill))
  
  # Plot vaccination completed
  polygon(as.Date(c(newdatT$jour, rev(newdatT$jour), newdatT[1, "jour"])), 
          c(newdatT$couv_complet, rep(0, nrow(newdatT)), newdatT[1, "couv_dose1"]), border = gray(0, 0), col = adjustcolor(colFin, opacityFill))
  
  }
  
  mtext("@flodebarre \nDonnées https://www.data.gouv.fr/en/datasets/r/54dd5f8d-1e2e-4ccb-8fb8-eac68245befd\nCode https://github.com/flodebarre/covid_vaccination", side = 1, line = 2, cex = 0.65, adj = 0)
  
  title(main = lejour)
  
  par(xpd = TRUE)
  legend(x = "topleft", col = c(1, colFin, col1D, col0), pch = c(0, rep(thepch, 3)), legend = c(paste(toString(scaleSquare), "personnes"), "Vaccination finie", "Vaccination en cours", "Non vacciné"), inset = c(-0., -0.05), bty = "n", cex = 0.8, pt.cex = 1, ncol = 4, text.width = c(5, 11, 12, 11))
  


}

for(i in dev.list())dev.off()

png("test.png", width = 650, height = 650, pointsize = 12)

#  png(filename = paste0("pics/vaccSqDyn_", i, ".png"), width = 800, height = 500, pointsize = 20)

plotSquaresDyn(max(dat$jour))
dev.off()
system("open test.png")
```


```{r}
# Vector of all vaccination days
vacDays <- sort(unique(dat$jour))

# Figure for each day
for(dd in vacDays){
  png(filename = paste0("pics/vacc-SquaresDyn_", dd, ".png"), 
      width = 650, height = 650, pointsize = 12)
  plotSquaresDyn(dd)
  dev.off()
}
    
# Convert into gif
system("convert -quality 100% -delay 0.1 -loop 0 pics/vacc-SquaresDyn_*.png pics/animationVaccination_squaresDyn.gif")

system(paste0("open pics/vacc-SquaresDyn_", dd, ".png"))

system(paste0("open pics/vacc-SquaresDyn_", dd, ".png"))
```