GIS Tutorial

This tutorial is part of Computerworld’s How to Make a Map with R

nhdatafile <- "NHD2016.xlsx"
nhdatafilecsv <- "NHD2016.csv"
usshapefile <- "cb_2014_us_county_5m/cb_2014_us_county_5m.shp" #keep entire name
nhfipscode <- "33"
scdatafile <- "SCGOP2016.csv"
scfipscode <- "45"

Run any of the install.packages() commands below for packages that are not yet on your system

Load the tmap, tmaptools, and leaflet packages into your working session:

library(tidyverse) 
## Warning: package 'tidyverse' was built under R version 4.2.3
## Warning: package 'purrr' was built under R version 4.2.3
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.0     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.1     ✔ tibble    3.1.8
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(tmap) 
## Warning: package 'tmap' was built under R version 4.2.3
library(tmaptools) 
## Warning: package 'tmaptools' was built under R version 4.2.3
library(leaflet) 
## Warning: package 'leaflet' was built under R version 4.2.3
library(sf) 
## Warning: package 'sf' was built under R version 4.2.3
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(leaflet.extras) 
## Warning: package 'leaflet.extras' was built under R version 4.2.3
library(rio) 
## Warning: package 'rio' was built under R version 4.2.3
library(sp) 
## Warning: package 'sp' was built under R version 4.2.3
getwd()
## [1] "C:/Users/aline/OneDrive/Desktop"

Step 1: Read in the NH election results file:

setwd("C:/Users/aline/Downloads/Gis")
nhdata <- import(nhdatafilecsv) 

Eliminate columns for minor candidates and just use County, Clinton and Sanders columns:

nhdata <- nhdata[,c("County", "Clinton", "Sanders")] 

Step 2: Decide what data to map

nhdata$SandersMarginVotes <- nhdata$Sanders - nhdata$Clinton
nhdata$SandersPct <- (nhdata$Sanders) / (nhdata$Sanders + nhdata$Clinton)

nhdata$ClintonPct <- (nhdata$Clinton) / (nhdata$Sanders + nhdata$Clinton)
nhdata$SandersMarginPCtgPoints <- nhdata$SandersPct - nhdata$ClintonPct
glimpse(nhdata)
## Rows: 10
## Columns: 7
## $ County                  <chr> "Belknap", "Carroll", "Cheshire", "Coos", "Gra…
## $ Clinton                 <int> 3495, 3230, 5132, 2013, 6918, 28147, 12250, 22…
## $ Sanders                 <int> 6005, 5638, 12441, 3639, 14245, 39245, 18107, …
## $ SandersMarginVotes      <int> 2510, 2408, 7309, 1626, 7327, 11098, 5857, 823…
## $ SandersPct              <dbl> 0.6321053, 0.6357691, 0.7079611, 0.6438429, 0.…
## $ ClintonPct              <dbl> 0.3678947, 0.3642309, 0.2920389, 0.3561571, 0.…
## $ SandersMarginPCtgPoints <dbl> 0.2642105, 0.2715381, 0.4159222, 0.2876858, 0.…

Step 3: Get geographic data files

Read in the shapefile for US states and counties:

#install.packages("raster")
#install.packages("rgdal")
setwd("C:/Users/aline/Downloads/Gis")
library(raster) 
## Warning: package 'raster' was built under R version 4.2.3
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
library(rgdal) 
## Warning: package 'rgdal' was built under R version 4.2.3
## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-5, (SVN revision 1199)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/aline/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/aline/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
usgeo <- shapefile("cb_2014_us_county_5m/cb_2014_us_county_5m.shp") 
## Warning: [vect] Z coordinates ignored

Do a quick plot (qtm stands for quick thematic map)of the shapefile and check its structure:

tmap_options(check.and.fix = TRUE)
qtm(usgeo) 
## Warning: The shape usgeo is invalid. See sf::st_is_valid

# (pause to wait for map to render, may take a few seconds)
#view(usgeo) 

Subset just the NH data from the US shapefile

nhgeo <- usgeo[usgeo$STATEFP==nhfipscode,] 

Tmap test plot of the New Hampshire data

qtm(nhgeo)

structure of the object

#str(nhgeo) 
#str(nhdata$County) 
# They're not. Change the county names to plain characters in nhgeo:
nhgeo$NAME <- as.character(nhgeo$NAME) 

Order each data set by county name

nhgeo <- nhgeo[order(nhgeo$NAME),] 
nhdata <- nhdata[order(nhdata$County),] 
# Are the two county columns identical now? They should be:
identical(nhgeo$NAME,nhdata$County) 
## [1] TRUE

Step 4: Merge geo data with results data using the merge function

library(sf) # sf stands for simple features# 
nhmap <- merge(nhgeo, nhdata, by.x = "NAME", by.y = "County") 
# See the new data structure with
#str(nhmap) 

Step 5: Create a static map with tmap’s qtm() function:

qtm(nhmap, "SandersMarginVotes")
## Some legend labels were too wide. These labels have been resized to 0.62, 0.62, 0.62, 0.57, 0.53. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

## Some legend labels were too wide. These labels have been resized to 0.62, 0.62, 0.62, 0.57, 0.53. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger. 
qtm(nhmap, "SandersMarginPCtgPoints")

For more control over look and feel, use the tm_shape()function:

tm_shape(nhmap) + 
 tm_fill("SandersMarginVotes", title="Sanders Margin, Total Votes", palette = "Spectral") + 
 tm_borders(alpha=.5) + 
 tm_text("NAME", size=0.6)
