How to Create a 3D Population Density Map in R

In this post, we’ll learn how to create a 3D population density map using R. Let’s learn to create engaging 3D population density maps in R, focusing on Bangladesh’s demographic data. This tutorial is ideal for data scientists and analysts, guiding you through advanced visualization techniques using R packages. Elevate your data presentation and analysis skills with this concise guide. Let’s dive in!

A Data Professional Guide to Advanced Visualization Techniques.

Population Density Map — Bangladesh

Step 1: Install Required Packages
First, we need to install the necessary packages. Run the following commands in your R console. Try to install them one by one, it might require restarting the R-session several times.

install.packages("sf", dependencies=TRUE)
install.packages("tmap", dependencies=TRUE)
install.packages("mapview", dependencies=TRUE)
install.packages("stars", dependencies=TRUE)
install.packages("rayshader", dependencies=TRUE)
install.packages("MetBrewer", dependencies=TRUE)
install.packages("rayrender")
install.packages("extrafont", dependencies=TRUE)
install.packages("magick", dependencies=TRUE)

Step 2: Load Packages and Set Options
Load the required libraries and set the RGL options:

options(rgl.useNULL = FALSE)

require(tidyverse)
require(sf)
require(tmap)
require(ggplot2)
require(mapview)
require(stars)
require(rayshader)
require(MetBrewer)
require(colorspace)
require(rayrender)
require(magick)
require(extrafont)

Step 3: Load and Transform Data
We’ll load the population data and administrative boundaries, transforming them into a suitable coordinate system. The data is downloaded from Kontur Population.

bd_hex <-
  st_read("Data/Kontur_Population_20231101.gpkg") %>%
  st_transform(3106)

bd_admin <-
  st_read("Data/Bangladesh_boundaries_20230628.gpkg") %>%
  st_transform(3106)

Step 4: Check and Create Boundaries
Inspect the ‘name_en’ column and create the boundary for Bangladesh.

• Use the filter option to plot specific districts and divisions on the map.
• Only keep the dhaka_hexpart if you want to plot the Dhaka Division only and plot only the dhaka_hex instead of bd_hex

distinct_names <- bd_admin %>%
  distinct(name_en)
print(distinct_names)


# Creating BD Boundary
bd_boundary <-
  bd_admin %>%
  # filter(name_en == 'Dhaka Division') %>% # Filtering Dhaka Only
  st_geometry %>%
  st_union %>%
  st_sf %>%
  st_make_valid()

# Only keep the part below if you want to plot the Dhaka Division only and plot only the dhaka_hex instead of bd_hex
# Otherwise comment it.

bd_boundary %>% 
  ggplot()+
  geom_sf()

dhaka_hex <- st_intersection(bd_hex, bd_boundary) %>% 
  st_transform(crs = 3106)

Step 5: Plot Boundaries for Verification
Visualize the hex data and boundaries to ensure accuracy.

# check the boundary plot
ggplot(bd_hex) +
  geom_sf(aes(fill = population),
          color = "gray66",
          linewidth = 0) +
  geom_sf(
    data = bd_boundary,
    fill = NA,
    color = "black",
    linetype = "dashed",
    linewidth = 1
  )

Step 6: Calculate Aspect Ratio
Determine the aspect ratio for the map based on the bounding box of the boundary.

# setting the bd boundary as a bounding box
bbox <- st_bbox(bd_boundary)

# finding the aspect ratio
bottom_left <- st_point(c(bbox[["xmin"]], bbox[["ymin"]])) %>%
  st_sfc(crs = 3106)
bottom_right <- st_point(c(bbox[["xmax"]], bbox[["ymin"]])) %>%
  st_sfc(crs = 3106)
top_left <- st_point(c(bbox[["xmin"]], bbox[["ymax"]])) %>%
  st_sfc(crs = 3106)
top_right <- st_point(c(bbox[["xmin"]], bbox[["ymax"]])) %>%
  st_sfc(crs = 3106)



width <- st_distance(bottom_left, bottom_right)
height <- st_distance(bottom_left, top_left)

if(width > height) {
  w_ratio = 1
  h_ratio = height / width
  
} else {
  h_ratio = 1.1
  w_ratio = width / height
}

