Introduction to Mapping Spatial Data in R
Part 5: Interactive Maps
Interactive Maps
# Load packages
library(rgdal)
library(dplyr)
library(ggmap)
library(leaflet)
library(RColorBrewer)
# Set working directory#### Import Data ####
# Options
options(stringsAsFactors = F)
# Geocoded ANC site data
ANC2012 <- read.csv("ANC2012_geocoded.csv")
# Provinces shapefile
provinces_shp <- readOGR("zwe_polbnda_adm1_250k_cso.shp")## OGR data source with driver: ESRI Shapefile
## Source: "zwe_polbnda_adm1_250k_cso.shp", layer: "zwe_polbnda_adm1_250k_cso"
## with 10 features
## It has 5 fieldsBelow is a screenshot of a map of HIV prevalence that I found in the 2015 DHS HIV Fact Sheet, a PDF. Let’s make these data map-able by adding the prevalence as a variable in our provinces shapefile’s data frame.
#### Add HIV prevalence data to provinces ####
# First check the ordering of provinces in the shapefile
provinces_shp$PROVINCE## [1] "Matabeleland North" "Matabeleland South" "Manicaland" "Mashonaland Central" "Masvingo" "Midlands" "Bulawayo" "Mashonaland West" "Mashonaland East" "Harare"# Now create a list of the HIV prevalence for each province, in the order above
HIVprev <- c(0.18, 0.22, 0.11, 0.12, 0.13, 0.15, 0.14, 0.13, 0.15, 0.14)
# Check that prevalence matches correct province
cbind(provinces_shp@data, HIVprev)## PROVINCE PROV_CODE HRName HRPCode HRParent HIVprev
## 0 Matabeleland North 15 Matabeleland North ZW15 ZW 0.18
## 1 Matabeleland South 16 Matabeleland South ZW16 ZW 0.22
## 2 Manicaland 11 Manicaland ZW11 ZW 0.11
## 3 Mashonaland Central 12 Mashonaland Central ZW12 ZW 0.12
## 4 Masvingo 18 Masvingo ZW18 ZW 0.13
## 5 Midlands 17 Midlands ZW17 ZW 0.15
## 6 Bulawayo 10 Bulawayo ZW10 ZW 0.14
## 7 Mashonaland West 14 Mashonaland West ZW14 ZW 0.13
## 8 Mashonaland East 13 Mashonaland East ZW13 ZW 0.15
## 9 Harare 19 Harare ZW19 ZW 0.14# Assign HIVprev as a new column in the provinces data
provinces_shp@data <- cbind(provinces_shp@data, HIVprev)
provinces_shp@data # See that provinces_shp now has an HIVprev variable## PROVINCE PROV_CODE HRName HRPCode HRParent HIVprev
## 0 Matabeleland North 15 Matabeleland North ZW15 ZW 0.18
## 1 Matabeleland South 16 Matabeleland South ZW16 ZW 0.22
## 2 Manicaland 11 Manicaland ZW11 ZW 0.11
## 3 Mashonaland Central 12 Mashonaland Central ZW12 ZW 0.12
## 4 Masvingo 18 Masvingo ZW18 ZW 0.13
## 5 Midlands 17 Midlands ZW17 ZW 0.15
## 6 Bulawayo 10 Bulawayo ZW10 ZW 0.14
## 7 Mashonaland West 14 Mashonaland West ZW14 ZW 0.13
## 8 Mashonaland East 13 Mashonaland East ZW13 ZW 0.15
## 9 Harare 19 Harare ZW19 ZW 0.14Leaflet Widget
#### Interactive Maps using Leaflet ####
# Read the R Documentation
?leaflet# Create a blank Widget in the Viewer tab
leaflet()# Create a Widget with a map of the world
leaflet() %>% addProviderTiles(providers$OpenStreetMap)Interactive Points
#### Point Data ####
# Create a Widget of the ANC data on a map
leaflet(ANC2012) %>% addProviderTiles(providers$OpenStreetMap) %>% addMarkers()# Represent the ANC data as circles
leaflet(ANC2012) %>% addProviderTiles(providers$OpenStreetMap) %>% addCircles(label = ~site)# Represent the ANC data with Circle Markers
leaflet(ANC2012) %>% addProviderTiles(providers$OpenStreetMap) %>%
addCircleMarkers(radius = ~prevalence,
color = c("red"),
label = ~as.character(prevalence),
popup = ~paste0("<b>", name,"</b>", "<br>",
"ANC HIV Prevalence: ", prevalence, "%"
)
)Interactive Polygons
Now let’s experiment with the province polygons.
#### Polygon Data ####
# Try out the provinces
leaflet(provinces_shp) %>% addPolygons(color = "blue", label = ~PROVINCE)# Map of HIV prevalence
leaflet(provinces_shp) %>%
addPolygons(color = c("grey"), weight = 1, smoothFactor = 0.5,
fillColor = ~colorQuantile("YlOrRd", HIVprev)(HIVprev),
opacity = 1.0, fillOpacity = 0.5,
label = ~as.character(HIVprev),
popup = ~PROVINCE
)# Change to equal interval, format to %, add legend
pal <- colorBin("Reds", provinces_shp$HIVprev*100, bins=4, pretty=F)
leaflet(provinces_shp) %>%
addPolygons(color = c("grey"), weight = 0.5, smoothFactor = 0.5,
fillColor = ~pal(HIVprev*100),
opacity = 1.0, fillOpacity = 0.7,
label = ~paste0(PROVINCE, ": ", HIVprev*100, "%"),
popup = ~paste0("<b>", PROVINCE, "</b>", "<br>",
"HIV Prevalence: ", HIVprev*100, "%"),
highlightOptions = highlightOptions(color = "white",
weight = 2,
bringToFront = TRUE)
) %>%
addLegend("bottomleft", pal = pal, values = ~HIVprev*100, opacity = 0.7,
title = "HIV Prevalence (2015)",
labFormat = labelFormat(suffix = "%")
)Layers Control
map <- leaflet(ANC2012) %>%
# Base groups
addProviderTiles(providers$CartoDB) %>%
# Overlay groups
addPolygons(data = provinces_shp, color = c("grey"),
weight = 0.5, smoothFactor = 0.5,
fillColor = ~pal(HIVprev*100),
opacity = 1.0, fillOpacity = 0.7,
label = ~paste0(PROVINCE, ": ", HIVprev*100, "%"),
popup = ~paste0("<b>", PROVINCE, "</b>", "<br>",
"HIV Prevalence: ", HIVprev*100, "%"),
highlightOptions = highlightOptions(color = "white",
weight = 2,
bringToFront = F),
group = "Provinces") %>%
addCircleMarkers(~lon, ~lat, radius = ~prevalence,
color = c("black"), stroke = F, fillOpacity = 0.5,
label = ~paste0(site, ": ", prevalence, "%"),
popup = ~paste0("<b>", name, "</b>", "<br>",
"ANC HIV Prevalence: ", prevalence, "%"),
group = "ANC Sites") %>%
# Layers control
addLayersControl(
overlayGroups = c("ANC Sites", "Provinces"),
options = layersControlOptions(collapsed = FALSE)
)
# Render the "map" Widget
mapAn Interactive HIV Map
#view the distribution of HIV prevalence for sites and provinces
summary(ANC2012$prevalence)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.10 12.60 15.90 16.06 19.30 26.00summary(provinces_shp$HIVprev)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.110 0.130 0.140 0.147 0.150 0.220#create a common color palette for sites and provinces
palette <- "OrRd"
range <- range(7.1,26)
cuts <- c(7.10,13,16,19,26)
pal <- colorBin(palette = palette, domain = range, bins = cuts)
# create the map
map <- leaflet(ANC2012) %>%
# Base groups
addProviderTiles(providers$CartoDB, group = "Canvas") %>%
addProviderTiles(providers$OpenStreetMap, group = "Features") %>%
addProviderTiles(providers$Esri.WorldImagery, group = "Satellite") %>%
# Overlay groups
addPolygons(data = provinces_shp, color = "grey", opacity = 1.0,
weight = 0.5, smoothFactor = 0.5,
fillColor = ~pal(HIVprev*100), fillOpacity = 0.8,
label = ~paste0(PROVINCE,": ", HIVprev*100, "%"),
popup = ~paste0("<b>", PROVINCE,"</b>", "<br>",
"HIV Prevalence: ", HIVprev*100, "%"),
highlightOptions = highlightOptions(color = "grey", weight = 2,
bringToFront = F),
group = "HIV Prevalence by Province (DHS, 2015)"
) %>%
addCircleMarkers(~lon, ~lat, radius = ~prevalence, color = "grey", weight = 1.5,
fillColor = ~pal(prevalence), opacity = 1.0, fillOpacity = 0.8,
label = ~paste0(site, ": ", prevalence, "%"),
popup = ~paste0("<b>", name,"</b>", "<br>",
"Setting: ", classification, "<br>",
"<u>", "HIV among ANC clients, 2012", "</u>", "<br>",
"Total women tested: ", tested.total, "<br>",
"Number tested positive: ", tested.positive, "<br>",
"Prevalence: ", prevalence, "% (CI: ",
prev.ci_lower, "-", prev.ci_upper, ")", "<br>",
"<u>", "Residence of ANC clients", "</u>", "<br>",
"Farm: ", res.farm_pct, "%", "<br>",
"Growth point: ", res.growthpt_pct,"%", "<br>",
"Mine: ", res.mine_pct,"%", "<br>",
"Rural: ", res.rural_pct,"%", "<br>",
"Town or City: ", res.urban_pct, "%"
),
group = "HIV Prevalence by ANC Sentinel Site (2012)") %>%
addLegend(position = "bottomleft", pal = pal, values = ~prevalence, opacity = 0.8,
title = "HIV Prevalence", labFormat = labelFormat(suffix = "%")
) %>%
# Layers control
addLayersControl(
overlayGroups = c("HIV Prevalence by ANC Sentinel Site (2012)",
"HIV Prevalence by Province (DHS, 2015)"),
baseGroups = c("Canvas", "Features", "Satellite"),
options = layersControlOptions(collapsed = TRUE)
)
mapView full screen! Click on the sites and provinces to see additional data!
Exporting to HTML
You can export your interactive map as an HTML document that can be opened outside of R. This can be shared with colleagues, for example via email, without being posted on the internet. The file will open in a web browser, but does not require an internet connection. It can only be viewed be people with access to the HTML file on their computer.
# Save as an HTML file for viewing outside R
library(htmlwidgets)
saveWidget(map, file="Zim_HIV_map.html")If you open your Spatial Workshop folder, you will see the Zim_HIV_map.html file has been saved there. Double-click on the file to open it in a web browser.