Working with maps

Part 1: Map presidential elections results with the maps package

Load the necessary packages and set plotting default

1.Download the elections file from Harvard database

df <- read.csv("./data/1976-2020-president.csv", 1)
colnames(df)

##  [1] "year"             "state"            "state_po"         "state_fips"      
##  [5] "state_cen"        "state_ic"         "office"           "candidate"       
##  [9] "party_detailed"   "writein"          "candidatevotes"   "totalvotes"      
## [13] "version"          "notes"            "party_simplified"

2.Data clean-up and creating new variable “party” that defines color according to the affiliation of the electoral party. Grouping necessary variables: Year-wise and state-wise winning party

win = df %>% 
  dplyr::select(year, state, state_fips, candidate,candidatevotes, party_detailed) %>% 
  group_by(year, state) %>% slice_max(candidatevotes) %>%  
  mutate(party = ifelse(party_detailed == "DEMOCRAT", "blue", 
                        ifelse(party_detailed == "REPUBLICAN", "red", "white")))

Mapping state-wise win in the year 2020

win_20 <- win[win$year == 2020, ]
win_20 <- win_20[match(state.fips$fips, win_20$state_fips),]
map('state', col = win_20$party, fill = TRUE)
legend("bottomright",bty = "n" , legend = c("Democrat", "Republican","Others"), fill = c("blue", "red","white"), cex = 0.7)
title("Presidential elections results by State in 2020", line = 1)

Mapping the election win state-wise over the period of last 44 years(1976 to 2020)

m <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12,13,13,13,13),nrow = 4,byrow = TRUE)
layout(m)


for (i in seq(min(win$year),max(win$year),4)) {
  temp <- win[win$year == i,]
  temp <- temp[match(state.fips$fips, temp$state_fips),]
  map('state', col = temp$party, fill = T)
  title(i, line =1)
}

plot(1, type = "n", axes=FALSE, xlab="", ylab="")


legend(x = "top",inset = 0, 
       legend = c("Democrat", "Republican","Others"), 
       fill = c("blue", "red","white"),
       title = "Presidential elections results by State & Year",
       xpd = TRUE, horiz = TRUE, cex = 1.3, seg.len=0.5, bty = 'n')

Part 2. Interactive Maps with Leaflet

Get familiar with the leaflet package

Esri.NatGeoWorldMap

myMap <- leaflet() %>%
addProviderTiles(providers$Esri.NatGeoWorldMap) %>%
  setView(lat=20, lng=-79, zoom = 7)
myMap

OpenStreetMap

myMap <- leaflet() %>%
addProviderTiles(providers$OpenStreetMap) %>%
  setView(lat=33.947474, lng=-83.373671, zoom = 12)
myMap

OpenTopoMap

myMap <- leaflet() %>%
addProviderTiles(providers$OpenTopoMap) %>%
  setView(lat=33.947474, lng=-100, zoom = 10)
myMap

Add the shapefile for Rohingya refugee camps. Note: Projection is necessary for R to place the coordinates correctly

sp <- shapefile("./data/rrc_al2/220130_RRC_Outline_block_AL2/T220130_RRC_Outline_Block_AL2.shp")
campShapeFile <- spTransform(sp, CRS("+proj=longlat +datum=WGS84 +no_defs"))
head(campShapeFile)

##   OBJECTID_1 OBJECTID Block_Let Camp_SSID Block_Name   Block_SSID SMSD_Cname
## 1          1        1         I   CXB-232     C04X_I CXB-232_I163   Camp 04X
## 2          2        2         B   CXB-232     C04X_B CXB-232_B165   Camp 04X
## 3          3        3         F   CXB-232     C04X_F CXB-232_F161   Camp 04X
## 4          4        4         C   CXB-232     C04X_C CXB-232_C166   Camp 04X
## 5          5        5         E   CXB-232     C04X_E CXB-232_E160   Camp 04X
## 6          6        6         H   CXB-232     C04X_H CXB-232_H162   Camp 04X
##         Camp_Alias         NPM_Cname Area_Acres         CampName
## 1 Camp 4 Extension Camp 04 Extension  17.579304 Camp 4 Extension
## 2 Camp 4 Extension Camp 04 Extension  19.796469 Camp 4 Extension
## 3 Camp 4 Extension Camp 04 Extension   8.892700 Camp 4 Extension
## 4 Camp 4 Extension Camp 04 Extension  40.189147 Camp 4 Extension
## 5 Camp 4 Extension Camp 04 Extension  17.429451 Camp 4 Extension
## 6 Camp 4 Extension Camp 04 Extension   8.238809 Camp 4 Extension
##           Area_SqM  Shape_Leng  Shape_Le_1   Shape_Area
## 1 71140.9196299387 0.012361753 0.012361755 6.192077e-06
## 2 80113.4668634516 0.010098287 0.010098287 6.973411e-06
## 3 35987.4806593324 0.007094379 0.007094378 3.132450e-06
## 4 162639.706686226 0.023266098 0.023266094 1.415641e-05
## 5 70534.4857807752 0.013253804 0.013253802 6.139350e-06
## 6 33341.2781127569 0.006681391 0.006681390 2.901994e-06

