This is a quick blog post to try and figure out which school boards are the “hottest zones” for covid.

In this blog post, I will do the following: - get the shapefiles for both the school boards (commission scolaires, or CS) and the smallest possible health districts (réseau local de santé, or RLS) - find the intersection between school boards and RLS using sf::st_intersection() - find the population of each intersection using cancensus::get_census() and tongfen::tongfen_estimate() - assign the number of cases in each health district (RLS) to each intersection proportionnaly to the population using tongfen::proportional_reaggregate() - sum the number of cases for each school board (commissions scolaires).

This post is massively dependent on Jens von Bergmann (@vb_jens) who created both the {cancensus} and the {tongfen} packages and also helped me getting started with tongfen. Thanks Jens!

First, download the number of covid cases at the RLS level from Claude Boucher’s google spreadsheet (and process data from my cronjob)

## claude data 
url_rls <- "https://docs.google.com/spreadsheets/d/17eJ05pg6fkzKlwwiWn0ztdMvzBX9pGx3q6w7h3_vJw4/"
gs4_deauth()
prout <- read_sheet(url_rls)

rls_population <- prout %>% select(RLS, Population) %>% filter(RLS != "Total") 

rls_claude  <- prout[1:which(prout$NoRLS == 1801), ] %>% 
  filter(!is.na(NoRLS)) %>% 
  select(-No,  -NoRLS,  -Population)%>% 
  gather(key=key, value=cumulative_cases, -RLS, -RSS) %>%
  mutate_at( vars(cumulative_cases) , as.numeric) %>% 
  filter(!is.na(cumulative_cases)) %>%
  #filter(NoRLS == "0112") %>%
  mutate(cumulative_cases = map_int(cumulative_cases, ~ as.integer(str_replace_all(.x, " ", "")))) %>% # remove spaces from numbers
  mutate(date_report = lubridate::ymd(key)) %>% 
  select(-key) %>%
  filter(!is.na(cumulative_cases))

# cronjo data


prepare_cronjob_data <- function(mypath = "~/cronjob/tableau-rls/"){
  pouet <- list.files(
    path =mypath,
    pattern = ".csv$", 
    full.names = TRUE)
  
  csvs <- purrr::map(pouet, read_csv)
  
  keep <- rep(TRUE, length(csvs))
  for (i in seq(from=2, to = length(csvs))){
    if ( identical(csvs[[i]],csvs[[i-1]])){ 
      keep[i]= FALSE
    }
  }
  dates <- lubridate::ymd(pouet %>% str_sub(start=-19, end = -12))
  
  
  datetimes <- lubridate::ymd_hms(
    paste0(
      pouet %>% str_sub(
        start=-19, end = -12),pouet %>% 
        str_sub(start=-10, end = -5) ))
  z1<- datetimes[keep]
  z2 <- dates[keep]
  z3 <- csvs[keep]
  
  dates_fixed <- tibble(download_datetime = z1) %>%
    mutate(download_hour = lubridate::hour(download_datetime),
           download_date = lubridate::date(download_datetime),
           report_date = if_else(download_hour >=8, download_date, download_date-1) # au milieu de la nuit tu es le rapport d'hier
    ) %>%
    arrange(desc(download_datetime)) %>%
    mutate(report_date_lag = lag(report_date),
           #si j'ai deux rapports différents dans la même journée, le premier est celui de la veille
           report_date_fixed = if_else(report_date == report_date_lag, report_date-1, report_date, report_date))  %>%
    arrange(download_datetime)
  
  gaa <- purrr::map2(z3, z1,
                     ~ .x %>% mutate(download_datetime = .y) )
  gaaa <- purrr::map2(gaa, dates_fixed$report_date_fixed,
                      ~ .x %>% mutate(date_report = .y) )
  
  rls_data <- gaaa %>%
    bind_rows() %>%
    rename(cumulative_cases = Cas) 
  
}

