addlayersControl displays only one base layer

jacmac · April 30, 2020, 2:01pm

Hello R Studio Community

I am very new to R and I am trying to create an interactive choropleth map. I have created two maps, one with 2020 data and one with 2016 data and would like to display both using the addLayersControl function. However when I click on the radio button, it only displays data for one year for each map instance.

2020 Map

2016 Map

I'm obviously going wrong somewhere and have looked at similar examples, without success. I would greatly appreciate any guidance as to where I am going wrong please.

Many thanks

Jacqueline

# Load required packages to handle shapefiles, transform the data, import from
# Excel and create the map. 
library(SPARQL)
library(dplyr)
library(readxl)
library(leaflet)
library(rgdal)

# Download the Local Authority District (LAD) shapefile. 
# NOTE: I store it in a local folder. You have to change that if needed.
download.file("https://opendata.arcgis.com/datasets/fab4feab211c4899b602ecfbfbc420a3_4.zip?outSR=%7B%22wkid%22%3A27700%2C%22latestWkid%22%3A27700%7D",destfile="LAD.zip")

# Unzip this file. You can do it with R (as below), or clicking on the object you downloaded.
system("unzip C:/Users/Jac/Downloads/Local_Authority_Districts_December_2017_Ultra_Generalised_Clipped_Boundaries_in_United_Kingdom_WGS84.zip")

# Read in historical data from Excel file
exceldata <- read_excel('C:/Users/Jac/Documents/BSc Computing Year 4 Pt2/Dissertation/RTest/Maps/SIMDAllTest.xlsx')

# SPARQL query variable to retrieve data zones and their respective SIMD rank values
query1 <- 'PREFIX qb: <http://purl.org/linked-data/cube#>
PREFIX rdfs: <http://www.w3.org/2000/01/rdf-schema#>
PREFIX sdmx: <http://purl.org/linked-data/sdmx/2009/concept#>
PREFIX data: <http://statistics.gov.scot/data/>
PREFIX sdmxd: <http://purl.org/linked-data/sdmx/2009/dimension#>
PREFIX mp: <http://statistics.gov.scot/def/measure-properties/>
PREFIX stat: <http://statistics.data.gov.uk/def/statistical-entity#>
PREFIX skos: <http://www.w3.org/2004/02/skos/core#>
SELECT ?dataZone ?SIMDrank
WHERE {
    ?indicator qb:dataSet data:scottish-index-of-multiple-deprivation;
              <http://statistics.gov.scot/def/dimension/simdDomain> <http://statistics.gov.scot/def/concept/simd-domain/simd>;
              mp:rank ?SIMDrank;
              sdmxd:refPeriod <http://reference.data.gov.uk/id/year/2020> ;
              sdmxd:refArea ?area.
    ?area rdfs:label ?dataZone.
}'

# SPARQL endpoint to retrieve the data
endpoint <- "http://statistics.gov.scot/sparql"

# Assign output of SPARQL query to 'qddata'
qddata <- SPARQL(endpoint, query1)

# Assign results of SPARQL query to data frame 'SIMDrank'
SIMDrank <- qddata$results

# SPARQL query to retrieve data zones, council areas and council area codes
query2 <- 'PREFIX qb: <http://purl.org/linked-data/cube#>
PREFIX rdfs: <http://www.w3.org/2000/01/rdf-schema#>
PREFIX sdmx: <http://purl.org/linked-data/sdmx/2009/concept#>
PREFIX data: <http://statistics.gov.scot/data/>
PREFIX sdmxd: <http://purl.org/linked-data/sdmx/2009/dimension#>
PREFIX mp: <http://statistics.gov.scot/def/measure-properties/>
PREFIX stat: <http://statistics.data.gov.uk/def/statistical-entity#>
PREFIX skos: <http://www.w3.org/2004/02/skos/core#>
SELECT ?dataZone ?councilArea ?councilAreaCode 
WHERE {
    ?dz <http://statistics.gov.scot/def/hierarchy/best-fit#council-area> ?ca.
    ?ca rdfs:label ?councilArea.
    ?ca <http://publishmydata.com/def/ontology/foi/code> ?councilAreaCode. 
    ?dz rdfs:label ?dataZone.
}'

# Assign output of SPARQL query to 'qddata2'
qddata2 <- SPARQL(endpoint, query2)

# Assign results of SPARQL query to data frame 'geo_lkp'
geo_lkp <- qddata2$results

# Join the 2 data frames to link SIMD to council areas
SIMD_ca_2020 <- inner_join(SIMDrank, geo_lkp, by="dataZone")

# Join the new data frame to the Excel data for 2016
SIMD_ca_2016 <- inner_join(exceldata, geo_lkp, by="dataZone")

