Block 4: Interactive Maps
Max Hofmann
Welcome to Block 4!
Here, we will build an interactive map with the package leaflet.
Please read the following code instructions carefully, try to understand the code that follows each instruction, execute it and see what happens. Do not hesitate to change the code or write your own code, and execute it as well!
1. Preparations & data inspection
As usual, we start by loading the packages and data we need, and by defining our working directory.
#install.packages("leaflet")
library(dplyr)##
## Attache Paket: 'dplyr'## Die folgenden Objekte sind maskiert von 'package:stats':
##
## filter, lag## Die folgenden Objekte sind maskiert von 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(leaflet)
library(sf)## Warning: Paket 'sf' wurde unter R Version 4.4.2 erstellt## Linking to GEOS 3.12.2, GDAL 3.9.3, PROJ 9.4.1; sf_use_s2() is TRUElibrary(terra)## terra 1.7.78library(viridis)## Lade nötiges Paket: viridisLitedir <- "C:/Users/hofmann/Desktop/IAMO_B/2024_GEO-WB6/WP2 - CAPACITY BUILDING/WP 2.1 - Training Courses/Course 04 - Visualizing and Analyzing Spatial Data in R/__SHARED_Participants/Block 4 - Interactive Maps/country_data/ALB/"
borders <- terra::vect(paste0(dir, "GADM_Administrative_Areas/gadm41_ALB_1.shp"))Again, let us clean up the attribute table a bit and only keep the column “NAME_1”, and then calculate the area in a new column with the expanse()-function:
values(borders) <- values(borders)[,4]
borders$area <- expanse(borders)
borders$area <- borders$area/1000000 ## convert from square meters to square kilometersWe could now plot the region borders in a static manner as we have done several times before:
terra::plot(borders, "area", col=viridis(5))
2. Leaflet - Step by step
However, we will use a different map plotting package now, which allows zooming, highlighting polygons, etc. - it is called leaflet. There is a very good overview about the leaflet package here: https://rstudio.github.io/leaflet/
For our “borders” shapefile to be rendered by leaflet, we have to project it to WGS84, and then convert it to an sf object:
borders <- project(borders, "EPSG:4326")
borders_sf <- st_as_sf(borders)Please note that the first column in the attribute table has now been renamed to “value” by default:
borders_sf| value | area | geometry |
|---|---|---|
| Berat | 1816.1049 | POLYGON ((20.3705 40.41808,… |
| Dibër | 2613.1600 | POLYGON ((20.30172 41.3728,… |
| Durrës | 779.7658 | POLYGON ((19.62202 41.40051… |
| Elbasan | 3167.3386 | POLYGON ((20.06979 40.81438… |
| Fier | 1875.7068 | POLYGON ((19.79524 40.42089… |
| Gjirokastër | 2896.8124 | POLYGON ((20.33032 39.91824… |
| Korçë | 3720.1986 | POLYGON ((20.58344 40.10232… |
| Kukës | 2408.3087 | POLYGON ((20.51898 41.84017… |
| Lezhë | 1665.0100 | MULTIPOLYGON (((19.56031 41… |
| Shkodër | 3504.0240 | POLYGON ((19.50945 41.85914… |
| Tiranë | 1636.3544 | POLYGON ((19.49299 41.03647… |
| Vlorë | 2644.8610 | MULTIPOLYGON (((20.17072 39… |
Similarly to ggplot(), we sequentially build a leaflet() map by adding layers to each other. However, instead of using the sign “+”, we here have to connect the layers with the operator “%>%”. A simple leaflet map would just consist of some Open Street Map tiles that you can add with the function addTiles(). For other tiles, you have to use the function addProviderTiles(). You can find previews of many available tiles here: https://leaflet-extras.github.io/leaflet-providers/preview/
leaflet() %>%
addTiles()leaflet() %>%
addProviderTiles(providers$CartoDB.Positron)You have seen that by default, we see the entire world. We can use the function setView() to zoom in and set the center coordinate of our window to Albania - in that way, we decide which part of the world should be seen when the map is loaded:
leaflet() %>%
setView(20,41, zoom = 7) %>%
addTiles()We can also add several tile sets and use addLayersControl() to create a little window in the upper right corner to select one of them, and add a mini world map in the lower right corner that helps us orientate ourselves on the map:
leaflet() %>%
setView(20,41, zoom = 7) %>%
addTiles(group = "OSM") %>%
addProviderTiles(providers$CartoDB.Positron, group = "Positron") %>%
addLayersControl(baseGroups = c("OSM", "Positron"),
options = layersControlOptions(collapsed = FALSE)) %>%
addMiniMap()Let us now add the shapefile “borders” to the map. We use addPolygons() for this. Whilst the argument “fillColor” determines the color of the area of the polygons, the argument “color” determines the color of the polygon borders. “Weight” determines the boldness of these borders. The argument “fillOpacity” determines the level of transparency of the polygon color:
leaflet() %>%
setView(20,41, zoom = 7) %>%
addTiles() %>%
addPolygons(data = borders_sf,
fillColor = "red",
weight = 2,
color = "black",
fillOpacity = 0.6)Since leaflet maps are interactive, we can define actions that happen when the mouse cursor hovers over a feature in the map. Let us first define some parameters to highlight polygons. The function highlightOptions() is inside addPolygons() but redefines the arguments defined therein for when the mouse cursor passes over a polygon:
leaflet() %>%
setView(20,41, zoom = 7) %>%
addTiles() %>%
addPolygons(data = borders_sf,
fillColor = "red",
weight = 1,
color = "black",
fillOpacity = 0.6,
highlight = highlightOptions(weight = 3, color = "white", fillOpacity = 0.9, bringToFront = FALSE))Let us now add labels to the polygons. We would like to show the name of each region when we pass over it - for this, we can use the column “value” in “borders_sf”. The argument “label” would be sufficient to achieve that, but we can further format the labels with “labelOptions”
leaflet() %>%
setView(20,41, zoom = 7) %>%
addTiles() %>%
addPolygons(data = borders_sf,
fillColor = "red",
weight = 1,
color = "black",
fillOpacity = 0.6,
highlight = highlightOptions(weight = 3, color = "white", fillOpacity = 0.9, bringToFront = FALSE),
label = borders_sf$value,
labelOptions = labelOptions(style = list("font-weight" = "bold",
padding = "3px 8px"),
textsize = "15px",
direction = "auto"))Coloring the polygons according to the values in the column “area” in the attribute table is a bit more complicated in leaflet compared to terra::plot. To do so, we first have to define the bins, i.e. the intervals that define which values should be displayed with the same color. Let us define 5 equal intervals from the minimum to the maximum value of the data range in the column “area”. The resulting object “bins” is a vector of six elements that determine the interval limits. We also then have to define a color palette function “pal” that is defined using “bins”.
bins <- seq(from = min(na.omit(borders_sf$area)),
to = max(na.omit(borders_sf$area))+0.001,
by = (max(na.omit(borders_sf$area))+0.001 - min(na.omit(borders_sf$area)))/5)
bins ## [1] 779.7658 1367.8526 1955.9393 2544.0261 3132.1128 3720.1996pal <- colorBin(magma(5), domain = na.omit(borders_sf$area), bins = bins)We now just have to supply the data column of the shapefile to “pal”, and supply this to the argument “fillColor” in addPolygons(). For each polygon, the palette function will allocate the data value of the polygon to one of the 5 intervals defined by “bins”, and then color the polygon according to color ramp defined in “pal”:
leaflet() %>%
setView(20,41, zoom = 7) %>%
addTiles() %>%
addPolygons(data = borders_sf,
fillColor = ~pal(borders_sf$area),
weight = 1,
color = "black",
fillOpacity = 0.6,
highlight = highlightOptions(weight = 3, color = "white", fillOpacity = 0.9, bringToFront = FALSE),
label = borders_sf$value,
labelOptions = labelOptions(style = list("font-weight" = "bold",
padding = "3px 8px"),
textsize = "15px",
direction = "auto"))We can also expand the labels to include the respective value in a second row, in the following way:
labels <- sprintf(paste0("<strong>", borders_sf$value, "</strong><br/>Area: %s ", "%s"), round(borders_sf$area, digits=2), "km²") %>% lapply(htmltools::HTML)
leaflet() %>%
setView(20,41, zoom = 7) %>%
addTiles() %>%
addPolygons(data = borders_sf,
fillColor = ~pal(borders_sf$area),
weight = 1,
color = "black",
fillOpacity = 0.6,
highlight = highlightOptions(weight = 3, color = "white", fillOpacity = 0.9, bringToFront = FALSE),
label = labels,
labelOptions = labelOptions(style = list("font-weight" = "normal",
padding = "3px 8px"),
textsize = "15px",
direction = "auto")) Let us add a legend to the map, and store the map as an object:
mymap <- leaflet() %>%
setView(20,41, zoom = 7) %>%
addTiles() %>%
addPolygons(data = borders_sf,
fillColor = ~pal(borders_sf$area),
weight = 1,
color = "black",
fillOpacity = 0.6,
highlight = highlightOptions(weight = 3, color = "white", fillOpacity = 0.9, bringToFront = FALSE),
label = labels,
labelOptions = labelOptions(style = list("font-weight" = "normal",
padding = "3px 8px"),
textsize = "15px",
direction = "auto")) %>%
addLegend(data = borders_sf,
pal = pal,
values = borders_sf$area,
opacity = 0.9,
title = "Size in square kilometers",
position = "bottomleft")
mymapLet us add the river networks and the hospitals. First, we have to pre-process the data:
rivers <- terra::vect(paste0(dir, "AQUASTAT_Rivers/Rivers.shp"))
rivers <- project(rivers, "EPSG:4326")
rivers_sf <- st_as_sf(rivers)
points <- terra::vect(paste0(dir, "OSM_Open_Street_Maps/Points_Interest.shp"))
hospitals <- points[points$fclass == "hospital",]
hospitals <- project(hospitals, "EPSG:4326")
hospitals_sf <- st_as_sf(hospitals)
mymap %>%
addCircleMarkers(data = hospitals_sf, radius=1, color='red') %>%
addPolylines(data = rivers, color = "skyblue", opacity = 1, weight = 1)3. Final Map
Let us add from above the tile selection and mini map:
myfinalmap <- mymap %>%
addCircleMarkers(data = hospitals_sf, radius=1, color='red') %>%
addPolylines(data = rivers, color = "skyblue", opacity = 1, weight = 1) %>%
addTiles(group = "OSM") %>%
addProviderTiles(providers$CartoDB.Positron, group = "Positron") %>%
addLayersControl(baseGroups = c("OSM", "Positron"),
options = layersControlOptions(collapsed = FALSE)) %>%
addMiniMap()
myfinalmapNow, please upload this script to RPubs!
Congratulations, you made it to the end of the script! You can now move on to the assignment. Good Luck!