myrls <- prepare_cronjob_data() %>% 
  filter(!is.na(NoRLS), RLS != "Total") %>%
  mutate(
    cumulative_cases = if_else(cumulative_cases != "n.d.", cumulative_cases, NA_character_),
    cumulative_cases = as.numeric(str_replace_all(cumulative_cases, "\\s+", "")),
    Taux = if_else(Taux != "n.d.", Taux, NA_character_)
  ) %>%
  select(-No, -NoRLS, -Population, -Taux, -download_datetime) %>%
  filter(!is.na(cumulative_cases))


# combine both sources, fill in the blanks

both <- bind_rows(
  rls_claude %>% 
    mutate(source = "bouchecl") %>% 
    filter(!(date_report %in% unique(myrls$date_report))) ,
  myrls %>% 
    mutate(source = "cronjob") 
)


rls <- both %>% 
  group_by(RSS,RLS) %>% 
  arrange( RSS,RLS, desc(date_report)) %>% ## descending date to fix cumulative 
  mutate(fix_cummin = cummin(replace_na(cumulative_cases, Inf))) %>% 
  mutate(cumulative_cases = if_else(cumulative_cases > fix_cummin, fix_cummin, cumulative_cases, cumulative_cases)) %>% ## cumulative can't be bigger than next day.. if so reduce to next day level.
  arrange( RSS,RLS, date_report) %>% 
  mutate(slope = (lead(cumulative_cases)- cumulative_cases) / as.numeric(lead(date_report) - date_report))  %>%
  complete(date_report = seq.Date(min(date_report), max(date_report), by = "day")) %>%
  fill(slope, .direction= "down") %>%  
  mutate(cumulative_cases = if_else(!is.na(cumulative_cases), cumulative_cases, floor(first(cumulative_cases) + cumsum(slope) - slope))    ) %>%   
  select(-slope) %>%
  mutate(cases = cumulative_cases - lag(cumulative_cases)) %>% 
  ungroup() %>%
  filter(!is.na(cases))%>%
  mutate(shortname_rls = str_replace(str_extract(str_replace(RLS, "RLS de ","RLS " ),"RLS.+"),"RLS ", "")) %>%
  left_join(rls_population) %>%
  group_by(RSS, RLS) %>%
  mutate(cases_per_100k = cases * 1e5 / Population,
         cases_last_7_days = (cumulative_cases - lag(cumulative_cases,7)),
         cases_last_7_days_per_100k = cases_last_7_days * 1e5 / Population,
         RLS_code = str_extract(RLS, "^\\d+")
  )%>%
  ungroup()




prep_data <- function(data, group, type){
  
  type_column <- enquo(type)   ## this has to be !!
  type_name <- quo_name(type_column) ## its a string, dont !!
  mean_name = paste0("avg_", type_name, "_last7")   
  mean_column <- sym(mean_name)
  gaa <-   data %>% 
    group_by( {{ group }} ) %>%
    arrange(date_report) %>% 
    mutate(!!mean_name := ( !!type_column + lag(!!type_column, 1) + lag(!!type_column, 2) + lag(!!type_column, 3) + lag(!!type_column, 4) + lag(!!type_column, 5) + lag(!!type_column, 6)) / 7)  %>%
    ungroup() 
  
  gaa1 <- gaa %>%
    group_by( {{ group }}) %>%
    summarise(total = sum(!!type_column),
              worst7 = max(!!mean_column, na.rm = TRUE),
              last7  = max(!!mean_column * (date_report == max(date_report)), na.rm = TRUE),
              ratio = last7 / worst7,
              winning = factor( 
                case_when(ratio < 0.33 ~ "Winning",
                          ratio < 0.67 ~ "Nearly there",
                          TRUE ~ "Needs action"
                )
                , levels = c("Winning", "Nearly there", "Needs action"))
    ) %>%
    ungroup() %>%
    
    mutate(
      group = fct_reorder( {{ group }}, total, .desc =TRUE)
    )
  
  gaa %>% 
    left_join(gaa1) 
  
}

rls_cases <- prep_data(rls, shortname_rls, type = cases)
firstdate <- min(rls_cases$date_report)
lastdate  <- max(rls_cases$date_report)

We want to assign the covid cases from these “RLS” (health districts)

