tmap Tutorial

library(devtools)
#install_github("mtennekes/tmaptools")
#install_github("mtennekes/tmap")
#devtools::install_github('walkerke/tigris')
library("tmap")
library("tmaptools")
library(tigris)
library(tidyverse)
#install.packages("rio")
library(rio)
library(stringr)

To import a shapefile of U.S. states into R with tigris as a simple features object, just run the command:

us_geo <- states(class = "sf")

trying URL 'http://www2.census.gov/geo/tiger/TIGER2015/STATE/tl_2015_us_state.zip'
Content type 'application/zip' length 8780814 bytes (8.4 MB)
==================================================
downloaded 8.4 MB

us_geo%>%data.frame()%>%head() 

Note that states() without class = “sf” will return a Spatial Polygons Data Frame and not a simple features object.

Then read the geographic information into an R object called us_geo with

#Alternatively
#download.file("http://www2.census.gov/geo/tiger/GENZ2015/shp/cb_2015_us_state_20m.zip", destfile = "states.zip")
#unzip("states.zip")
#us_geo <-read_shape("cb_2015_us_state_20m.shp", as.sf = TRUE, stringsAsFactors = FALSE)

Import the data you want to map. I went to the U.S. Bureau of Labor Statistics Occupational Employment Statistics query page, chose “One occupation for multiple geographical areas,” and then picked “Computer and Information Systems Managers” by geographic type State, “All states in this list,” “Annual median wage” as the data type, and Excel as the output. If you open the file, you’ll see there are 3 rows of meta data at the top (R will count merged rows 2, 3 and 4 as a single row) and 4 rows at the bottom. I moved the footnotes to a separate tab and skipped the first 3 rows when importing (there’s no easy way to skip rows at the end).

#load excel file with rio
Employment<- rio::import("/Users/nanaakwasiabayieboateng/Documents/memphisclassesbooks/DataMiningscience/Maps/OES_Report.xlsx", skip = 3)
names(Employment)

[1] "Period:May 2016" "X__1"           

ncol(Employment)

[1] 2

names(Employment) <-c("Area", "Employment")
Employment=Employment%>%slice(-c(1,2))
Employment

More importantly, the wage column imports as character strings, not numbers, I’m guessing because of extra spaces in the data. I trimmed white space and then converted Median.Wages into numbers with

# remove whitespaces on the both sides
Employment$Employment=stringr::str_trim(Employment$Employment, side = "both")
#Alternatively
#1 Employment$Employment <- trimws(Employment$Employment)
#2 Employment$State <- str_sub(Employment$Area, 1, nchar(Employment$Area) - 9)
Employment$Employment<-as.numeric(Employment$Employment)
Employment%>%head()

There are several ways of getting rid of the (01000000) and similar text in the state names. The most robust way is with a regular expression – delete anything that’s not an alpha character. This code does so:

#  Employment$State <-stringr::str_replace_all(Employment$Area, pattern = "[^[:alpha:]]", "")
# 
# Employment$State%>%head()
# 
# #Alternatively
# #Employment$State <-gsub("[^[:alpha:]]", "", Employment$Area)
# 
#  Employment%>%head()
 
 Employment= Employment%>%mutate(State=stringr::str_replace_all(Area, pattern = "[^[:alpha:]]", ""))%>%select(-Area) 
 
  Employment%>%head()

 Employment= Employment%>%filter(State!="NA")%>%slice(1:54)
 Employment%>%head()

NA

Merge (join) the geospatial and data files. This is easy with the tmaptools package’s append_data() function

 Employmentmap<-append_data(us_geo,  Employment, key.shp = "NAME", key.data = "State")

Under coverage: 15 out of 56 shape features did not get appended data. Run under_coverage() to get the corresponding feature id numbers and key values.
Over coverage: 13 out of 54 data records were not appended. Run over_coverage() to get the corresponding data row numbers and key values.

 Employmentmap%>%data.frame()%>%head() 

Finally, Create your map. This part is incredibly easy with tmap’s qtm() (quick theme map) function:

qtm( Employmentmap, fill = "Employment")

tm_shape(Employmentmap) +
  tm_polygons("Employment")

To visualize a bigger map without Alaska

contig_48 <- subset(Employmentmap, !(NAME %in% c("Alaska", "Hawaii", "Puerto Rico")))

tm_shape(contig_48) +
  tm_polygons("Employment", id = "NAME") 

To make this interactive, you just need to switch tmap’s mode from “plot,” which is static, to “view”, which is interactive, using the tmap_mode() function:

tmap_mode("view")

tmap mode set to interactive viewing

interactive_map <- tm_shape(Employmentmap) +
  tm_polygons("Employment", id = "NAME")
interactive_map<-tmap_leaflet(interactive_map)
interactive_map %>%
leaflet::setView(-96, 37.8, zoom = 3)

With the htmlwidgets package, save the map with saveWidget(my_interactive_map, “mymap.html”). You can see a version of the map below, which I centered and zoomed using leaflet’s setView() function after loading the leaflet package:

interactive_map<-interactive_map %>%
leaflet::setView(-96, 37.8, zoom = 4)
htmlwidgets::saveWidget(interactive_map, "mymap.html")

unemp <- read.csv("http://datasets.flowingdata.com/unemployment09.csv",
                  header = FALSE, stringsAsFactors = FALSE)
names(unemp) <- c("id", "state_fips", "county_fips", "name", "year",
                  "?", "?", "?", "rate")
unemp$county <- tolower(gsub(" County, [A-Z]{2}", "", unemp$name))
unemp$county <- gsub("^(.*) parish, ..$","\\1", unemp$county)
unemp$state <- gsub("^.*([A-Z]{2}).*$", "\\1", unemp$name)
unemp%>%head()

county_df <- map_data("county", projection = "albers", parameters = c(39, 45))
names(county_df) <- c("long", "lat", "group", "order", "state_name", "county")
county_df$state <- state.abb[match(county_df$state_name, tolower(state.name))]
county_df$state_name <- NULL
state_df <- map_data("state", projection = "albers", parameters = c(39, 45))
choropleth <- merge(county_df, unemp, by = c("state", "county"))
choropleth <- choropleth[order(choropleth$order), ]
choropleth %>%head()

ggplot(choropleth, aes(long, lat, group = group)) +
  geom_polygon(aes(fill = rate), colour = alpha("white", 1 / 2), size = 0.2) +
  geom_polygon(data = state_df, colour = "white", fill = NA) +
  coord_fixed() +
  theme_minimal() +
  ggtitle("US unemployment rate by county") +
  theme(axis.line = element_blank(), axis.text = element_blank(),
        axis.ticks = element_blank(), axis.title = element_blank()) +
  scale_fill_viridis(option="magma")