Drawing maps with the package covid19br

Fábio N. Demarqui

When dealing with epidemiological data, we are often interested in visualizing the data through different types of maps. Drawing maps in R nowadays is a simple task provided the necessary geometry is available along with the data.

The function covid19br::add_geo() can be used to add the geometry to the downloaded data as follows:


library(covid19br)
library(tidyverse)

# downloading data at state level:
cities <- downloadCovid19("cities") 

# adding the geometry to the data:
cities_geo <- cities %>%
  filter(date == max(date)) %>%
  add_geo()

# looking at the data:
glimpse(cities_geo)
#> Rows: 5,570
#> Columns: 28
#> $ region            <chr> "North", "North", "North", "North", "North", "North"…
#> $ state             <chr> "RO", "RO", "RO", "RO", "RO", "RO", "RO", "RO", "RO"…
#> $ city              <chr> "Alta Floresta D'Oeste", "Ariquemes", "Cabixi", "Cac…
#> $ state_code        <dbl> 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, 11, …
#> $ city_code         <dbl> 110001, 110002, 110003, 110004, 110005, 110006, 1100…
#> $ healthRegion_code <int> 11005, 11001, 11006, 11002, 11006, 11006, 11006, 110…
#> $ healthRegion      <chr> "ZONA DA MATA", "VALE DO JAMARI", "CONE SUL", "CAFE"…
#> $ date              <date> 2021-10-16, 2021-10-16, 2021-10-16, 2021-10-16, 202…
#> $ epi_week          <int> 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, …
#> $ pop               <dbl> 22945, 107863, 5312, 85359, 16323, 15882, 7391, 1833…
#> $ accumCases        <int> 4183, 22769, 927, 15011, 2594, 2386, 861, 2001, 3777…
#> $ newCases          <int> 8, 0, 0, 2, 14, 0, 3, 6, 4, 0, 11, 4, 0, 0, 0, 0, 36…
#> $ accumDeaths       <int> 66, 500, 21, 306, 64, 45, 22, 38, 76, 224, 184, 613,…
#> $ newDeaths         <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0…
#> $ newRecovered      <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
#> $ newFollowup       <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, …
#> $ metrop_area       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0…
#> $ capital           <lgl> FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FAL…
#> $ DHI               <dbl> 0.641, 0.702, 0.650, 0.718, 0.692, 0.685, 0.613, 0.6…
#> $ EDHI              <dbl> 0.526, 0.600, 0.559, 0.620, 0.602, 0.584, 0.473, 0.4…
#> $ LDHI              <dbl> 0.763, 0.806, 0.757, 0.821, 0.799, 0.814, 0.774, 0.7…
#> $ IDHI              <dbl> 0.657, 0.716, 0.650, 0.727, 0.688, 0.676, 0.630, 0.6…
#> $ region_code       <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
#> $ mesoregion_code   <int> 1102, 1102, 1102, 1102, 1102, 1102, 1102, 1101, 1102…
#> $ microregion_code  <int> 11006, 11003, 11008, 11006, 11008, 11008, 11008, 110…
#> $ area              [km^2] 7112.842 [km^2], 4439.795 [km^2], 1319.708 [km^2], …
#> $ demoDens          [1/km^2] 3.225855 [1/km^2], 24.294590 [1/km^2], 4.025134 […
#> $ geometry          <MULTIPOLYGON [°]> MULTIPOLYGON (((-62.05044 -..., MULTIPO…

The map of the accumulated number of deaths by city can be easily drawn using the package ggplot2 as illustrated below:

ggplot(cities_geo) +
  geom_sf(aes(fill = accumDeaths)) 

Suppose now that we want to draw a map with the incidence of the COVID-19 in the cities belonging to Minas Gerais (MG) state. This can be done as follows:


mg <- cities_geo %>%
  filter(state == "MG") %>%
  add_epi_rates()

ggplot(mg) +
  geom_sf(aes(fill = incidence)) 

Next, we will show how to draw interactive maps using the package leaflet. As an illustration, we will consider the lethality by states:

library(leaflet)

# downloading data at state level:
states <- downloadCovid19("states") 

# adding the geometry to the data:
states_geo <- states %>%
  filter(date == max(date)) %>%
  add_geo() %>%
  add_epi_rates()

# looking at the data:
glimpse(states_geo)
#> Rows: 27
#> Columns: 23
#> $ region       <chr> "North", "North", "North", "North", "North", "North", "No…
#> $ state        <chr> "RO", "AC", "AM", "RR", "PA", "AP", "TO", "MA", "PI", "CE…
#> $ date         <date> 2021-10-16, 2021-10-16, 2021-10-16, 2021-10-16, 2021-10-…
#> $ epi_week     <int> 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 41, 4…
#> $ newCases     <int> 121, 4, 44, 5, 406, 30, 156, 96, 186, 115, 373, 234, 408,…
#> $ accumCases   <int> 267500, 88005, 427147, 126305, 595130, 123256, 226727, 35…
#> $ newDeaths    <int> 1, 0, 2, 0, 4, 1, 3, 1, 6, 29, 2, 3, 8, 3, 0, 10, 74, 7, …
#> $ accumDeaths  <int> 6551, 1842, 13756, 2016, 16706, 1988, 3831, 10211, 7064, …
#> $ newRecovered <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
#> $ newFollowup  <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
#> $ pop          <dbl> 1777225, 881935, 4144597, 605761, 8602865, 845731, 157286…
#> $ state_code   <dbl> 11, 12, 13, 14, 15, 16, 17, 21, 22, 23, 24, 25, 26, 27, 2…
#> $ DHI          <dbl> 33.490, 12.894, 35.037, 9.153, NA, 10.285, 88.950, 125.03…
#> $ EDHI         <dbl> 26.854, 9.949, 27.090, 7.488, NA, 8.799, 75.863, 106.031,…
#> $ LDHI         <dbl> 41.019, 16.865, 47.464, 11.971, NA, 12.543, 109.778, 160.…
#> $ IDHI         <dbl> 34.225, 12.880, 33.797, 8.668, NA, 9.902, 84.859, 115.371…
#> $ region_code  <dbl> 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, …
#> $ area         [km^2] 238719.84 [km^2], 164806.00 [km^2], 1565986.97 [km^2], 2…
#> $ demoDens     [1/km^2] 7.444815 [1/km^2], 5.351353 [1/km^2], 2.646636 [1/km^2…
#> $ geometry     <GEOMETRY [°]> POLYGON ((-60.70931 -13.693..., POLYGON ((-68.38…
#> $ incidence    <dbl> 15051.555, 9978.627, 10306.117, 20850.633, 6917.812, 1457…
#> $ lethality    <dbl> 2.45, 2.09, 3.22, 1.60, 2.81, 1.61, 1.69, 2.85, 2.19, 2.5…
#> $ mortality    <dbl> 368.6084, 208.8589, 331.9020, 332.8045, 194.1911, 235.062…

reds = colorBin("Reds", domain = states_geo$lethality, bins = 5)
mymap <- states_geo %>%
  leaflet() %>%
  addPolygons(fillOpacity = 1, 
              weight = 1,
              smoothFactor = 0.2,
              color = "gray",
              fillColor = ~ reds(lethality),
              popup = paste0(states_geo$state, ":  ",  states_geo$lethality, 2)
  ) %>%
  addLegend(position = "bottomright", 
            pal = reds, values = ~states_geo$lethality, 
            title = "lethality")
mymap