rls_shp <- read_sf(here::here("content/post/data/downloads/Territoires_RLS_2020.shp")) %>%
  st_transform(crs=4326) %>%
  rmapshaper::ms_simplify() %>%
  left_join(rls %>% filter(date_report ==max(date_report)) %>% select(date_report, RLS_code, cases_last_7_days_per_100k, cumulative_cases, cases_last_7_days, Population )) %>%
  mutate(dailycases_per_1M_avg_7_days = round(cases_last_7_days * 1e6 /7 / Population,1))

ggplot(data = rls_shp )+
  geom_sf(aes(fill=dailycases_per_1M_avg_7_days))+
  labs(title = paste0("There are ", nrow(rls_shp), " RLS"))

rls_shp %>% arrange(desc(dailycases_per_1M_avg_7_days)) %>%  
  mutate(rang = row_number())%>%
  select(rang, RLS_nom, dailycases_per_1M_avg_7_days, cases_last_7_days, Population) %>% 
  st_drop_geometry %>%
  kable(
    caption  = paste0("Nombre moyen de nouveaux cas de covid19 par million d'habitants \n pour les 7 derniers jours pour les RLS en date du  ", lastdate),
    col.names = c("Rang", "Nom RLS", "Nombre moyen de cas par million dans les 7 derniers jours", "Nombre de nouveaux cas dans les 7 derniers jours", "Population 2016")
  ) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE, position = "center") %>% 
  row_spec(0, background = "#555555", color = "#fff")
