Geospatial Visualization

Both the coronavirus and covid19_vaccine datasets provide country-level information on the covid19 cases and vaccination progress, respectively. One common approach for communicating country-level data is with the use of choropleth maps. The focus of this vignette is to demonstrate how to merge the dataset with the geometric metadata and plot it. We will use the following packages:

• rnaturalearth - Provides geo-spatial metadata from the Natural Earth dataset. More details available here

• sf - A package that provides simple features access for R

• mapview - A wrapper for the leaflet library

• tmap - A package for creating a thematic maps

• ggplot2 - Is a system for declaratively creating graphics

• viridis - A package that provide a series of color maps

library(rnaturalearth)
library(sf)
library(mapview)
library(tmap)
library(ggplot2)
library(viridis)

Note: This vignette is not available on the CRAN version (due to size limitation). Therefore, as the packages above are not on the dependencies list of the coronavirus package, you may need to install them before.

Additional setting: following changes in the default options of the sf package from version 1.0-1 by default option is to use s2 spherical geometry as default when coordinates are ellipsoidal. That cause some issues with the tmap package, therefore we will set this functionality as FALSE:

sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off

More details available on this issue and follow-up issue.

Let’s get start by loading the data:

library(coronavirus)

data("covid19_vaccine")

head(covid19_vaccine)
#>   country_region       date doses_admin people_partially_vaccinated
#> 1         Canada 2020-12-14           5                           0
#> 2          World 2020-12-14           5                           0
#> 3         Canada 2020-12-15         723                           0
#> 4          China 2020-12-15     1500000                           0
#> 5         Russia 2020-12-15       28500                       28500
#> 6          World 2020-12-15     1529223                       28500
#>   people_fully_vaccinated report_date_string uid province_state iso2 iso3 code3
#> 1                       0         2020-12-14 124           <NA>   CA  CAN   124
#> 2                       0         2020-12-14  NA           <NA> <NA> <NA>    NA
#> 3                       0         2020-12-15 124           <NA>   CA  CAN   124
#> 4                       0         2020-12-15 156           <NA>   CN  CHN   156
#> 5                       0         2020-12-15 643           <NA>   RU  RUS   643
#> 6                       0         2020-12-15  NA           <NA> <NA> <NA>    NA
#>   fips      lat     long combined_key population continent_name continent_code
#> 1 <NA> 60.00000 -95.0000       Canada   37855702  North America             NA
#> 2 <NA>       NA       NA         <NA>         NA           <NA>           <NA>
#> 3 <NA> 60.00000 -95.0000       Canada   37855702  North America             NA
#> 4 <NA> 35.86170 104.1954        China 1404676330           Asia             AS
#> 5 <NA> 61.52401 105.3188       Russia  145934460         Europe             EU
#> 6 <NA>       NA       NA         <NA>         NA           <NA>           <NA>

We will use the ne_countries function from the rnaturalearth package to pull the country geometric data:

library(dplyr)

map <- ne_countries(returnclass = "sf") %>%
  dplyr::select(name, iso2 = iso_a2, iso3 = iso_a3, geometry)

head(map)
#> Simple feature collection with 6 features and 3 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -73.41544 ymin: -55.25 xmax: 75.15803 ymax: 42.68825
#> CRS:           +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
#>                   name iso2 iso3                       geometry
#> 0          Afghanistan   AF  AFG MULTIPOLYGON (((61.21082 35...
#> 1               Angola   AO  AGO MULTIPOLYGON (((16.32653 -5...
#> 2              Albania   AL  ALB MULTIPOLYGON (((20.59025 41...
#> 3 United Arab Emirates   AE  ARE MULTIPOLYGON (((51.57952 24...
#> 4            Argentina   AR  ARG MULTIPOLYGON (((-65.5 -55.2...
#> 5              Armenia   AM  ARM MULTIPOLYGON (((43.58275 41...

df <- map %>% left_join(
   covid19_vaccine %>%
    filter(date == max(date),
           is.na(province_state)) %>%
    mutate(perc = round(100 * people_fully_vaccinated / population, 2)) %>%
    select(country_region, iso2, iso3, people_fully_vaccinated, perc, continent_name),
    by = c("iso2", "iso3")
)

class(df)
#> [1] "sf"         "data.frame"

After we merged the country data with the corresponding geometry data it is straightforward to plot the data as sf object.

Choropleth maps with the mapview package

The mapview package, a wrapper for the leaflet library, enables to plot sf objects seamlessly. Let’s start by plotting the percentage of the population that fully vaccinated by country using the perc variable:

df  %>%
  mapview::mapview(zcol = "perc")

By default, the function is using continues color scale for the objects (in this case countries) color. We can modify it and set color buckets by using the at argument. Also, we can define the legend title with the use of the layer.name argument:

df  %>%
  mapview::mapview(zcol = "perc", 
                   at = seq(0, max(df$perc, na.rm = TRUE), 10), 
                   legend = TRUE,
                   layer.name = "Fully Vaccinated %")

Some of the missing values in the plot are un-populated areas such as Antarctica or a territory that is count under different state such as Greenland. We can remove those and re-plot the map:

df1 <- df %>% 
  filter(!name %in% c("Greenland", "Antarctica"))

df1  %>%
  mapview::mapview(zcol = "perc", 
                   at = seq(0, max(df1$perc, na.rm = TRUE), 10), 
                   legend = TRUE,
                   layer.name = "Fully Vaccinated %")