Step 7: Rasterize Population Data
Convert the population data into a raster format suitable for 3D rendering.

• For interactively checking the 3D plot setting the size low will help render in real time.
• To improve the quality of the 3D image when saving, change the settings to a higher resolution.

# convert to raster to convert to matrix
# size = 100
size = 1000 * 3.5

pop_raster <- st_rasterize(
  bd_hex,
  nx = floor(size * w_ratio) %>% as.numeric(),
  ny = floor(size * h_ratio) %>% as.numeric()
)

pop_matrix <- matrix(pop_raster$population,
                     nrow = floor(size * w_ratio),
                     ncol = floor(size * h_ratio))

Step 8: Define Color Palette
Select a color palette from the MetBrewer library and customize it for your map.

color <- MetBrewer::met.brewer(name="Benedictus", direction = -1)

tx <- grDevices::colorRampPalette(subset_colors, bias = 4.5)(256)
swatchplot(tx)
swatchplot(subset_colors)

Step 9: Render 3D Map
Use Rayshader to create a 3D representation of the population density.

# Close any existing 3D plot before plotting another
rgl::close3d()

pop_matrix %>%
  height_shade(texture = tx) %>%
  plot_3d(heightmap = pop_matrix,
          zscale = 250 / 4.5,
          solid = F,
          shadowdepth = 0)

# Adjusting Camera Angle
render_camera(theta = 0,
              phi = 70,
              zoom = 0.55,
              fov = 100
)

# To interactively view the 3D plot
rgl::rglwidget()

Step 10: Render in high-quality and Save Image
Fine-tune the camera angle and render a high-quality image of the 3D map.

outfile <- glue::glue("Plots/Dhaka_Benedictus_4.png")

{
  start_time <- Sys.time()
  cat(crayon::cyan(start_time), "\n")
  if(!file.exists(outfile)) {
    png::writePNG(matrix(1), target = outfile)
  }
  
  render_highquality(
    filename = outfile,
    interactive = F,
    lightdirection = 55, #Degree
    lightaltitude = c(30, 80),
    lightcolor = c("white", "white"),  # Set both lights to white
    lightintensity = c(600, 100),
    width = 1400,
    height = 1580,
    samples = 500
    #samples = 2
  )
  
  end_time <- Sys.time()
  diff <- end_time - start_time
  cat(crayon::cyan(diff), "\n")
}

Population Density Map — Bangladesh

Step 11: Annotate the image

pop_raster <- image_read("final_plot_bd_Benedictus_3.png")

text_color <- darken(subset_colors[3], .4)
swatchplot(text_color)


# Automatically enable font support
showtext_auto()

# Download and register the Philosopher font from Google Fonts
font_add_google("Philosopher", regular = "400", bold = "700")

pop_raster %>%
  image_annotate("Bangladesh",
                 gravity = "northeast",
                 location = "+50+50",
                 color = text_color,
                 size = 120,
                 font = "Philosopher",
                 weight = 700,
                 # degrees = 0,
  ) %>%
  image_annotate("POPULATION DENSITY MAP",
                 gravity = "northeast",
                 location = "+50+175",
                 color = text_color,
                 size = 28.5,
                 font = "Philosopher",  # Corrected font name
                 weight = 500,
                 # degrees = 0,
  ) %>%
  image_annotate("Visualization by: Niloy Biswas with Rayshader\nData: Kontur Population 2023",
                 gravity = "southwest",
                 location = "+20+20",
                 color = alpha(text_color, .8),
                 font = "Philosopher",  # Corrected font name
                 size = 25,
                 # degrees = 0,
  ) %>%
  image_write("Annotated_plot_bd_Benedictus_3.png", format = "png", quality = 100)

Annotated Population Density Map — Bangladesh

After filtering each division and adjusting the color palette, you will obtain images like below.

Dhaka Division

Population Density Map — Dhaka

Chattogram Division

Population Density Map — Chattogram

Barishal Division

Population Density Map — Barishal

Conclusion

Congratulations! You’ve successfully created a 3D population density map. This map can provide valuable insights into population distribution and can be a powerful tool for presentations and analysis.

• For detailed code, explore the GitHub Repository.
• To delve into more projects, visit my portfolio.
• For any inquiries, feel free to reach out to me on LinkedIn.