Introduction to Mapping Spatial Data in R

Part 5: Interactive Maps

<< Workshop Home

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")

Below 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.

http://dhsprogram.com/publications/publication-hf58-hiv-fact-sheets.cfm

#### Add HIV prevalence data to provinces ####

# First check the ordering of provinces in the shapefile
provinces_shp$PROVINCE

# 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)

# 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

Leaflet 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
map

An Interactive HIV Map

#view the distribution of HIV prevalence for sites and provinces
summary(ANC2012$prevalence)
summary(provinces_shp$HIVprev)

#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)
  )
map

View 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.

<< Workshop Home