# Calculate mean SIMD rank per council area for 2020
SIMD_mean_2020 <- SIMD_ca_2020 %>% 
  group_by(councilAreaCode, councilArea) %>% 
  summarise(meanSIMDrank=mean(SIMDrank))

# Calculate mean SIMD rank per council area for 2016
SIMD_mean_2016 <- SIMD_ca_2016 %>% 
  group_by(councilAreaCode, councilArea) %>% 
  summarise(meanSIMDrank2=mean(SIMD2016))

# Load shapefile into R as spatial polygons data frame
boundary <- readOGR(dsn="C:/Users/Jac/Documents/BSc Computing Year 4 Pt2/Dissertation/RTest/Maps", layer="Local_Authority_Districts_December_2017_Ultra_Generalised_Clipped_Boundaries_in_United_Kingdom_WGS84")

# Merge spatial polygons data frame with data frame containing mean for 2020
merged_2020 <- merge(boundary, SIMD_mean_2020, by.x = "lad17nm", 
                by.y = "councilArea", all.x = FALSE)

# Merge spatial polygons data frame with data frame containing mean for 2016
merged_2016 <- merge(boundary, SIMD_mean_2016, by.x = "lad17nm", 
                     by.y = "councilArea", all.x = FALSE)

# Project 2020 data to World Geodetic System 1984 using spTransform to ensure correct plotting
project_2020 <- spTransform(merged_2020, 
                         CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

# Project 2016 data to World Geodetic System 1984 using spTransform to ensure correct plotting
project_2016 <- spTransform(merged_2016, 
                            CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))

# Create bins and palette for mean SIMD rank
bins <- c(2000, 2500, 3000, 3500, 4000, 4500, 5000, 6000)
pal2020 <- colorBin("YlOrRd", domain = project_2020$meanSIMDrank, bins = bins)
pal2016 <- colorBin("Blues", domain = project_2016$meanSIMDrank, bins = bins)

# Plot mean 2020 SIMD rank for each council area
map2020 <- leaflet(project_2020) %>% 
  addProviderTiles("CartoDB.Positron", 
                   options= providerTileOptions(opacity = 0.99)) %>% 
  addPolygons(fillColor = ~pal2020(meanSIMDrank),
              weight = 2,
              opacity = 1,
              color = "white",
              dashArray = "3",
              fillOpacity = 0.7,
              highlight = highlightOptions(
                weight = 2,
                color = "#666",
                dashArray = "",
                fillOpacity = 0.7,
                bringToFront = TRUE),
              group = "SIMD2020",
              label=~paste(project_2020$lad17nm, 
                           round(project_2020$meanSIMDrank)),
              labelOptions = labelOptions(textsize = "15px",
                                          direction = "auto")) %>% 
  addLegend(pal = pal2020, 
            values = ~meanSIMDrank, 
            opacity = 0.7, 
            title = "Mean SIMD Rank 2020",
            position = "bottomright")

addLayersControl(map2020,
                 baseGroups = c("SIMD2016", "SIMD2020"),
                 position = "bottomleft",
                 options = layersControlOptions(collapsed = FALSE)
)

# Plot mean 2016 SIMD rank for each council area
map2016 <- leaflet(project_2016) %>% 
  addProviderTiles("CartoDB.Positron", 
                   options= providerTileOptions(opacity = 0.99)) %>% 
  addPolygons(fillColor = ~pal2016(meanSIMDrank2),
              weight = 2,
              opacity = 1,
              color = "white",
              dashArray = "3",
              fillOpacity = 07,
              highlight = highlightOptions(
                weight = 2,
                color = "#666",
                dashArray = "",
                fillOpacity = 0.7,
                bringToFront = TRUE),
              group = "SIMD2016",
              label=~paste(project_2016$lad17nm, 
                           round(project_2016$meanSIMDrank2)),
              labelOptions = labelOptions(textsize = "15px",
                                          direction = "auto"))%>% 
  addLegend(pal = pal2016, 
            values = ~meanSIMDrank2, 
            opacity = 0.7, 
            title = "Mean SIMD Rank 2016",
            position = "bottomright")

addLayersControl(map2016,
  baseGroups = c("SIMD2016", "SIMD2020"),
  position = "bottomleft",
  options = layersControlOptions(collapsed = FALSE)
  )

jlacko · April 30, 2020, 3:04pm

Hi @jacmac , and welcome!

The example you are sharing is not exactly reproducible (you know what is liable to happen to you for using path like C:/Users/Jac/Documents/ ? your on !!! ).

So just some general comments:

