人口密度地图：如何处理rayshader高度矩阵中的NA值？

问题描述 #

复刻3D人口密度地图时，plot_3d()能在rgl窗口正常显示地图和柱状图，但执行render_highquality()时触发错误：

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

错误源于高度矩阵（heightmap matrix）中存在大量无数据区域的NA值。临时将所有NA替换为0可解决报错，但会在地图底部生成多余的平坦图层。

可复现代码 #

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)) # 个人习惯统一列名

# install.packages("tigris")
# install.packages("tidyverse")
library(tigris)
library(tidyverse)  

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)

# 获取目标县/市
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"
)

# 提取DMV区域并转换CRS
df_dmv_st <- df_counties_st |>
  filter(statefp %in% statefps$statefp, namelsad %in% counties_list) |>
  filter(statefp != '51' | namelsad != 'Montgomery County') |> # 排除弗吉尼亚州的Montgomery County
  st_transform(crs=st_crs(df_pop_st))

df_dmv_st %>% 
  ggplot() +
  geom_sf()

# 裁剪人口数据到目标区域
df_pop_dmv_st <- st_intersection(df_pop_st, df_dmv_st)

# 计算边界框与宽高比
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))

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)

# 设置宽高比
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))

#---------
# 临时 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() # 关闭已有窗口

# 绘制3D地图
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() # 在RStudio Notebook中显示窗口

# 高质量渲染报错
render_highquality(
  filename = "my_first_plot.png"
)

解决方案 #

方案1：保留NA值，通过渲染参数隐藏底部图层 #

不需要替换NA为0，直接调整plot_3d和render_highquality的参数：

1. 删除mat[is.na(mat)] <- 0这行代码，保留原始高度矩阵的NA值。
2. 修改plot_3d的参数，启用底部平面并将其置于所有柱状图下方，同时匹配背景色：

mat %>% 
  height_shade(texture = texture) %>% 
  plot_3d(heightmap = mat,
          zscale = 20,
          solid = TRUE,
          soliddepth = -max(mat, na.rm = TRUE)/20 - 1, # 确保底部平面低于所有数据柱
          solidcolor = "white", # 与背景色一致
          shadowdepth = 0)

3. 调用render_highquality时指定背景色：

render_highquality(
  filename = "my_first_plot.png",
  background = "white"
)

此时NA区域会被背景色填充，不会出现多余的底部图层。

方案2：手动处理网格，移除NA对应的面 #

如果保留NA值后仍报错，可直接处理rayshader生成的3D网格：

1. 删除mat[is.na(mat)] <- 0，保留原始NA。
2. 在plot_3d之后添加以下代码，移除包含NA的网格面：

# 获取当前场景
scene <- rayvertex::get_scene()
# 移除NA对应的网格面
scene$objects[[1]]$mesh <- rayvertex::remove_na_faces(scene$objects[[1]]$mesh)
# 更新场景
rayvertex::set_scene(scene)

3. 再执行render_highquality即可正常渲染，无多余图层。

方案3：优化栅格化步骤，减少NA区域 #

在栅格化时限定范围，减少不必要的NA：
修改st_rasterize的调用，添加extent参数限定为目标区域的边界框：

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

这样生成的高度矩阵NA区域更少，再结合方案1或2处理即可。

内容的提问来源于stack exchange，提问作者Silvan