In our previous post, we introduced the leaflet package along with the GHS-SMOD data source for measuring population density. When we left off, we had loaded the GHS data and prepared it for mapping with leaflet.

This post will use the prepared data (ghs_smod) and the basemap we built in the last post to explore the functionality of leaflet in more depth.

Needless to say, if you haven’t yet taken a look at our previous post, you’ll want to check that out before continuing with this one!

# Load relevant libraries
library(leaflet)
library(terra)
library(dplyr)
library(sf)

Building out our leaflet map

Last time, we showed how to plot raster data on a leaflet map, but we used the default color scheme, which was a little garish.

Fortunately, Leaflet has built-in functions to help you generate your own color palettes for your maps.

Color palettes

Leaflet’s palette functions depend on the type of data your raster contains. Because our data are categorical, we’ll use colorFactor() (categorical levels are often referred to as factors in R).

For numeric data, you might instead consider colorNumeric() or colorBin(). See the leaflet documentation for the available options.

To build a palette (which we call ghs_pal), we provide the set of colors that we want to use along with the domain of values that should be matched to each of these colors. This produces a palette function that takes a raster value as input and produces a corresponding color as output.

In our case, we’ll use an 8-color Inferno palette from the colorspace package. We then match each of these colors to one of the unique values in our raster:

# Get color codes for the Inferno palette from colorspace
colors <- colorspace::sequential_hcl(palette = "inferno", n = 8)

# Map these colors to each unique raster value
ghs_pal <- colorFactor(
  colors,
  c(10, 11, 12, 13, 21, 22, 23, 30),
  na.color = "transparent" # Make missing values transparent
)

Now, we just provide our palette to the colors argument of our raster image layer:

basemap |> 
  addRasterImage(ghs_smod[[1]], colors = ghs_pal)

This is starting to look more consistent with the style of our basemap.

However, based on our color palette, areas of our raster that are coded as water are colored black. This is somewhat confusing, because our basemap uses black to represent land areas. (You can see in the Red Sea that the water transitions from gray to black at our raster boundary.) We can fix this by making the water areas of our raster transparent, so they show the basemap color instead.

We’ll also do the same for the “Very-low” density areas on land, making them show the black basemap color. This will help make the geographic context provided by our basemap more apparent and will help the data blend into the overall map.

Caution

Note that this conversion is intended for visualization purposes. For an analysis, we may not want to convert these values to NA!

We can use terra’s subst() to substitute these categories with missing values:

ghs_smod[[1]] <- subst(ghs_smod[[1]], c("Water", "Very low"), NA)
ghs_smod[[2]] <- subst(ghs_smod[[2]], c("Water", "Very low"), NA)

Finally, let’s mask out the raster data outside of the Ethiopia borders to focus on the data we’re most interested in. We’ve downloaded 2016 Ethiopia boundary files from the DHS Program’s Spatial Data Repository and stored them in our data directory.

First, we’ll load our data with sf:

et_border <- st_read(
  "data/sdr_subnational_boundaries_2024-08-26/shps/sdr_subnational_boundaries.shp",
  quiet = TRUE
)

Then, we’ll dissolve away the internal boundaries to get a single country-level border and project to Web Mercator for consistency with our raster data.

# Dissolve internal boundaries to get external border only
et_border <- st_union(et_border)

# Transform to Web Mercator for consistency with basemap and raster
et_border <- st_transform(et_border, "epsg:3857")

Now, we can mask out the raster values outside the country border with terra’s mask():

# Mask out raster values outside of border
ghs_smod <- mask(ghs_smod, vect(et_border))

basemap |> 
  addRasterImage(ghs_smod[[1]], colors = ghs_pal)

Legends

Our map looks nice, but it’s hard to interpret without a legend. We can add one with addLegend().

To produce a legend, we again need to use our color palette to indicate the colors and labels we want to include in the legend. Since we’ve excluded water and very-low density areas from the map, we don’t want these to show in the legend. We’ll use our palette (ghs_pal) to generate the colors for the other density categories. (Recall that ghs_pal is a function that takes raster values in its domain and outputs the corresponding colors that we set earlier.)

Then, we’ll manually input the desired labels that should correspond to these colors.

legend_colors <- ghs_pal(c(12, 13, 21, 22, 23, 30))
legend_labels <- c("Low", "Rural", "Suburban", "Semi-dense urban", "Dense urban", "Urban center")

Now, we provide this info to addLegend():

basemap |> 
  addRasterImage(ghs_smod[[1]], colors = ghs_pal) |> 
  addLegend(
    title = "GHS-SMOD Population Density",
    position = "bottomright",
    colors = legend_colors,
    labels = legend_labels,
    opacity = 0.7
  )