Dynamic mapping

Objectives of the fourth section

This fourth section presents a method for generating dynamic maps for publication on the web. By dynamic, we mean the ability to zoom in and out, and to bring up information when hovering over or clicking on the map. By the end of this sixth section, you will be able to perform the following operations:

• generate a choropleth map in leaflet format
• adjust mouse interactions

Data provided

For this exercise, we will use the following data provided:

• the layer of divisions with population densities (divisions_PK_pop.gpkg see first part)

Interactive choropleth map

We begin by loading the necessary libraries. We will use the sf library rather than terra to load our vector, as leaflet does not support terra’s spatVector format. The RColorBrewer library allows us to use preconfigured colour palettes that are relevant to our project.

library(sf)
library(leaflet)
library(RColorBrewer)
library(dplyr)

We then load the vector layer of divisions containing an attribute relating to population density. We take this opportunity to check the data structure.

# import of the layer of the divisions with densities in sf format
pk_div_pop <- st_read("./Data_bonus/divisions_PK_pop.gpkg")

Reading layer `divisions_PK_pop' from data source 
  `/home/poulpe/Enseignement/hors_UP/LCWU_online/R_for_Geomatics/Data_bonus/divisions_PK_pop.gpkg' 
  using driver `GPKG'
Simple feature collection with 31 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -283875.3 ymin: 2621633 xmax: 1303036 ymax: 4120742
Projected CRS: WGS 84 / UTM zone 42N

# we check the data
print(head(pk_div_pop))

Simple feature collection with 6 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -283875.3 ymin: 2754998 xmax: 1093622 ymax: 3904815
Projected CRS: WGS 84 / UTM zone 42N
        NAME_2  COUNTRY       NAME_1   TYPE_2 Census_2017  area_km2 density_inh
1 Azad Kashmir Pakistan Azad Kashmir Division     4032363  13897.56   290.14899
2        Kalat Pakistan  Balochistan Division     2174722 137881.15    15.77244
3       Makran Pakistan  Balochistan Division     1484788  54564.87    27.21143
4    Nasirabad Pakistan  Balochistan Division     1661077  16799.92    98.87409
5       Quetta Pakistan  Balochistan Division     3764730  63238.25    59.53248
6         Sibi Pakistan  Balochistan Division      963941  26537.85    36.32326
                            geom
1 MULTIPOLYGON (((950408.5 36...
2 MULTIPOLYGON (((235244.5 28...
3 MULTIPOLYGON (((-19025.22 2...
4 MULTIPOLYGON (((457008.1 31...
5 MULTIPOLYGON (((-111280.5 3...
6 MULTIPOLYGON (((430331.3 32...

print(summary(pk_div_pop$density_inh))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  15.77  135.32  310.65  672.74  698.83 4475.84 

The leaflet library only handles layers whose reference coordinate system is WGS84. We therefore change the SCR of our layer if necessary.

# We reproject the layer into WGS84 if necessary (required by Leaflet)
if(st_crs(pk_div_pop)$input != "EPSG:4326") {
  pk_div_pop <- st_transform(pk_div_pop, 4326)
}

We then create a colour palette based on thresholds, referred to here as breaks. We make breaks in 20% quantiles on the population density column. We then apply the colour palette named YlOrRd from RColorBrewer to the population densities according to the thresholds defined above.

# we create a color palette
# we define the breaks from 0 to 100% every 20% (O.2)
breaks <- quantile(pk_div_pop$density_inh, probs = seq(0, 1, 0.2), na.rm = TRUE)

# we ceate the color palette
pal <- colorBin("YlOrRd", domain = pk_div_pop$density_inh, bins = breaks)

We will now create the labels that will be displayed when hovering over a division. We will display the name of the division in bold, followed on the line below by the density rounded to one decimal place.

# we define the labels for the popups
labels <- sprintf(
  "<strong>%s</strong><br/>
  Density : %g inh/km²",
  pk_div_pop$NAME_2,
  round(pk_div_pop$density_inh, 1)
) %>% lapply(htmltools::HTML)

Finally, we create the map itself by adding our layer of divisions on top of an OpenStreetMap background. We adjust various formatting options, add labels, a legend and a zoom level.

# we create the leaflet map
carte <- leaflet(pk_div_pop) %>%
  # we add a basemap
  addProviderTiles(providers$CartoDB.Positron) %>%
  
  # we add the polygons of the divisions
  addPolygons(
    fillColor = ~pal(density_inh),
    weight = 1,
    opacity = 1,
    color = "white",
    dashArray = "3",
    fillOpacity = 0.7,
    highlight = highlightOptions(
      weight = 3,
      color = "#666",
      dashArray = "",
      fillOpacity = 0.8,
      bringToFront = TRUE
    ),
    label = labels,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "15px",
      direction = "auto"
    )
  ) %>%
  
  # we add the legend
  addLegend(
    pal = pal,
    values = ~density_inh,
    opacity = 0.7,
    title = "Population density<br>(inh/km²)",
    position = "bottomright"
  ) %>%
  
  # we center the map on Pakistan
  setView(lng = 69.5, lat = 29.7, zoom = 5)

# we display the map
carte

If we wish to export the HTML file corresponding to this interactive map, simply run the following command once the map has been generated.

htmlwidgets::saveWidget(carte, "./Data_processed/map_divisions_densities.html")

Your turn

Complete the following two tasks:

• Change the various layout settings (colour palette, opacity, zoom level, etc.) to understand them better.
• Repeat the same process for the divisions in a district of your choice.

Adding to external data

Here, we will display a tooltip containing a link to the search results for the division we clicked on, using the Duckduckgo search engine. The script is identical to the previous one, except that a new variable popup is created and called. In this pop-up, we retrieve the name of the division and create the generic search URL with this name.

# loading of the necessary libraries
library(sf)
library(leaflet)
library(RColorBrewer)
library(dplyr)

# we imoort the layer of the divisions with densities in sf format
pk_div_pop <- st_read("./Data_bonus/divisions_PK_pop.gpkg")

Reading layer `divisions_PK_pop' from data source 
  `/home/poulpe/Enseignement/hors_UP/LCWU_online/R_for_Geomatics/Data_bonus/divisions_PK_pop.gpkg' 
  using driver `GPKG'
