Geographic Information Science and Systems

library(tmap)
library(tmaptools)
library(sf)
library(leaflet)
library(dplyr)
library(raster)
library(rio)
library(RColorBrewer)
library(scales)
library(leaflet.extras)

Step 1: Get election results data

nhdatafile <- "data/NHD2016.xlsx"
nhdata <- rio::import(nhdatafile)
nhdata <- nhdata[,c("County", "Clinton", "Sanders")]

Step 2: Decide what data to map

# Add columns for percents and margins
nhdata$SandersMarginVotes <- nhdata$Sanders - nhdata$Clinton
nhdata$SandersPct <- (nhdata$Sanders - nhdata$Clinton) / (nhdata$Sanders + nhdata$Clinton) # Will use formatting later to multiply by a hundred 
nhdata$ClintonPct <- (nhdata$Clinton - nhdata$Sanders) / (nhdata$Sanders + nhdata$Clinton)
nhdata$SandersMarginPctgPoints <- nhdata$SandersPct - nhdata$ClintonPct

Step 3: Get your geographic data

usshapefile <- "data/cb_2014_us_county_5m/cb_2014_us_county_5m.shp"
usgeo <- read_shape(file=usshapefile, as.sf = TRUE)
#qtm(usgeo)
#str(usgeo)
nhgeo <- filter(usgeo, STATEFP=="33")
qtm(nhgeo)

Step 4: Merge spatial and results data

str(nhgeo$NAME)
##  Factor w/ 1921 levels "Abbeville","Acadia",..: 684 791 416 138 1470 334 1653 1131 282 1657
str(nhdata$County)
##  chr [1:10] "Belknap" "Carroll" "Cheshire" "Coos" "Grafton" ...
nhgeo$NAME <- as.character(nhgeo$NAME)
nhgeo <- nhgeo[order(nhgeo$NAME),]
nhdata <- nhdata[order(nhdata$County),]
identical(nhgeo$NAME,nhdata$County )
## [1] TRUE
nhmap <- append_data(nhgeo, nhdata, key.shp = "NAME", key.data="County")
str(nhmap)
## Classes 'sf' and 'data.frame':   10 obs. of  16 variables:
##  $ STATEFP                : Factor w/ 56 levels "01","02","04",..: 30 30 30 30 30 30 30 30 30 30
##  $ COUNTYFP               : Factor w/ 328 levels "001","003","005",..: 1 2 3 5 6 8 10 12 14 15
##  $ COUNTYNS               : Factor w/ 3233 levels "00023901","00025441",..: 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505
##  $ AFFGEOID               : Factor w/ 3233 levels "0500000US01001",..: 1765 1766 1767 1768 1769 1770 1771 1772 1773 1774
##  $ GEOID                  : Factor w/ 3233 levels "01001","01003",..: 1765 1766 1767 1768 1769 1770 1771 1772 1773 1774
##  $ NAME                   : chr  "Belknap" "Carroll" "Cheshire" "Coos" ...
##  $ LSAD                   : Factor w/ 11 levels "00","03","04",..: 5 5 5 5 5 5 5 5 5 5
##  $ ALAND                  : Factor w/ 3233 levels "1000508842","1001064387",..: 67 1991 1424 2526 2494 1831 1997 1382 3172 730
##  $ AWATER                 : Factor w/ 3233 levels "0","10017640",..: 778 638 2460 3097 55 1942 2454 1327 1852 1817
##  $ Clinton                : num  3495 3230 5132 2013 6918 ...
##  $ Sanders                : num  6005 5638 12441 3639 14245 ...
##  $ SandersMarginVotes     : num  2510 2408 7309 1626 7327 ...
##  $ SandersPct             : num  0.264 0.272 0.416 0.288 0.346 ...
##  $ ClintonPct             : num  -0.264 -0.272 -0.416 -0.288 -0.346 ...
##  $ SandersMarginPctgPoints: num  0.528 0.543 0.832 0.575 0.692 ...
##  $ geometry               :sfc_POLYGON of length 10; first list element: List of 1
##   ..$ : num [1:33, 1:2] -71.7 -71.7 -71.7 -71.7 -71.7 ...
##   ..- attr(*, "class")= chr  "XY" "POLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
##   ..- attr(*, "names")= chr  "STATEFP" "COUNTYFP" "COUNTYNS" "AFFGEOID" ...