Table 1: Nombre moyen de nouveaux cas de covid19 par million d’habitants pour les 7 derniers jours pour les RLS en date du 2020-08-10
Rang Nom RLS Nombre moyen de cas par million dans les 7 derniers jours Nombre de nouveaux cas dans les 7 derniers jours Population 2016
1 RLS de Deux-Montagnes - Mirabel-Sud 50.4 44 124641
2 RLS du Haut-Saint-François 44.0 7 22749
3 RLS de la Haute-Gaspésie 40.5 3 10582
4 RLS des Faubourgs - Plateau-Mont-Royal - St-Louis-du-Parc 38.1 43 161317
5 RLS des Etchemins 34.8 4 16426
6 RLS de la Petite Patrie - Villeray 34.6 28 115485
7 RLS de Hochelaga - Mercier-Ouest - Rosemont 31.0 41 189095
8 RLS de Saint-Léonard - Saint-Michel 28.7 28 139162
9 RLS du Suroît 28.3 12 60523
10 RLS d’Ahuntsic - Montréal-Nord 27.3 33 172517
11 RLS de Verdun - Côte St-Paul - St-Henri - Pointe-St-Charles 25.8 30 166225
12 RLS de Dorval - Lachine - Lasalle 23.9 25 149522
13 RLS de Côte-Saint-Luc - NDG - Montréal-Ouest 23.9 22 131424
14 RLS de Lanaudière-Sud 23.1 48 296993
15 RLS du Nord de l’Île - Saint-Laurent 22.9 26 161894
16 RLS de Laval 22.9 71 442745
17 RLS Pierre-Boucher 22.9 42 262440
18 RLS de Côte-des-Neiges - Métro - Parc-Extension 22.7 39 245703
19 RLS de Richelieu-Yamaska 22.6 35 220957
20 RLS de la Rivière-du-Nord - Mirabel-Nord 21.7 27 177517
21 RLS de Champlain 21.5 34 225442
22 RLS de Lanaudière-Nord 20.5 32 223099
23 RLS de Vaudreuil-Soulanges 19.6 22 160002
24 RLS de Thérèse-De Blainville 18.4 21 163430
25 RLS de la Pommeraie 17.8 7 56125
26 RLS de Jardins-Roussillon 17.2 28 232373
27 RLS d’Argenteuil 17.1 4 33422
28 RLS de La Côte-de-Gaspé 15.2 2 18753
29 RLS du Haut-Richelieu - Rouville 14.6 20 195297
30 RLS de la Haute-Yamaska 13.5 10 105861
31 RLS de Grande-Rivière - Hull - Gatineau 13.5 25 264517
32 RLS des Pays-d’en-Haut 12.9 4 44362
33 RLS de Rivière-des-Prairies - Anjou - Montréal-Est 12.3 18 208717
34 RLS des Îles-de-la-Madeleine 11.3 1 12606
35 Nord-du-Québec 10.5 1 13558
36 RLS de Pierrefonds - Lac Saint-Louis 10.2 16 224596
37 RLS de la région de Thetford 10.1 3 42500
38 RLS d’Asbestos 10.0 1 14257
39 RLS de la Vallée-de-la-Lièvre et de la Petite-Nation 10.0 4 56920
40 RLS de la Baie-des-Chaleurs 8.9 2 32176
41 RLS Pierre-De Saurel 8.3 3 51371
42 RLS d’Antoine-Labelle 8.0 2 35580
43 RLS de Beauce 7.8 4 73173
44 Terres-Cries-de-la-Baie-James 7.8 1 18385
45 RLS de Coaticook 7.6 1 18741
46 RLS de Sherbrooke 7.6 9 170075
47 RLS de l’Abitibi-Ouest 7.0 1 20528
48 RLS de Kamouraska 6.9 1 20780
49 RLS du Granit 6.7 1 21351
50 RLS de Rouyn-Noranda 6.6 2 43180
51 RLS de Maskinongé 6.2 1 22961
52 RLS de Québec-Sud 5.8 13 321016
53 RLS de la Vallée-de-la-Batiscan 5.6 1 25562
54 RLS de Rivière-du-Loup 4.1 1 34800
55 RLS Alphonse-Desjardins 3.9 7 258118
56 RLS des Collines-de-l’Outaouais 3.7 1 38891
57 RLS de Québec-Nord 3.3 8 351429
58 RLS de la Vallée-de-l’Or 3.3 1 43551
59 RLS de Bécancour–Nicolet-Yamaska 3.2 1 44619
60 RLS des Laurentides 3.0 1 48016
61 RLS d’Arthabaska-et-de-L’Érable 2.9 2 98244
62 RLS de Lac-Saint-Jean-Est 2.7 1 52298
63 RLS de Memphrémagog 2.7 1 52514
64 RLS de Portneuf 2.6 1 54895
65 RLS de Rimouski 2.5 1 57630
66 RLS du Centre-de-la-Mauricie 2.2 1 64641
67 RLS de Jonquière 2.1 1 67417
68 RLS de Chicoutimi 1.8 1 79201
69 RLS de Drummond 1.3 1 107872
70 RLS de Trois-Rivières 1.0 1 143422
71 RLS de Témiscouata 0.0 0 19117
72 RLS de Matane 0.0 0 20726
73 RLS du Domaine-du-Roy 0.0 0 31078
74 RLS de Maria-Chapdelaine 0.0 0 25064
75 RLS de La Baie 0.0 0 22574
76 RLS de Charlevoix 0.0 0 28119
77 RLS de Val Saint-François 0.0 0 31490
78 RLS de la Vallée-de-la-Gatineau 0.0 0 20552
79 RLS de l’Abitibi 0.0 0 24729
80 RLS du Témiscaming 0.0 0 15628
81 RLS de Caniapiscau 0.0 0 3230
82 RLS de la Haute-Côte-Nord - Manicouagan 0.0 0 40494
83 RLS de Port-Cartier 0.0 0 7314
84 RLS de Sept-Îles 0.0 0 27344
85 RLS de la Minganie 0.0 0 6428
86 RLS de la Basse-Côte-Nord 0.0 0 4582
87 RLS de Kawawachikamach 0.0 0 641
88 RLS du Rocher-Percé 0.0 0 15881
89 RLS de Montmagny-L’Islet 0.0 0 40092
90 RLS du Haut-Saint-Laurent 0.0 0 24346
91 Nunavik 0.0 0 14260
92 RLS des Basques NA NA NA
93 RLS de La Mitis NA NA NA
94 RLS de La Matapédia NA NA NA
95 RLS du Haut-Saint-Maurice NA NA NA
96 RLS du Pontiac NA NA NA

To these “commission scolaires”.

