Making 3D population maps of Argentina

rayshader_argentina.R

# ADAPTED FROM:
# - Spencer Schien's rayshader+kontur tutorial: https://justjensen.co/making-population-density-maps-with-rayrender-in-r/
# - Erik Jensen's adaptation of it: https://justjensen.co/making-population-density-maps-with-rayrender-in-r/ 

# INSTALL NECESSARY PACKAGES

install.packages("sf")
install.packages("rayrender")
install.packages("rayshader")
install.packages("R.utils")
install.packages('tidyverse')
install.packages("stars")
install.packages("RColorBrewer")
install.packages("colorspace")
install.packages("OpenStreetMap")

# LOAD PACKAGES

library(RColorBrewer)
library(colorspace)
library(sf)
library(R.utils)
library(tidyverse)
library(rayrender)
library(rayshader)
library(stars)
library(OpenStreetMap)

# LOAD KONTUR DATA: POPULATION & SHAPEFILES BY
# DOWNLOADING FROM https://data.humdata.org/dataset

# load population data for Argentina 
df_pop_st <- st_read(gunzip("Downloads/kontur_population_AR_20231101.gpkg.gz", remove=FALSE, skip=TRUE)) 

# transform to WGS84
df_pop_transformed <- st_transform(df_pop_st,"WGS84")

#### ASPECT RATIO

# define aspect ratio based on bounding box

bb <- st_bbox(df_pop_transformed)

bottom_left <- st_point(c(bb[["xmin"]], bb[["ymin"]])) |> 
  st_sfc(crs = st_crs(df_pop_transformed))

bottom_right <- st_point(c(bb[["xmax"]], bb[["ymin"]])) |> 
  st_sfc(crs = st_crs(df_pop_transformed))

width <- st_distance(bottom_left, bottom_right)

top_left <- st_point(c(bb[["xmin"]], bb[["ymax"]])) |> 
  st_sfc(crs = st_crs(df_pop_transformed))

height <- st_distance(bottom_left, top_left)

# handle conditions of width or height being the longer side

if (width > height) {
  w_ratio <- 1
  h_ratio <- as.numeric(height / width)
} else {
  h_ratio <- 1
  w_ratio <- as.numeric(width / height)
}

# convert to raster so we can then convert to matrix
# set size to something larger after initial testing, but will increase render time/resources

size <- 1000

rast_arg <- st_rasterize(subset(df_pop_transformed,
                                   df_pop_transformed$population > 10), 
                            nx = floor(size * w_ratio),
                            ny = floor(size * h_ratio))


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

# set color palette
colors = brewer.pal(n=9, name = "PuRd")
texture <- grDevices::colorRampPalette(colors, bias = 3)(256)
swatchplot(texture)

rgl::close3d() # Close any existing renders

# make temp plot

mat_arg |> 
  height_shade(texture = texture) |> 
  plot_3d(heightmap = mat_arg,
          zscale = 20,
          solid = FALSE,
          shadowdepth = 10,
          windowsize = c(500,750),
          background = 'grey80')

# set camera
render_camera(theta = -15, phi = 45, zoom = 0.5) 

# render high-quality output, increase samples for better quality on same x/y resolution image
render_highquality(
  filename = "argentina.png",
  interactive = FALSE,
  light = TRUE,
  width = 2000,
  height = 3000,
  samples = 256
)

### RENDER ONLY PROVINCE, SANTA FE IN THIS CASE

# READ SHAPE FILE
shape <- read_sf("Downloads/arg_adm_unhcr2017_shp/arg_admbnda_adm1_unhcr2017.shp")

# SUBSET TO SANTA FE & TRANSFORM TO WGS84
shape_sf <- subset(shape,shape$ADM1_ES == "Santa Fe")
shape_sf <- st_transform(shape_sf,"WGS84")
df_pop_st <- st_transform(df_pop_st,"WGS84")

# INTERSECT TO LIMIT DATA TO SHAPE
stfe_df <- st_intersection(df_pop_st,shape_sf)

# make sure intersected data remains in wgs84
stfe_df_transformed <- st_transform(stfe_df,"WGS84")

