Interactive Maps in R

This post shows you how to create interactive maps in R using the highcharter package.

Why Interactive Maps?

Static maps are an effective visual tool that communicate geographic data in an interpretive way that is generally lost if that data is only viewed in a spreadsheet.

Adding interactivity to a map further improves data interpretability by allowing users to:

  • Explore the data by zooming in on areas of interest
  • Choose what data values to be displayed or excluded
  • Hover over an area of interest and get additional info/the exact value that’s being displayed

Interactive Choropleth Map

The first interactive map I will show you how to create is a population density map by county for the state of Texas:


Prerequisites

Before we get started, you will need to load the following packages:

library(data.table)
library(highcharter)
library(dplyr)
library(tidyr)
library(RColorBrewer)

Data Prep

The dataset we’ll be using is from the USDA ERS. The data is available in both xlsx and csv format. I downloaded the data in csv format and loaded the People.csv file which contains the population estimates:

# Load USDA ERS dataset
county_df <- fread("C:/Users/Kyle/Downloads/People.csv") %>% 
  filter(State == 'TX')

The dataset contains FIPS codes for each county. In order to map this data, we will need to join it to the dataset containing the geographic information for each county from the highcharter package.

# Load Texas county map
tx_counties <- get_data_from_map(download_map_data("countries/us/us-tx-all"))

glimpse(tx_counties)
## Observations: 254
## Variables: 7
## $ `hc-group`    <chr> "admin2", "admin2", "admin2", "admin2", "admin2", "ad...
## $ `hc-middle-x` <dbl> 0.50, 0.50, 0.50, 0.50, 0.50, 0.52, 0.36, 0.50, 0.50,...
## $ `hc-middle-y` <dbl> 0.50, 0.50, 0.50, 0.50, 0.77, 0.50, 0.36, 0.50, 0.50,...
## $ `hc-key`      <chr> "us-tx-179", "us-tx-393", "us-tx-311", "us-tx-131", "...
## $ `hc-a2`       <chr> "GR", "RO", "MC", "DU", "LO", "LE", "HO", "LA", "BA",...
## $ fips          <chr> "48179", "48393", "48311", "48131", "48297", "48289",...
## $ name          <chr> "Gray", "Roberts", "McMullen", "Duval", "Live Oak", "...

Before joining the datasets, we need to calculate the population density for 2018. We can do this by using mutate from the dplyr package to create a column called ‘density’ and divide the ‘TotalPopEst2018’ column by the ‘LandAreaSQMiles2010’ column as shown below:

# Calculate population density
density_df <- county_df %>% 
  select(c(fips=FIPS,County,TotalPopEst2018,LandAreaSQMiles2010)) %>%
  mutate(density = round(TotalPopEst2018/LandAreaSQMiles2010,0))

head(density_df)
##    fips   County TotalPopEst2018 LandAreaSQMiles2010 density
## 1 48000    Texas        28701845           261231.71     110
## 2 48001 Anderson           58057             1062.60      55
## 3 48003  Andrews           18128             1500.71      12
## 4 48005 Angelina           87092              797.78     109
## 5 48007  Aransas           23792              252.07      94
## 6 48009   Archer            8786              903.11      10

Creating the Map

Using the hcmap function from highcharter, we can create a basic interactive map like this one:

# Create interactive map of Texas counties
density_map <- hcmap(map = "countries/us/us-tx-all", 
      data = density_df, 
      value = "density", 
      joinBy = c("fips"))  %>%
  hc_mapNavigation(enabled = TRUE)


Customizing the Map

In the map above, the counties with the highest population densities clearly stand out. However, it is difficult to distinguish the differences between the counties with lower population densities.

In order to correct this, we can assign color breaks to the data by utilizing the hc_colorAxis function and assigning a color palette from the RColorBrewer package:

# Add color classes and legend to map
density_map <- hcmap(map = "countries/us/us-tx-all", 
      data = density_df, 
      value = "density", 
      joinBy = c("fips"),
      borderWidth = 0.3)  %>%
  hc_mapNavigation(enabled = TRUE) %>% 
  hc_legend(layout = "vertical", 
            align = "right",
            valueDecimals = 0) %>% 
  hc_colorAxis(dataClasses = color_classes(breaks = c(0,10,25,50,100,250,500,1000,2500,max(density_df$density)),
              colors = brewer.pal(name="YlOrRd",n=9)))

The next thing we will need to do is modify what is displayed when you hover (or click if you are viewing on a mobile device) on a particular county.

By creating a JavaScript function within the hc_tooltip option, we can more clearly display the name and population density for each county.