CS_FRA_SDA <- read_sf(here::here("content/post/data/downloads/CS", "CS_FRA_SDA.shp"))%>%
  rmapshaper::ms_simplify()

ggplot(data = CS_FRA_SDA) +
  geom_sf()+
  labs(title = paste0("There are ", nrow(CS_FRA_SDA), " commissions scolaires"))

First, we create the smallest common denominators (intersections) between the two geographies using sf::st_intersection() and st_collection_extract()

intersections <- st_intersection(
  rls_shp %>% select(RLS_code),  ## RLS = région locale de service = health region,   RLS_code = health region id
  CS_FRA_SDA %>% select(CD_CS, NOM_CS)) %>%  ## CS = school board, CD_CS = school board id
  st_collection_extract(type="POLYGON") ## drop polylines

ggplot(data = intersections) +
  geom_sf()+
  labs(title = paste0("There are ", nrow(intersections), " intersections"))

Then, we will allocate the cases at the RLS to the lower geography level “intersections” proportionnally to the population of the intersections.

To do this, we must find the population of each intersection.
First, we download census population and shapefile using cancensus::get_census() Table below shows the population for a sample of the dissemination areas (DA)

census_data_da <- get_census(dataset='CA16', regions=list(PR="24"), vectors=c("pop2016"= "v_CA16_401"), level='DA',geo_format="sf") 

head(census_data_da %>% st_drop_geometry)
##   Shape Area Type Households CD_UID Dwellings   GeoUID Population CSD_UID
## 1   40.78502   DA        213   2401       254 24010018        465 2401042
## 2    7.91099   DA        230   2401       261 24010019        496 2401023
## 3   33.38888   DA        224   2401       248 24010020        523 2401023
## 4    2.94327   DA        212   2401       235 24010021        524 2401023
## 5    3.98348   DA        251   2401       285 24010022        530 2401023
## 6    5.70935   DA        251   2401       301 24010023        542 2401023
##   CMA_UID CT_UID              Region Name Area (sq km) pop2016
## 1    <NA>   <NA>               Grosse-Île     40.78502     465
## 2    <NA>   <NA> Les Îles-de-la-Madeleine      7.91099     496
## 3    <NA>   <NA> Les Îles-de-la-Madeleine     33.38888     523
## 4    <NA>   <NA> Les Îles-de-la-Madeleine      2.94327     524
## 5    <NA>   <NA> Les Îles-de-la-Madeleine      3.98348     530
## 6    <NA>   <NA> Les Îles-de-la-Madeleine      5.70935     542

Then we estimate the population of the “intersections” using tongfen::tongfen_estimate(). The estimates are “perfect” if the borders of the intersections line up with statcan’s “dissemination area”. When they dont, the population of a dissemination area is spread proportionally to the area covered by each intersection.

Table below shows the population for a sample of the intersections

intersections_populations <- tongfen_estimate(
  intersections,
  census_data_da,
  meta_for_additive_variables("census_data_da", "Population")
)

head(intersections_populations)
## Simple feature collection with 6 features and 4 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -72.09877 ymin: 45.00646 xmax: -71.39747 ymax: 45.7684
## geographic CRS: WGS 84
##     Population RLS_code  CD_CS                NOM_CS
## 1 1.501631e+01     0513 751000 CSS des Hauts-Cantons
## 2 3.087784e-01     0513 751000 CSS des Hauts-Cantons
## 3 2.717933e-01     0513 751000 CSS des Hauts-Cantons
## 4 1.667248e+04     0514 751000 CSS des Hauts-Cantons
## 5 1.741468e-01     0516 751000 CSS des Hauts-Cantons
## 6 1.966473e-02     0517 751000 CSS des Hauts-Cantons
##                         geometry
## 1 MULTIPOLYGON (((-71.90947 4...
## 2 MULTIPOLYGON (((-71.92083 4...
## 3 MULTIPOLYGON (((-71.91972 4...
## 4 MULTIPOLYGON (((-71.81586 4...
## 5 MULTIPOLYGON (((-71.71563 4...
## 6 MULTIPOLYGON (((-71.47856 4...

