Prerequisites

• Intermediate level R
• Basic level HTML
• Basic understanding of how FIPS code is structured
• Basic understanding of how shapefiles are structured

In this tutorial, we will go step by step on how to build the interactive map presented in the burden of asthma due to TRAP study using R. See the maps and associated data sets in the inserts below.

Step 1: Load the required libraries

library(data.table)
library(dplyr)
library(tidyr)
library(RColorBrewer)
library(leaflet)
library(htmltools)
library(htmlwidgets)

Step 2: Set the working directory to the required folder.

setwd("C:/SET/DATA/FOLDER/PATH/HERE")  

In this step we will set the working directory to the location of the files we downloaded, we will also use this directory to store the final interactive map as an HTML file.

Step 3: Read the dataset

load("county_wide.Rda")
load("counties_shp.Rda")

Step 4: Exploring the data sets

> head(county_wide) 
# A tibble: 6 x 8 
# Groups:   GEOID, Pollutant [3]   GEOID Pollutant year  total_children    AC cases     AF  conc   <fct> <fct>     <fct>          <int> <dbl> <dbl>  <dbl> <dbl> 1 01001 NO2       2000           12494  20.9  137. 0.153  12.9  2 01001 NO2       2010           14613  15.2  158. 0.0967  7.92 3 01001 PM10      2000           12494  55.1  137. 0.403  20.9  4 01001 PM10      2010           14613  53.2  158. 0.338  16.4  5 01001 PM2.5     2000           12494  49.8  137. 0.364  15.4  6 01001 PM2.5     2010           14613  36.4  158. 0.231   8.91 

"county_wide" is a data frame with 18,654 rows and 8 columns, each row represents a county. 'GEOID' is the unique identifier for each county, it's taken from the FIPS code in which the first two digits represent the state and the following three digits represent the county. 'Pollutant' is a factor with three levels (NO2, PM10, and PM2.5) we will only be working with 'NO2' in this tutorial (in the next step we will filter out the data). 'AF' is the percentage of childhood asthma incident cases due to the pollutant measured at the county level, we will be mapping 'AF' to the interactive map as our main variable of interest. 'year' is a factor with two levels (2000 and 2010), we will be working with data for the year 2010 only.

Step 5: Filter the "county_wide" data set

county_NO2 <- county_wide %>%     filter(Pollutant == 'NO2', year == 2010) 

We only need rows with "NO2" as the "Pollutant" and rows with "2010" as the "year". Using the pipe operator "%>%" and the filter() function we filter the rows we need then save the filtered data set to a new dataset naming it "county_NO2". The filtered data set "county_NO2" now should contain 3,109 rows and 8 columns.

Step 6: Examining the "counties" shape file

str(counties@data) 'data.frame':	3109 obs. of  10 variables:  $ GEO_ID    : chr  "0500000US01001" "0500000US01009" "0500000US01017" "0500000US01021" ...  $ STATE     : Factor w/ 49 levels "1","4","5","6",..: 1 1 1 1 1 1 1 1 1 1 ...  $ COUNTY    : chr  "001" "009" "017" "021" ...  $ NAME      : chr  "Autauga" "Blount" "Chambers" "Chilton" ...  $ LSAD      : chr  "County" "County" "County" "County" ...  $ CENSUSAREA: num  594 645 597 693 593 ...  $ COUNTYFP  : chr  "001" "009" "017" "021" ...  $ STATEFP   : chr  "01" "01" "01" "01" ...  $ GEOID     : chr  "01001" "01009" "01017" "01021" ...  $ state_name: Factor w/ 51 levels "Alabama","Alaska",..: 1 1 1 1 1 1 1 1 1 1 ...

By using the "@data" we can access the data portion of a shapefile. The "counties" shapefile has 3,109 rows and 17 columns. The "GEOID" column is the unique identifier for each county.

Step 7: Merge the data set "County_NO2" with the Shapefile "counties"

counties@data <- left_join(counties@data,county_NO2, by = "GEOID")  

Step 8: Creating color palettes