# Add custom tooltip to map
density_map <- density_map %>%
    hc_tooltip(formatter = JS("function() {
  return ('<br><b>County:</b> ' + this.point.County +
          '<br><b>Population Density:</b> ' + this.point.density + ' people per sq mi'
  )}"))

Now, the only thing left is to add a title and source to the map which we can do with the hc_title and hc_credits options:

# Add title and credits to map
density_map <- density_map %>% 
  hc_title(text = "Population Density by County, 2018") %>%
  hc_credits(enabled = TRUE,
             text = "Author: Kyle Cuilla, Data: USDA ERS",
             href = "https://www.ers.usda.gov/data-products/atlas-of-rural-and-small-town-america/download-the-data/")

And here is our final result!


Animated Choropleth Map

So now we have our map that displays the population density by county in 2018.

Let’s say that we want to see how the population density has changed over time. How would we go about doing this?

Well, we could create nine separate maps (one for each year from 2010 to 2018), but this would take up a lot of space and because the maps would each be separate, and because each map would be separate, it may be difficult to detect subtle difference between each year.

To solve these issues, we can create an animated map instead:


Data Prep

The county_df dataset we’ve been using contains estimated populations for each year.

We can calculate the population densities for each of these years by creating a function called ‘pop_density’ and applying it to each population estimate.

We can then use the gather function from the tidyr package to put all of the population densities into a single column called ‘density’ and all of the years into a single column called ‘years’.

# Calculate population density for each year in dataset
pop_density <- function(x){
  round(x/county_df$LandAreaSQMiles2010,0)
}

density_df_by_year <- county_df %>% 
  select(c(FIPS,State,County,
           '2010'=TotalPopEst2010,
           '2011'=TotalPopEst2011,
           '2012'=TotalPopEst2012,
           '2013'=TotalPopEst2013,
           '2014'=TotalPopEst2014,
           '2015'=TotalPopEst2015,
           '2016'=TotalPopEst2016,
           '2017'=TotalPopEst2017,
           '2018'=TotalPopEst2018)) %>%
  mutate_at(vars(matches("201")),pop_density) %>%
  filter(State == 'TX') %>%
  gather(year,density,-c(FIPS,State,County)) %>%
  mutate(fips = ifelse(nchar(FIPS)<5,paste0("0",FIPS),FIPS)) %>%
  filter(!grepl('000',FIPS),
         !State == 'US')

head(density_df_by_year)
##    FIPS State    County year density  fips
## 1 48001    TX  Anderson 2010      55 48001
## 2 48003    TX   Andrews 2010      10 48003
## 3 48005    TX  Angelina 2010     109 48005
## 4 48007    TX   Aransas 2010      92 48007
## 5 48009    TX    Archer 2010      10 48009
## 6 48011    TX Armstrong 2010       2 48011

The animated highcarter map needs the population densities in a single list called ‘sequence’ in order to work properly. We can create the list of densities by using the list_parse function:

# Create list column containing population densities by year 
density_df_seq <- density_df_by_year %>%
  group_by(fips) %>%
  do(sequence = list_parse(select(., value = density)))

head(density_df_seq)
## Source: local data frame [6 x 2]
## Groups: <by row>
## 
## # A tibble: 6 x 2
##    fips sequence  
##   <int> <list>    
## 1 48001 <list [9]>
## 2 48003 <list [9]>
## 3 48005 <list [9]>
## 4 48007 <list [9]>
## 5 48009 <list [9]>
## 6 48011 <list [9]>

You can see in the output we have a column containing the FIPS codes for each county and a list of length 9 which contains one population density value for each year from 2010 to 2018.

Next, we need to join this dataset back to the original dataset so that we have the county names, years, and population densities all in one dataset:

# Join with original dataset
density_df_by_year <- left_join(density_df_by_year,density_df_seq)
## Joining, by = "fips"
head(density_df_by_year)
##    FIPS State    County year density  fips
## 1 48001    TX  Anderson 2010      55 48001
## 2 48003    TX   Andrews 2010      10 48003
## 3 48005    TX  Angelina 2010     109 48005
## 4 48007    TX   Aransas 2010      92 48007
## 5 48009    TX    Archer 2010      10 48009
## 6 48011    TX Armstrong 2010       2 48011
##                                      sequence
## 1          55, 55, 55, 55, 54, 54, 54, 55, 55
## 2          10, 10, 11, 11, 12, 12, 12, 12, 12
## 3 109, 109, 110, 109, 110, 110, 110, 110, 109
## 4        92, 92, 93, 95, 97, 98, 100, 101, 94
## 5          10, 10, 10, 10, 10, 10, 10, 10, 10
## 6                   2, 2, 2, 2, 2, 2, 2, 2, 2

Creating the Map

To create the animated map, all we need to do is take the existing density_map that we created and update the dataset from density_df to density_df_by_year

# Create interactive map of Texas counties
animated_map <- hcmap(map = "countries/us/us-tx-all", 
      data = density_df_by_year, 
      value = "density", 
      joinBy = c("fips"),
      borderWidth = 0.3)  %>%
  hc_mapNavigation(enabled = TRUE) %>% 
   hc_colorAxis(dataClasses = color_classes(breaks = c(0,10,25,50,100,250,500,1000,2500,max(density_df_by_year$density)),
              colors = brewer.pal(name="YlOrRd",n=9))) %>%
  hc_legend(layout = "vertical", 
            align = "right",
            valueDecimals = 0) %>% 
    hc_tooltip(formatter = JS("function() {
  return ('<br><b>County:</b> ' + this.point.County +
          '<br><b>Population Density:</b> ' + this.point.density + ' people per sq mi'
  )}")) %>%
  hc_title(text = "Population Density by County, 2010 to 2018") %>%
  hc_credits(enabled = TRUE,
             text = "Author: Kyle Cuilla, Data: USDA ERS",
             href = "https://www.ers.usda.gov/data-products/atlas-of-rural-and-small-town-america/download-the-data/")  

And then add the hc_motion option to the map:

# Add animation to map
animated_map <- animated_map %>% 
      hc_motion(enabled = TRUE, 
          series = 0, 
          autoPlay = TRUE,
          loop = TRUE,
          labels = unique(density_df_by_year$year))


Avatar
Kyle Cuilla
Lead Quantitative Analyst
MPS Data Analytics