Assign the health districts covid cases to the lower geography level (intersections) using proportional_reaggregate(). This is population-weighted. Then, sum the number of cases and population in each intersection to get the number of cases and population at the school board (Commission scolaire) level.

The end result is this map and table:

cases_intersections <- tongfen::proportional_reaggregate(
  data = intersections_populations,
  parent_data = rls_shp ,
  geo_match = c("RLS_code" = "RLS_code"),
  categories = c("cases_last_7_days"),
  base = "Population"
)


cases_cs <- cases_intersections %>% 
  group_by(CD_CS, NOM_CS) %>% 
  summarise(
    cases_last_7_days = sum(cases_last_7_days, na.rm = TRUE),
    Population = sum(Population, na.rm = TRUE),
    
  ) %>%
  ungroup() %>%
  mutate(
    dailycases_per_1M_avg_7_days = round(cases_last_7_days * 1e6 /7 / Population,1),
    
    cases_last_7_days_per_100k = round(cases_last_7_days * 100000 / Population,1),
    cases_last_7_days = round(cases_last_7_days),
    Population = round(Population)
    
  )

ggplot(data = cases_cs) +
  geom_sf(aes(fill=dailycases_per_1M_avg_7_days))+ 
  #scale_fill_viridis_c() + 
  scale_fill_gradient(
    
    low = "white",
    high = "red",
    space = "Lab",
    na.value = "grey50",
    guide = "colourbar",
    aesthetics = "fill"
  ) + 
  labs(title = paste0("Nombre moyen de nouveaux cas de covid19 par million d'habitants pour les 7 derniers jours pour les commissions scolaires"),
       fill = "Cas par 1M habitants",
       caption = paste0("En date du " , format(lastdate, format="%B %d %Y"), ". Les données sont publiées au niveau des RLS (réseau local de service) et distribuées aux CSS proportionnellement à la population lorsque le RLS couvre plus d'une CSS")
  )+
  theme_bw() +
  theme(panel.grid.major = element_line(colour = "transparent"))

cases_cs %>% 
  ungroup %>% 
  st_drop_geometry() %>%
  arrange(-dailycases_per_1M_avg_7_days) %>%
  mutate(rang = row_number())%>%
  select(rang, NOM_CS, dailycases_per_1M_avg_7_days, cases_last_7_days, Population) %>% 
  kable(
    caption  = paste0("Nombre moyen de nouveaux cas de covid19 par million d'habitants \n pour les 7 derniers jours pour les commissions scolaires en date du  ", lastdate),
    col.names = c("Rang", "Nom Commission Scolaire", "Nombre moyen de cas par million dans les 7 derniers jours", "Nombre de nouveaux cas dans les 7 derniers jours", "Population 2016")) %>% 
  kable_styling(bootstrap_options = c("striped", "hover"), full_width = FALSE, position = "center") %>% 
  row_spec(0, background = "#555555", color = "#fff")