Step 5: Create a static map

qtm(nhmap, "SandersMarginVotes")

qtm(nhmap, "SandersMarginPctgPoints")

tm_shape(nhmap) +
tm_fill("SandersMarginVotes", title="Sanders Margin, Total Votes", palette = "PRGn") +
tm_borders(alpha=.5) +
tm_text("NAME", size=0.8)

tm_shape(nhmap) +
  tm_fill("SandersMarginVotes", title="Sanders Margin, Total Votes", palette = "PRGn") +
  tm_borders(alpha=.5) +
  tm_text("NAME", size=0.8) + 
tm_style_classic()

nhstaticmap <- tm_shape(nhmap) +
  tm_fill("SandersMarginVotes", title="Sanders Margin, Total Votes", palette = "PRGn") +
  tm_borders(alpha=.5) +
tm_text("NAME", size=0.8)
save_tmap(nhstaticmap, filename="nhdemprimary.jpg")

Step 6: Create palette and pop-ups for interactive map

Create a Leaflet palette with this syntax:

#mypalette <- colorFunction(palette = "colors I want", domain = mydataframe$dataColumnToMap)
clintonPalette <- colorNumeric(palette = "Blues", domain=nhmap$ClintonPct)
#install.packages("scales")
nhpopup <- paste0("County: ", nhmap$County,
"Sanders ", percent(nhmap$SandersPct), " - Clinton ", percent(nhmap$ClintonPct))
nhmap <- rename(nhmap, County = NAME)

Step 7: Generate an interactive map

Map code:

leaflet(nhmap) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(stroke=FALSE, 
              smoothFactor = 0.2,
              fillOpacity = .8, 
              popup=nhpopup,
              color= ~clintonPalette(nhmap$ClintonPct)
)

Step 8: Add palettes for a multi-layer map

scdatafile <- "data/SCGOP2016.csv"
scfipscode <- "45"
# South Carolina data
scdata <- rio::import(scdatafile)

# South Carolina shapefile:
scgeo <- dplyr::filter(usgeo, STATEFP==scfipscode)

# Quick plot of scgeo SC geospatial object:
qtm(scgeo)

# Add a column with percent of votes for each candidate. Candidates are in columns 2-7:
candidates <- colnames(scdata[2:7])
for(i in 2:7){
  j = i + 7
  temp <- scdata[[i]] / scdata$Total
  scdata[[j]] <- temp
  colnames(scdata)[j] <- paste0(colnames(scdata)[i], "Pct")
}  
  
# Get winner in each precinct
for(i in 1:nrow(scdata)){
  scdata$winner[i] <- names(which.max(scdata[i,2:7]))
}

# Import spreadsheet with percent of adult population holding at least a 4-yr college degree
sced <- rio::import("data/SCdegree.xlsx")


# Check if county names are in the same format in both files
str(scgeo$NAME)
##  Factor w/ 1921 levels "Abbeville","Acadia",..: 554 995 810 35 1073 523 1662 359 100 331 ...
str(scdata$County)
##  chr [1:46] "Abbeville" "Aiken" "Allendale" "Anderson" "Bamberg" ...
# Change the county names to plain characters in scgeo:
scgeo$NAME <- as.character(scgeo$NAME)

# Order each data set by county name
scgeo <- scgeo[order(scgeo$NAME),]
scdata <- scdata[order(scdata$County),]