Create the interactive map

OpenTopoMap

myMap <- leaflet() %>%
addProviderTiles(providers$OpenTopoMap) %>%
  setView(lng = 92.14871, lat = 21.18780, zoom = 12) %>%
  addPolygons(data=campShapeFile, fill=TRUE, stroke=T, weight=1, highlight = highlightOptions(fillOpacity = 0.7), label = campShapeFile$Block_No)

myMap

OpenStreetMap

myMap <- leaflet() %>%
addProviderTiles(providers$OpenStreetMap) %>%
  setView(lng = 92.14871, lat = 21.18780, zoom = 12) %>%
  addPolygons(data=campShapeFile, fill=TRUE, stroke=T, weight=1, highlight = highlightOptions(fillOpacity = 0.7), label = campShapeFile$Block_No)

myMap

Add tiles from the web

Overlayed a raster image of current rainfall in the USA. It updates every few minutes (almost on a real time bases).

mapUSA <- leaflet() %>%
addProviderTiles(providers$OpenStreetMap.Mapnik) %>%
  setView(lat=33.947474, lng=-83.373671, zoom = 12) %>%
  addWMSTiles("http://mesonet.agron.iastate.edu/cgi-bin/wms/nexrad/n0r.cgi",
layers = "nexrad-n0r-900913",
options = WMSTileOptions(format = "image/png", transparent = TRUE),
attribution = "Weather data 2012 IEM Nexrad"
)

mapUSA

For more excellent leaflet features

Reference and disclaimer: This is an attempt to reproduce parts of the article from github repository for the purpose of learning excellent leaflet features.

states <- geojsonio::geojson_read("https://rstudio.github.io/leaflet/json/us-states.geojson", what = "sp")
#class(states)
names(states)

## [1] "id"      "name"    "density"

Basic states map

m <- leaflet(states) %>%
  setView(-96, 37.8, 4) %>%
  addProviderTiles("MapBox", options = providerTileOptions(
    id = "mapbox.light",
    accessToken = Sys.getenv('MAPBOX_ACCESS_TOKEN')))
m %>% addPolygons()

Adding colour

bins <- c(0, 10, 20, 50, 100, 200, 500, 1000, Inf)
pal <- colorBin("YlOrRd", domain = states$density, bins = bins)
m %>% addPolygons(
  fillColor = ~pal(density),
  weight = 2,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7)

Adding Interaction

m %>% addPolygons(
  fillColor = ~pal(density),
  weight = 2,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlightOptions = highlightOptions(
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE))

Adding Custom info

labels <- sprintf(
  "<strong>%s</strong><br/>%g people / mi<sup>2</sup>",
  states$name, states$density
) %>% lapply(htmltools::HTML)

m %>% addPolygons(
  fillColor = ~pal(density),
  weight = 2,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlightOptions = highlightOptions(
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto"))

Final interactive map with legend

states <- geojsonio::geojson_read("https://rstudio.github.io/leaflet/json/us-states.geojson", what = "sp")

bins <- c(0, 10, 20, 50, 100, 200, 500, 1000, Inf)
pal <- colorBin("YlOrRd", domain = states$density, bins = bins)

labels <- sprintf(
  "<strong>%s</strong><br/>%g people / mi<sup>2</sup>",
  states$name, states$density
) %>% lapply(htmltools::HTML)

leaflet(states) %>%
  setView(-96, 37.8, 4) %>%
  addProviderTiles("MapBox", options = providerTileOptions(
    id = "mapbox.light",
    accessToken = Sys.getenv('MAPBOX_ACCESS_TOKEN'))) %>%
  addPolygons(
    fillColor = ~pal(density),
    weight = 2,
    opacity = 1,
    color = "white",
    dashArray = "3",
    fillOpacity = 0.7,
    highlightOptions = highlightOptions(
      weight = 5,
      color = "#666",
      dashArray = "",
      fillOpacity = 0.7,
      bringToFront = TRUE),
    label = labels,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "15px",
      direction = "auto")) %>%
  addLegend(pal = pal, values = ~density, opacity = 0.7, title = NULL,
    position = "bottomright")

Working with maps

Presidential elections

Cliff

11/03/2022

Part 1: Map presidential elections results with the maps package

Part 2. Interactive Maps with Leaflet

Esri.NatGeoWorldMap

OpenStreetMap

OpenTopoMap

OpenTopoMap

OpenStreetMap

Basic states map

Adding colour

Adding Interaction

Adding Custom info