Simple feature collection with 31 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -283875.3 ymin: 2621633 xmax: 1303036 ymax: 4120742
Projected CRS: WGS 84 / UTM zone 42N

# We reproject the layer into WGS84 if necessary (required by Leaflet)
if(st_crs(pk_div_pop)$input != "EPSG:4326") {
  pk_div_pop <- st_transform(pk_div_pop, 4326)
}

# we create a color palette
# we define the breaks from 0 to 100% every 20% (O.2)
breaks <- quantile(pk_div_pop$density_inh, probs = seq(0, 1, 0.2), na.rm = TRUE)

# we create the color palette
pal <- colorBin("YlOrRd", domain = pk_div_pop$density_inh, bins = breaks)

# we define the labels for the popups
labels <- sprintf(
  "<strong>%s</strong><br/>
  Density : %g inh./km²",
  pk_div_pop$NAME_2,
  round(pk_div_pop$density_inh, 1)
) %>% lapply(htmltools::HTML)

# we create popups with a clickable link
labels_links <- sprintf(
  "<strong>%s</strong><br/>
  Density : %g inh./km²<br/>
  <a href='https://duckduckgo.com/?q=%s&kp=-1&kl=us-en' target='_blank'>
  📊 See on the web</a>",
  pk_div_pop$NAME_2, # we use the name of the division
  round(pk_div_pop$density_inh, 1),
  pk_div_pop$NAME_2
) %>% lapply(htmltools::HTML)

# we create the leaflet map
carte <- leaflet(pk_div_pop) %>%
  # we add the OSM basemap
  addProviderTiles(providers$CartoDB.Positron) %>%
  
  # we add the polygons of the divisions with the relevant symbology
  addPolygons(
    fillColor = ~pal(density_inh),
    weight = 1,
    opacity = 1,
    color = "white",
    dashArray = "3",
    fillOpacity = 0.7,
    highlight = highlightOptions(
      weight = 3,
      color = "#666",
      dashArray = "",
      fillOpacity = 0.8,
      bringToFront = TRUE
    ),
    popup = labels_links,  # we add the clickable links
    popupOptions = popupOptions(
      maxWidth = 250,
      closeButton = TRUE
    ),
    label = labels,
    labelOptions = labelOptions(
      style = list("font-weight" = "normal", padding = "3px 8px"),
      textsize = "13px",
      direction = "auto"
    )
  ) %>%
  
  # we add the legend
  addLegend(
    pal = pal,
    values = ~density_inh,
    opacity = 0.7,
    title = "Population density<br>(inh./km²)",
    position = "bottomright"
  ) %>%
  
  # we center the map on Pakistan
  setView(lng = 69.5, lat = 29.7, zoom = 5)

# we display the map
carte

We now have an interactive map that easily redirects us to an external page.

Your turn

Change the label to display the Goole page for the division you clicked on.