Population Density Maps: How to deal with NAs in heightmap matrix for rayshader?

Question

I am replicating a 3D population density map. Everything seems so work fine and when calling the plot_3d() function, my map and its barcharts properly in an rgl window. However, when I run the render_highquality() function, I get

Error in rayvertex::validate_mesh(mesh) : 
  !any(is.na(normals)) is not TRUE

This is due to the fact that my heightmap matrix contains many NAs in the areas where there is no data to be plotted. Replacing all NAs by 0s fixes that error. However, this workaround produces an unwanted visual layer at the base:.

Obviously, I am asking myself how to get rid of this layer or how to deal with the NAs in my heightmap matrix in general. All online sources I've searched through seem to not get this error. For reference, I'm trying to replicate this tutorial and map by justjensen.

As a reproducible example, here is my full code:

remotes::install_github("dmurdoch/rgl", force = T)
remotes::install_github("tylermorganwall/rayshader", force = T)
remotes::install_github("tylermorganwall/rayrender", force = T)

url <- 'https://geodata-eu-central-1-kontur-public.s3.amazonaws.com/kontur_datasets/kontur_population_US_20220630.gpkg.gz'
destination_file <- 'kontur_population_US_20220630.gpkg.gz'
download.file(url, destination_file, 'curl')

library(sf)
library(R.utils)

df_pop_st <- st_read(gunzip(destination_file, remove=FALSE, skip=TRUE))
colnames(df_pop_st) <- tolower(colnames(df_pop_st)) # this is just personal preference

# install.packages("tigris")
# install.packages("tidyverse")
library(tigris)
library(tidyverse)  # need this for filter and ggplot later, but could
                    # just use dplyr for filter and distinct

df_states_st <- states()
colnames(df_states_st) <- tolower(colnames(df_states_st))
statefps <- df_states_st |>
  filter(name %in% c("District of Columbia", "Virginia", "Maryland")) |>
  distinct(statefp)

# Now grab the counties
df_counties_st <- counties()
colnames(df_counties_st) <- tolower(colnames(df_counties_st))
counties_list <- c("Montgomery County", "Alexandria city", "District of Columbia",
                   "Fairfax County", "Loudoun County", "Arlington County",
                   "Prince George's County", "Falls Church city",
                   "Fairfax city"
)

# Pull out DC Area counties and cities and set up the correct CRS for
# mapping later
df_dmv_st <- df_counties_st |>
  filter(statefp %in% statefps$statefp, namelsad %in% counties_list) |>
  filter(statefp != '51' | namelsad != 'Montgomery County') |> # exclude Montgomery County, Virginia
  st_transform(crs=st_crs(df_pop_st))

df_dmv_st %>% 
  ggplot() +
  geom_sf()

# do intersection on data to limit kontur to states
df_pop_dmv_st <- st_intersection(df_pop_st, df_dmv_st)

# define aspect ratio based on bounding box
bb <- st_bbox(df_pop_dmv_st)

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

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

# check by plotting points

df_dmv_st %>% 
  ggplot() +
  geom_sf() +
  geom_sf(data = bottom_left) +
  geom_sf(data = bottom_right, color = "red")

width <- st_distance(bottom_left, bottom_right)

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

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 <- height / width
} else {
  h_ration <- 1
  w_ratio <- width / height
}

library(stars)
size <- 1000

dmv_rast <- stars::st_rasterize(df_pop_dmv_st[,"population", "geom"],
                         nx = floor(size * w_ratio),
                         ny = floor(size * h_ratio))

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

#---------
#here is my workaround
mat[is.na(mat)] <- 0
#---------

# install.packages("RColorBrewer")
# install.packages("colorspace")
library(RColorBrewer)
library(colorspace)
colors = brewer.pal(n=9, name = "PuRd")

texture <- grDevices::colorRampPalette(colors, bias = 3)(256)
swatchplot(texture)

library(rayshader)
library(rayrender)
library(rgl)

rgl::close3d() # Close 

mat %>% 
  height_shade(texture = texture) %>% 
  plot_3d(heightmap = mat,
          zscale = 20,
          solid = F,
          shadowdepth = 0)

render_camera(theta = -15, phi = 50, zoom = .7)
rgl::rglwidget() # Required to show the window in an RStudio Notebook

render_highquality(
  filename = "my_first_plot.png"
)

Answer 1

Here is my improved workaround: Setting NAs to a negative value actually renders the base properly. I'm still not too sure whether there are more elegant versions of tackling this behaviour of rayshader though.

#---------
mat[is.na(mat)] <- -1
#---------

Results in