# Are the two county columns identical now? They should be:
identical(scgeo$NAME,scdata$County )
## [1] TRUE
# Add the election results and rename county column
scmap <- append_data(scgeo, scdata, key.data = "County", key.shp = "NAME")
scmap <- rename(scmap, County = NAME)
scmap <- append_data(scmap, sced, key.shp = "County", key.data = "County")
minpct <- min(c(scmap$Donald.J.TrumpPct, scmap$Marco.RubioPct , scmap$Ted.CruzPct))
maxpct <- max(c(scmap$Donald.J.TrumpPct, scmap$Marco.RubioPct , scmap$Ted.CruzPct))
trumpPalette <- colorNumeric(palette = "Purples", domain=c(minpct, maxpct))
rubioPalette <- colorNumeric(palette = "Reds", domain = c(minpct, maxpct))
cruzPalette <- colorNumeric(palette = "Oranges", domain = c(minpct, maxpct))
winnerPalette <- colorFactor(palette=c("#984ea3", "#e41a1c"), domain = scmap$winner)
edPalette <- colorNumeric(palette = "Blues", domain=scmap$PctCollegeDegree)
scpopup <- paste0("<b>County: ", scmap$County, "<br />Winner: ", scmap$winner, "</b><br /><br />Trump: ", percent(scmap$Donald.J.TrumpPct), "<br />Rubio: ", percent(scmap$Marco.RubioPct), "<br />Cruz: ", percent(scmap$Ted.CruzPct), "<br /><br />Pct w college ed: ", scmap$PctCollegeDegree, "% vs state-wide avg of 25%")
scmap <- sf::st_transform(scmap, "+proj=longlat +datum=WGS84")
leaflet(scmap) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(stroke=TRUE,
              weight=1,
              smoothFactor = 0.2,
              fillOpacity = .75,
              popup=scpopup, 
              color= ~winnerPalette(scmap$winner),
              group="Winners"
 ) %>%
    addLegend(position="bottomleft", colors=c("#984ea3", "#e41a1c"), labels=c("Trump", "Rubio"))

Step 9: Add map layers and controls

scGOPmap <- leaflet(scmap) %>%
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(stroke=TRUE,
              weight=1,
              smoothFactor = 0.2,
              fillOpacity = .75,
              popup=scpopup, 
              color= ~winnerPalette(scmap$winner),
              group="Winners"
  ) %>% 
    addLegend(position="bottomleft", colors=c("#984ea3", "#e41a1c"), labels=c("Trump", "Rubio"))  %>%

  addPolygons(stroke=TRUE,
     weight=1,
     smoothFactor = 0.2, 
     fillOpacity = .75, 
     popup=scpopup, 
     color= ~trumpPalette(scmap$Donald.J.TrumpPct),
     group="Trump"
    ) %>%

  addPolygons(stroke=TRUE,
              weight=1,
              smoothFactor = 0.2, 
              fillOpacity = .75, 
              popup=scpopup, 
              color= ~rubioPalette(scmap$Marco.RubioPct),
              group="Rubio"
  ) %>%

  addPolygons(stroke=TRUE,
              weight=1,
              smoothFactor = 0.2, 
              fillOpacity = .75, 
              popup=scpopup, 
              color= ~cruzPalette(scmap$Ted.CruzPct),
              group="Cruz"
  ) %>%

  addPolygons(stroke=TRUE,
              weight=1,
              smoothFactor = 0.2, 
              fillOpacity = .75, 
              popup=scpopup, 
              color= ~edPalette(scmap$PctCollegeDegree),
              group="College degs"
  ) %>%

  addLayersControl(
      baseGroups=c("Winners", "Trump", "Rubio", "Cruz", "College degs"),
      position = "bottomleft",
      options = layersControlOptions(collapsed = FALSE)
      ) 
scGOPmap

Step 10: Save your interactive map

library("htmlwidgets")
saveWidget(widget=scGOPmap, file="scGOPprimary.html")

Extra: Add address search to an R map

scGOPmap  %>% addSearchOSM() %>% addSearchGoogle()