#### ASPECT RATIO
# define aspect ratio based on bounding box
bb <- st_bbox(stfe_df_transformed)
bottom_left <- st_point(c(bb[["xmin"]], bb[["ymin"]])) |> 
  st_sfc(crs = st_crs(stfe_df_transformed))
bottom_right <- st_point(c(bb[["xmax"]], bb[["ymin"]])) |> 
  st_sfc(crs = st_crs(stfe_df_transformed))

width <- st_distance(bottom_left, bottom_right)
top_left <- st_point(c(bb[["xmin"]], bb[["ymax"]])) |> 
  st_sfc(crs = st_crs(stfe_df_transformed))
height <- st_distance(bottom_left, top_left)

# handle conditions of width or height being the longer side
if (width > height) {
  w_ratio <- 1
  h_ratio <- as.numeric(height / width)
} else {
  h_ratio <- 1
  w_ratio <- as.numeric(width / height)
}

# convert to raster so we can then convert to matrix
size <- 1000 # again, increase for real plot at end

rast_stfe <- st_rasterize(subset(stfe_df_transformed,
                                   stfe_df_transformed$population > 10), 
                            nx = floor(size * w_ratio),
                            ny = floor(size * h_ratio))


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


rgl::close3d() # Close exisitng windows

mat_stfe |> 
  height_shade(texture = texture) |> 
  plot_3d(heightmap = mat_stfe,
          zscale = 25,
          solid = FALSE,
          shadowdepth = 10,
          windowsize = c(400,600),
          background = 'grey80')
render_camera(theta = -15, phi = 45, zoom = .5)

# render output
render_highquality(
  filename = "santa_fe.png",
  interactive = FALSE,
  light = TRUE,
  width = 2000,
  height = 3000,
  samples = 256
)

### MAP OF DEPARTMENT (GENERAL LOPEZ)

# load pop data again
df_pop_st <- st_read(gunzip("Downloads/kontur_population_AR_20231101.gpkg.gz", remove=FALSE, skip=TRUE))
# load shape
shape <- read_sf("Downloads/arg_adm_unhcr2017_shp/arg_admbnda_adm2_unhcr2017.shp")
shape_gl <- subset(shape,shape$ADM2_REF == "General Lopez")
shape_gl <- st_transform(shape_gl,"WGS84")

# combine
df_pop_st <- st_transform(df_pop_st,"WGS84")
gl_dt <- st_intersection(df_pop_st,shape_gl)
gl_dt_transformed <- st_transform(gl_dt,"WGS84")

#### ASPECT RATIO

# define aspect ratio based on bounding box

bb <- st_bbox(gl_dt_transformed)

bottom_left <- st_point(c(bb[["xmin"]], bb[["ymin"]])) |> 
  st_sfc(crs = st_crs(gl_dt_transformed))

bottom_right <- st_point(c(bb[["xmax"]], bb[["ymin"]])) |> 
  st_sfc(crs = st_crs(gl_dt_transformed))

width <- st_distance(bottom_left, bottom_right)
top_left <- st_point(c(bb[["xmin"]], bb[["ymax"]])) |> 
  st_sfc(crs = st_crs(gl_dt_transformed))

height <- st_distance(bottom_left, top_left)

if (width > height) {
  w_ratio <- 1
  h_ratio <- as.numeric(height / width)
} else {
  h_ratio <- 1
  w_ratio <- as.numeric(width / height)
}

# convert to raster 4 matrix

size <- 1000 # don't forget to set higher

rast_gl <- st_rasterize(subset(gl_dt_transformed,
                                   gl_dt_transformed$population > 10), 
                            nx = floor(size * w_ratio),
                            ny = floor(size * h_ratio))


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

rgl::close3d() 

mat_gl |> 
  height_shade(texture = texture) |> 
  plot_3d(heightmap = mat_gl,
          zscale = 25,
          solid = FALSE,
          shadowdepth = 10,
          windowsize = c(1000,700),
          background = 'grey80')

render_camera(theta = -15, phi = 40, zoom = .5)

render_highquality(
  filename = "test_plot2.png",
  interactive = FALSE,
  light = TRUE,
  width = 3000,
  height = 2100,
  samples = 256
)