First, we install the packages that we need to use.

library(pacman)
pacman::p_load(leaflet, tidyverse, htmltools, terra, sf, scales, glue, readxl, stringi)

General example 

content <- paste(sep = "<br/>",
  "<b><a href='http://www.samurainoodle.com'>Samurai Noodle</a></b>",
  "606 5th Ave. S",
  "Seattle, WA 98138"
)

leaflet() %>% addTiles() %>%
  addPopups(-122.327298, 47.597131, content,
    options = popupOptions(closeButton = FALSE)
  )

Customized map 

## Create and Customize the Base Map and the Initial View
leaflet() %>%
  setView(lng = 13.3892, lat = 52.5128, zoom = 12) %>% 
  addTiles()

Adding layers 

## Add data layers, Customize
leaflet() %>%
  setView(lng = 13.3892, lat = 52.5128, zoom = 12) %>% 
  addTiles() %>% 
  addCircles(
  lng = 13.3892,   # Longitude of the circle's center
  lat = 52.5128,   # Latitude of the circle's center
  radius = 300,   # Radius of the circle in meters
  color = "blue",   # Color of the circle
)

Adding popups 

# But, we can create popups as well
content <- paste(sep = "<br/>",
                 "<b><a href='https://www.hertie-school.org/en/'>Hertie School</a></b>",
                 "Friedrichstraße 180",
                 "10117 Berlin"
)
  
leaflet() %>% addTiles() %>%
  addPopups(13.389223597194448, 52.51298665365591, content,
            options = popupOptions(closeButton = FALSE)
  )

Practical exercise

Cleaning Data sets

Now, we import the data sets that we are going to use. The deforestation data was collected by the Global Forest Watch and the specific data for Colombia can be found here. The geographic files for the municipalities of Colombia are shared by Laboratorio Urbano - Bogotá and they can be accessed here. Although, the files come in different formats, we decided to used the Shapefile format for this exercise.

deforestation <- read_xlsx("COL.xlsx", 
                           sheet = "Subnational 2 tree cover loss") %>%
  distinct(subnational2, .keep_all = TRUE) %>%
  mutate(subnational2 = toupper(subnational2),
         subnational1 = toupper(subnational1),
         pct_gain = `gain_2000-2020_ha`/area_ha) %>%
  rename(nombre_mpi = subnational2,
         nombre_dpt = subnational1) %>%
  arrange(nombre_mpi)

map <- sf::st_read("shapes/shapes.shp", 
                   stringsAsFactors = FALSE) %>%
  arrange(nombre_mpi)
Reading layer `shapes' from data source 
  `C:\Users\cpedr\OneDrive - Hertie School\GitHub\09-leaflet-satam-jimenez-hossain\shapes\shapes.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 1122 features and 9 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -81.73899 ymin: -4.2473 xmax: -66.874 ymax: 13.39744
Geodetic CRS:  WGS 84

Since we are working with data from Colombia, some of the municipalities for the “deforestation” data set include the accent. Thus, we use the stri_trans_general() function from the package stringi. However, since the Spanish language also uses the Ñ and the Ü, we do not want ot get rid of those. Therefore, we need to create a function:

# Function to remove accents
remove_accents <- function(input_string) {
  paste(sapply(strsplit(input_string, "")[[1]], function(x) ifelse(x %in% c("Ñ", "Ü"), x, stringi::stri_trans_general(x, "Latin-ASCII"))), collapse = "")
}

# Apply the function to the vector
deforestation$nombre_mpi <- sapply(deforestation$nombre_mpi, remove_accents)
deforestation$nombre_dpt <- sapply(deforestation$nombre_dpt, remove_accents)

Merging Data sets

After cleaning the data, we need to merge the different files. There are 91 municipalities with no data available in the deforestation data set

#Merging data
deforestation_map <- full_join(map, deforestation, by = c("nombre_mpi", "nombre_dpt"))
save(deforestation_map, file = "deforestation_map.RData")

#No matches
not_m <- map %>%
  anti_join(deforestation, by = c("nombre_mpi", "nombre_dpt")) %>%
  arrange(nombre_mpi)

Getting creative with the map

Now that we have the data needed for our map, we can create a palette of colors and some tags

#If you have not followed the previous steps, you can import the .RData file
load("deforestation_map.RData")

#The palette of colors is based on a scale of greens that changes based on the gained ha of trees between 2002 and 2020 by municipality
palette <- colorNumeric(palette = "Greens", domain = deforestation_map$`gain_2000-2020_ha`)

#We create some tags by municipality using the glue and the htmltools packages
def_popup <- glue("<strong>{deforestation_map$nombre_mpi}</strong><br />
                    <strong>Department: {deforestation_map$nombre_dpt}</strong><br />
                    Total area (ha): {scales::comma(deforestation_map$area_ha)}<br /> 
                    Gained ha of trees (2002-2020): {scales::comma(deforestation_map$`gain_2000-2020_ha`, accuracy = 1)}<br />
                    Tree cover loss 2022 (ha): {scales::comma(deforestation_map$tc_loss_ha_2022, accuracy = 1)}<br />")  %>%   
  lapply(htmltools::HTML)

#Final map
leaflet() %>%
  #addTiles() %>% #For the default map
  addProviderTiles("CartoDB.Positron") %>% #For one of the options included in the addProviderTiles() function
  #addTiles(urlTemplate = "https://server.arcgisonline.com/ArcGIS/rest/services/Ocean/World_Ocean_Base/MapServer/tile/{z}/{y}/{x}") %>% #For an example of third party tiles
  addPolygons(
    data = deforestation_map,
    fillColor = ~palette(deforestation_map$`gain_2000-2020_ha`),
    label = def_popup,
    fillOpacity = 1,
    color = "grey",
    weight = 1
  )