Table 2: Nombre moyen de nouveaux cas de covid19 par million d’habitants pour les 7 derniers jours pour les commissions scolaires en date du 2020-08-10
Rang Nom Commission Scolaire Nombre moyen de cas par million dans les 7 derniers jours Nombre de nouveaux cas dans les 7 derniers jours Population 2016
1 CSS de la Seigneurie-des-Mille-Îles 33.3 68 291498
2 CSS de Montréal 30.3 215 1014236
3 CSS des Chic-Chocs 24.3 5 29447
4 CSS des Affluents 24.2 45 267561
5 CSS de Laval 24.0 71 423003
6 CSS de Saint-Hyacinthe 23.3 17 102679
7 CSS de la Rivière-du-Nord 23.1 31 190460
8 CSS Marie-Victorin 23.1 56 347095
9 CSS des Patriotes 22.3 42 270108
10 CSS des Samares 21.6 31 207792
11 CSS de la Vallée-des-Tisserands 21.3 13 86777
12 CSS de la Pointe-de-l’Île 21.2 47 316733
13 CSS des Trois-Lacs 21.0 22 149319
14 CSS Marguerite-Bourgeoys 20.3 87 610999
15 CSS des Hauts-Cantons 19.6 7 51887
16 CSS des Grandes-Seigneuries 19.4 27 197898
17 CSS de la Baie-James 17.7 2 16177
18 CSS du Val-des-Cerfs 15.8 16 146620
19 CSS des Hautes-Rivières 15.5 17 155632
20 CSS des Portages-de-l’Outaouais 12.7 14 153452
21 CSS des Draveurs 12.6 12 140717
22 CSS des Îles 11.5 1 12475
23 CSS au Coeur-des-Vallées 10.5 4 53662
24 CSS de la Beauce-Etchemin 10.4 9 123825
25 CSS des Appalaches 9.9 3 45953
26 CSS des Laurentides 8.7 6 92009
27 CSS de la Région-de-Sherbrooke 8.5 11 179142
28 CSS de Sorel-Tracy 8.4 3 50995
29 CSS Pierre-Neveu 8.1 2 35249
30 CSS du Lac-Abitibi 7.0 1 20533
31 CSS de Rouyn-Noranda 6.7 2 41784
32 CSS des Découvreurs 6.0 6 138839
33 CSS René-Lévesque 5.9 2 48363
34 CSS de Kamouraska-Rivière-du-Loup 5.0 2 55735
35 CSS de la Capitale 5.0 10 278887
36 CSS des Navigateurs 4.0 5 183344
37 CSS des Premières-Seigneuries 3.4 5 230723
38 CSS de l’Or-et-des-Bois 3.4 1 43767
39 CSS des Sommets 3.3 2 84188
40 CSS de la Riveraine 3.3 1 43263
41 CSS des Bois-Francs 3.0 2 94346
42 CSS du Lac-Saint-Jean 2.7 1 52082
43 CSS de Portneuf 2.7 1 53072
44 CSS de l’Énergie 2.6 2 97132
45 CSS De La Jonquière 2.1 1 66717
46 CSS des Phares 1.9 1 74864
47 CSS du Chemin-du-Roy 1.9 2 168945
48 CSS de la Côte-du-Sud 1.8 1 70534
49 CSS des Rives-du-Saguenay 1.4 1 100969
50 CSS des Chênes 1.4 1 103398
51 CSS du Fleuve-et-des-Lacs 0.2 0 29130
52 CSS des Monts-et-Marées 0.0 0 39227
53 CSS du Pays-des-Bleuets 0.0 0 56735
54 CSS de Charlevoix 0.0 0 28363
55 CSS des Hauts-Bois-de-l’Outaouais 0.0 0 33871
56 CSS du Lac-Témiscamingue 0.0 0 15948
57 CSS Harricana 0.0 0 24652
58 CSS de l’Estuaire 0.0 0 41865
59 CSS du Fer 0.0 0 39323
60 CSS de la Moyenne-Côte-Nord 0.0 0 6424
cases_cs %>% 
  arrange(-dailycases_per_1M_avg_7_days)  %>%
  mutate(NOM_CS = as.factor(NOM_CS) )%>%
  ggplot(aes(x= fct_reorder(NOM_CS, dailycases_per_1M_avg_7_days), y= dailycases_per_1M_avg_7_days )) + 
  scale_y_continuous(expand = c(0, 0), limits = c(0, NA))+
  geom_col() +
  geom_hline(yintercept = 20, color = "red", size =2 ) +
  geom_text(x=1, y=1, hjust="left", label= "Selon le critère établi par le ministre Dubé, la courbe est aplatie sous 20 cas / million")+
  
  coord_flip() +
  
  cowplot::theme_half_open() +
  cowplot::background_grid() +
  colorblindr::scale_color_OkabeIto( ) +
  labs(
    title = "Nombre estimé de cas par 1 million d'habitant au cours des 7 derniers jours",
    subtitle = paste0("En date du " , format(lastdate, format="%B %d %Y")), 
    caption = paste0(". Les données sont publiées au niveau des RLS (réseau local de service) et distribuées aux CSS proportionnellement à la population lorsque le RLS couvre plus d'une CSS"),
    x = "Commission scolaire",
    y= "Nombre de cas par million"
  )