you are creating two leaflet objects - map2020 and map2016; this is not how addLayersControl works; you are supposed to create two polygon layers in one leaflet object, and use the layers control to switch between them (via baseGroup = radio button behaviour) or turn each on and off (via overlayGroups = checkbox behaviour).
while the addLayersControl call will work separately when used with a map as a first argument it is not a particularly good practice to do so; it is much better to attach it via pipe (%>%) to an existing leaflet pipeline.

For an example of a reproducible layers control call consider this code, taken from my Leaflet in R walkhrough.

Note how two groups are defined in two separate addCircleMarkers calls (the rest of the map will remain the same) and how they are utilized in the addLayersControls.

For points it makes sense to add the groups one by one (checkbox like) but polygons will likely require radio button functrionality (baseGroups instead of overlayGroups)

library(sf)             # for working with spatial data
library(dplyr)          # data frame manipulation
library(leaflet)        # because leaflet :)
library(htmltools)      # tools to support html workflow
library(leaflet.extras) # extending the leaflet.js


# set up a data frame of points
body <- data.frame(name = c("Kramářova vila", "Pražský hrad", "Strakova akademie", "Nejvyšší soud", "Ústavní soud", "Sněmovna", "Senát"),
                   branch = c("executive", "executive", "executive", "judiciary", "judiciary", "legislature", "legislature"),
                   link = c("https://en.wikipedia.org/wiki/Kram%C3%A1%C5%99%27s_Villa",
                            "https://en.wikipedia.org/wiki/Prague_Castle",
                            "https://en.wikipedia.org/wiki/Straka_Academy",
                            "https://en.wikipedia.org/wiki/Supreme_Court_of_the_Czech_Republic",
                            "https://en.wikipedia.org/wiki/Constitutional_Court_of_the_Czech_Republic",
                            "https://en.wikipedia.org/wiki/Chamber_of_Deputies_of_the_Czech_Republic",
                            "https://en.wikipedia.org/wiki/Senate_of_the_Czech_Republic"),
                   lat = c(14.4104392, 14.3990089, 14.4117831, 16.6021958, 16.6044039, 14.4039458, 14.4053489),
                   lon = c(50.0933681, 50.0895897, 50.0920997, 49.2051925, 49.1977642, 50.0891494, 50.0900269))


# transform the data frame from plain vanilla one to spatial
body <- body %>% 
  sf::st_as_sf(coords = c("lat", "lon"), # columns with geometry
               crs = 4326) # WGS84 is a sensible default...

# prepare a palette - manual colors according to branch column
palPwr <- leaflet::colorFactor(palette = c("executive" = "red", 
                                           "judiciary" = "goldenrod", 
                                           "legislature" = "steelblue"), 
                               domain = body$branch)


# data frame or Prague points
praha <- body %>% 
  dplyr::filter(branch != "judiciary") # a shortcut - all except judiciary are in Prague

# data frame of Brno points
brno <- body %>% 
  filter(branch == "judiciary") # only judiciary happens to be in Brno (and only in Brno)

# first prepare a leaflet plot ...
lplot <- leaflet(data = body) %>% # data = original body - to get the zoom right
  addProviderTiles("CartoDB.Positron") %>% 
  addCircleMarkers(data = praha, # first group
                   radius = 10,
                   fillOpacity = .7,
                   stroke = FALSE,
                   popup = ~htmlEscape(name),
                   color = palPwr(praha$branch), # using already created palette
                   clusterOptions = markerClusterOptions(),
                   group = "Prague") %>% 
  addCircleMarkers(data = brno, # second group
                   radius = 10,
                   fillOpacity = .7,
                   stroke = FALSE,
                   popup = ~htmlEscape(name),
                   color = palPwr(brno$branch), # using already created palette
                   clusterOptions = markerClusterOptions(),
                   group = "Brno") %>% 
  addLegend(position = "bottomright",
            values = ~branch,
            opacity = .7,
            pal = palPwr, # palette declared previously
            title = "Branch") %>% 
  leaflet::addLayersControl(overlayGroups = c("Prague", "Brno"),
                   options = layersControlOptions(collapsed = FALSE)) %>% 
  addResetMapButton()

lplot #  ... then display it

jacmac · May 2, 2020, 2:05pm

Hello and apologies for the delay in responding.

Thank you so much for the guidance and for sharing your example. I'm pleased to say the issue is now resolved and everything works are expected. I can now expand the map to include more historical data and hopefully some machine learning predictions if time permits.

Once again a heartfelt thanks for explaining where I was going wrong so clearly, I really appreciate it!

jlacko · May 2, 2020, 4:02pm

Glad to be of service, and best of luck with your project!

system · May 9, 2020, 4:02pm

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.