## Some legend labels were too wide. These labels have been resized to 0.62, 0.62, 0.62, 0.57, 0.53. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

Same code as above, but store the static map in avariable, and change the theme to “classic” style:

nhstaticmap <- tm_shape(nhmap) + 
tm_fill("SandersMarginVotes", title="Sanders Margin, Total Votes", palette = "magma") +
#I like viridis #I like viridis too Aline
tm_borders(alpha=.5) + 
tm_text("NAME", size=0.6) + 
tm_style("classic") 

View the map

nhstaticmap 
## Some legend labels were too wide. These labels have been resized to 0.62, 0.62, 0.62, 0.57, 0.53. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output 
## Some legend labels were too wide. These labels have been resized to 0.62, 0.62, 0.62, 0.57, 0.53. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger. 

save the map to a jpg file with tmap’s tmap_save():

tmap_save(nhstaticmap, filename="nhdemprimary.jpg") 
## Map saved to C:\Users\aline\OneDrive\Desktop\nhdemprimary.jpg
## Resolution: 1501.336 by 2937.385 pixels
## Size: 5.004452 by 9.791282 inches (300 dpi)

Part 6

Next up: Code for a basic interactive map, this time for Clinton percentages in NH

Create a palette

clintonPalette <- colorNumeric(palette = "Blues", domain=nhmap$ClintonPct) 

and a pop-up window

library(scales) 
## 
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
## 
##     discard
## The following object is masked from 'package:readr':
## 
##     col_factor
#library(raster) 
nhpopup <- paste0(nhmap$NAME, " County: ",  "Sanders ",
                  percent(nhmap$SandersPct, accuracy = 0.01), " - Clinton ",
                  percent(nhmap$ClintonPct, accuracy = 0.01))

Step 7: Now generate the interactive map:

# re-project
# add the project, WGS84
nhmap_projected <- sp::spTransform(nhmap, "+proj=longlat +datum=WGS84")
leaflet(nhmap_projected) %>% 
  addProviderTiles("CartoDB.Positron") %>%
  addPolygons(stroke=FALSE,
              smoothFactor = 0.2,
              fillOpacity = 0.8,
              popup=nhpopup,
              color= ~clintonPalette(nhmap$ClintonPct)) 

South Carolina data

setwd("C:/Users/aline/Downloads/Gis") 
scdata <- rio::import(scdatafile) 

South Carolina shapefile and Quick plot of scgeo SC geospatial object:

scgeo <- usgeo[usgeo@data$STATEFP=="45",] 
qtm(scgeo) 

Add a column with percent of votes for each candidate. Candidates are in columns 2-7:

#ALine notes from data scdata column2 Jeb Bush - column7 Donald Thump candidates chuck name reads from colnames (data[column numbers]) 

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")
}
winner <- colnames(scdata[2:7])

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

setwd("C:/Users/aline/Downloads/Gis") 
sced <- rio::import("SCdegree.xlsx") 

Check if county names are in the same format in both files

#str(scgeo$NAME) 
## chr [1:46] "Edgefield" "Lee" "Horry" "Allendale" "Marion" "Dorchester" ... 
#str(scdata$County) 
## chr [1:46] "Abbeville" "Aiken" "Allendale" "Anderson" "Bamberg" "Barnwell" ... 
# 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
## [1] TRUE 

Add the election results and rename county column

scmap <- merge(scgeo, scdata, by.x = "NAME", by.y = "County") 

Instead of just coloring the winner, let’s color by strength of win with multiple layers

#Use same intensity for all - get minimum and maximum for the top 3 combined
minpct <- min(c(scdata$`Donald J TrumpPct`, scdata$`Marco RubioPct`, scdata$`Ted CruzPct`)) 
maxpct <- max(c(scdata$`Donald J TrumpPct`, scdata$`Marco RubioPct`, scdata$`Ted CruzPct`)) 

Create leaflet palettes for each layer of the map:

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("#006699", "#ff00cc"), domain = scmap$winner) 
edPalette <- colorNumeric(palette = "Blues", domain=scmap$PctCollegeDegree) 

Create a pop-up:

scpopup <- paste0("<b>County: ", scmap$NAME, "<br />Winner: ", scmap$winner, "<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: ", sced$PctCollegeDegree, "% vs state-wide avg of 25%") 

Add the projection we know from the NH map we’ll need for this data on a Leaflet map:

scmap <- sp::spTransform(scmap, "+proj=longlat +datum=WGS84") 

Basic interactive map showing winner in each county:

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="bottomright", colors=c("#006699", "#ff00cc"),
            labels=c("Trump", "Rubio"))

Put top 3 candidates in their own layers and add education layer, store in scGOPmap variable

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="bottomright", colors=c("#006699", "#ff00cc"), 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(sced$PctCollegeDegree), 
 #this data is in the sced table,not scmaps
 group="College degs") %>%
  
 addLayersControl(
 baseGroups=c("Winners", "Trump", "Rubio", "Cruz", "College degs"),
 position = "bottomright",
 options = layersControlOptions(collapsed = FALSE))
# Now display the map
scGOPmap  %>%
addSearchOSM() 

Save as a self-contained HTML file

htmlwidgets::saveWidget(scGOPmap, file="scGOPwidget.html")
# save as an HTML file with dependencies in another directory:
htmlwidgets::saveWidget(widget=scGOPmap, file="scGOPprimary_withdependencies.html", selfcontained=FALSE, libdir = "js")