# Brewing colors
rc1 <- colorRampPalette(colors = c("#41980a", "#ffe950" ), space = "Lab")(10)
rc2 <- colorRampPalette(colors = c("#ffe950", "#ff4300" ), space = "Lab")(7)
rc3 <- colorRampPalette(colors = c("#ff4300", "#8d1111" ), space = "Lab")(15)
rc4 <- colorRampPalette(colors = c("#8d1111", "black" ), space = "Lab")(8)
rampcols <- c(rc1, rc2, rc3, rc4)
previewColors(colorNumeric(palette = rampcols , domain = NULL), values = 0:55)

AF_pal <- colorNumeric(rampcols, domain = county_NO2$AF*100)

We will manually build a color palette, the colors will be mapped to the distribution of our main variable "AF - the percentage of childhood asthma incident cases due to NO2". We will be using four colors to map the distribution of "AF" using Hex color codes (you can use google to look up the colors i.e. #ffe950), the colors are green, yellow, red and black. The previewColors() function allows you to examine the distribution of colors, in this case, we want to see how are colors are distributed for values between 0 and 55. The color palette is saved to "rampcols". We then use the colorNumeric() function with the color palette we created to assign color values to our main variable "AF", however, here we actually want to assign the values to the percentage so we multiply "AF" with 100. The assigned color values are stored in a variable named "AF_pal"

Step 9: Creating labels

labels <- sprintf(
    "<strong>%s</strong><strong>, %s</strong> <br/>
    Population of children  %g<br/>
    Attributable case %g<br/>
    Percent of all cases %g<br/>
    Mean %s  (ug/m<sup>3</sup>) %g",
    counties@data$NAME,
    counties@data$state_name,
    counties@data$total_children, 
    round(counties@data$AC), 
    round(counties@data$AF*100),
    'NO2',
    round(counties@data$conc, digits = 1)) %>% 
    lapply(htmltools::HTML)
  

label_legend <- HTML("Percentage")
label_title <- HTML('<p><strong><b>Burden of Childhood Asthma Due to NO<sub>2</sub></b>')

In this step we will be creating three things: "labels" that will present county data when we hover the mouse over the map "label_legend" the legend for the map "label_Title" the title of the map The code will allow us to observe specific data when we hover the mouse over a county, the data presented will be the county names "NAME" and the state abbreviation "state_name" separated by and ",", the population of children "total_children", the attributable number of cases "AC", the percentage of childhood asthma cases due to NO2 "AF" and the average NO2 concentration for the county "conc". The sprintf() function allows us to combine formatted text and variable values. while the HTML() function from "htmltools" package marks a given text as HTML format.

Step 10: Creating the interactive map

interactiveMap <- counties %>%
    leaflet(options = leafletOptions(minZoom = 4, maxZoom = 9)) %>% 
    setMaxBounds( -60 , 25 , -130 , 50 ) %>% 
    addProviderTiles(providers$Esri.WorldGrayCanvas) %>%
    addPolygons(data = counties,
                color = "White", weight = 0.15,  smoothFactor = 0.2,
                opacity = 1, fillOpacity = 0.9, 
                fillColor = ~AF_pal(AF*100),
                highlightOptions = highlightOptions(color = "White", weight = 2, bringToFront = TRUE),
                label = ~labels) %>%
    addLegend(pal = AF_pal, values = ~county_NO2$AF*100,
              labFormat = labelFormat(suffix = "%", between = "- "), 
              title = label_legend, position = "bottomright") %>% 
    addControl(label_title, position = "bottomleft")

Using the leaflet() function from the "leaflet" package, we put together a map using the previous data sets, color palettes and labels we created. The map is assigned to a variable "interactiveMap" as a leaflet object.

Step 11: Saving the map as an HTML file

saveWidget(interactiveMap, file = "interactiveMap.html")

Using the saveWidget() function from the "htmlwidgets" package we can transform our leaflet object "interactiveMap" into an HTML object and save it to our working directory. The final map is presented below. Percentage of childhood asthma incident cases due to NO2 from all causes in the year 2010.