使用crop和mask后SpatRaster值全为1，地图颜色显示异常

环境信息 #

• 系统：MacOS
• R版本：4.4.0
• 相关包版本：rayshader 0.35.7，terra 1.7-78

输出对比 #

当前输出（颜色单一） #

预期输出（多分类颜色） #

完整原始代码 #

library(terra)
library(giscoR)
library(sf)
library(tidyverse)
library(ggtern)
library(elevatr)
library(png)
library(rayshader)
library(magick)

# 2. COUNTRY BORDERS ----

country_sf <- giscoR::gisco_get_countries(
  country = "BA",
  resolution = "1"
)

plot(sf::st_geometry(country_sf))

png("bih-borders.png")
plot(sf::st_geometry(country_sf))
dev.off()

# 3. DOWNLOAD ESRI LAND COVER TILES ----

urls <- c(
  "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2022/33T_20220101-20230101.tif",
  "https://lulctimeseries.blob.core.windows.net/lulctimeseriesv003/lc2022/34T_20220101-20230101.tif"
)

for(url in urls){
  download.file(
    url = url,
    destfile = basename(url),
    mode = "wb"
  )
}

# 4. LOAD TILES ----

raster_files <- list.files(
  path = getwd(),
  pattern = "20230101.tif$",
  full.names = T
)

crs <- "EPSG:4326"

for(raster in raster_files){
  rasters <- terra::rast(raster)
  
  country <- country_sf |>
    sf::st_transform(
      crs = terra::crs(
        rasters
      )
    )
  
  land_cover <- terra::crop(
    rasters,
    terra::vect(country),
    snap = "in",
    mask = T
  ) |>
    terra::aggregate(
      fact = 5,
      fun = "modal"
    ) |>
    terra::project(crs)
  
  terra::writeRaster(
    land_cover,
    overwrite=TRUE,
    paste0(
      raster,
      "_bosnia",
      ".tif"
    )
  )
}

# 5. LOAD VIRTUAL LAYER ----

r_list <- list.files(
  path = getwd(),
  pattern = "_bosnia",
  full.names = T
)

land_cover_vrt <- terra::vrt(
  r_list,
  "bosnia_land_cover_vrt.vrt",
  overwrite = T
)

# 6. FETCH ORIGINAL COLOURS ----

ras <- terra::rast(
  raster_files[[1]]
)

raster_color_table <- do.call(
  data.frame,
  terra::coltab(ras)
)

head(raster_color_table)

hex_code <- ggtern::rgb2hex(
  r = raster_color_table[,2],
  g = raster_color_table[,3],
  b = raster_color_table[,4]
)

# 7. ASSIGN COLORS TO RASTER----

cols <- hex_code[c(2:3, 5:6, 8:12)]

from <- unique_vals
to <- t(col2rgb(cols))
land_cover_vrt <- na.omit(land_cover_vrt)

land_cover_bosnia <- terra::subst(
  land_cover_vrt,
  from = from,
  to = to,
  names = cols
)

terra::plotRGB(land_cover_bosnia)

# 8. DIGITAL ELEVATION MODEL ----

elev <- elevatr::get_elev_raster(
  locations = country_sf,
  z = 9, clip = "locations"
)

crs_lambert <-
  "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +datum=WGS84 +units=m +no_frfs"

land_cover_bosnia_resampled <- terra::resample(
  x = land_cover_bosnia,
  y = terra::rast(elev),
  method = "near"
) |>
  terra::project(crs_lambert)

terra::plotRGB(land_cover_bosnia_resampled)

img_file <- "land_cover_bosnia.png"

terra::writeRaster(
  land_cover_bosnia_resampled,
  img_file,
  overwrite = T,
  NAflag = 255
)

img <- png::readPNG(img_file)


# 9. RENDER SCENE ----


elev_lambert <- elev |>
  terra::rast() |>
  terra::project(crs_lambert)

elmat <- rayshader::raster_to_matrix(
  elev_lambert
)

h <- nrow(elev_lambert)
w <- ncol(elev_lambert)

elmat |>
  rayshader::height_shade(
    texture = colorRampPalette(
      cols[9]
    )(256)
  ) |>
  rayshader::add_overlay(
    img,
    alphalayer = 1
  ) |>
  rayshader::plot_3d(
    elmat,
    zscale = 12,
    solid = F,
    shadow = T,
    shadow_darkness = 1,
    background = "white",
    windowsize = c(
      w / 5, h / 5
    ),
    zoom = .5,
    phi = 85,
    theta = 0
  )

rayshader::render_camera(
  zoom = .58
)

# 10. RENDER OBJECT ----

u <- "https://dl.polyhaven.org/file/ph-assets/HDRIs/hdr/4k/air_museum_playground_4k.hdr"
hdri_file <- basename(u)

download.file(
  url = u,
  destfile = hdri_file,
  mode = "wb"
)

filename <- "3d_land_cover_bosnia-dark.png"

rayshader::render_highquality(
  filename = filename,
  preview = F,
  light = F,
  environment_light = hdri_file,
  intensity_env = 1,
  rotate_env = 90,
  interactive = F,
  parallel = T,
  width = w * 1.5,
  height = h * 1.5
)

# 11. PUT EVERYTHING TOGETHER ----

c(
  "#419bdf", "#397d49", "#7a87c6", 
  "#e49635", "#c4281b", "#a59b8f", 
  "#a8ebff", "#616161", "#e3e2c3"
)

legend_name <- "land_cover_legend.png"
png(legend_name)
par(family = "mono")

plot(
  NULL,
  xaxt = "n",
  yaxt = "n",
  bty = "n",
  ylab = "",
  xlab = "",
  xlim = 0:1,
  ylim = 0:1,
  xaxs = "i",
  yaxs = "i"
)
legend(
  "center",
  legend = c(
    "Water",
    "Trees",
    "Crops",
    "Built area",
    "Rangeland"
  ),
  pch = 15,
  cex = 2,
  pt.cex = 1,
  bty = "n",
  col = c(cols[c(1:2, 4:5, 9)]),
  fill = c(cols[c(1:2, 4:5, 9)]),
  border = "grey20"
)
dev.off()

# filename <- "land-cover-bih-3d-b.png"

lc_img <- magick::image_read(
  filename
)

my_legend <- magick::image_read(
  legend_name
)

my_legend_scaled <- magick::image_scale(
  magick::image_background(
    my_legend, "none"
  ), 2500
)

p <- magick::image_composite(
  magick::image_scale(
    lc_img, "x7000" 
  ),
  my_legend_scaled,
  gravity = "southwest",
  offset = "+100+0"
)

magick::image_write(
  p, "3d_bosnia_land_cover_final.png"
)

问题核心 #

1. 在第4部分使用terra::crop()时设置了mask=T，该参数会直接将裁剪后的栅格转换为二值掩码栅格（仅保留1和NA），完全丢失原始土地覆盖的分类值，导致后续颜色映射失效。
2. 第7部分的unique_vals变量未定义，会导致subst()操作无法正确匹配分类值与颜色。

修复方案 #

1. 拆分裁剪与掩码操作 #

修改第4部分的循环代码，将crop()和mask()分开执行，保留原始分类值：

# 4. LOAD TILES ----

raster_files <- list.files(
  path = getwd(),
  pattern = "20230101.tif$",
  full.names = T
)

crs <- "EPSG:4326"

for(raster in raster_files){
  rasters <- terra::rast(raster)
  
  country <- country_sf |>
    sf::st_transform(
      crs = terra::crs(rasters)
    )
  
  # 先裁剪到边界范围，再单独应用mask，保留原始分类值
  land_cover <- terra::crop(
    rasters,
    terra::vect(country),
    snap = "in"
  ) |>
    terra::mask(terra::vect(country)) |>  # 单独执行mask
    terra::aggregate(
      fact = 5,
      fun = "modal"
    ) |>
    terra::project(crs)
  
  terra::writeRaster(
    land_cover,
    overwrite=TRUE,
    paste0(raster, "_bosnia.tif")
  )
}

2. 定义unique_vals变量 #

修改第7部分代码，从原始栅格提取对应分类值：

# 7. ASSIGN COLORS TO RASTER----

# 从原始栅格提取唯一分类值
unique_vals <- terra::unique(ras)[[1]]
# 筛选需要映射的分类值（对应cols的索引）
from <- unique_vals[c(2:3,5:6,8:12)]
cols <- hex_code[c(2:3,5:6,8:12)]

to <- t(col2rgb(cols))
land_cover_vrt <- na.omit(land_cover_vrt)

land_cover_bosnia <- terra::subst(
  land_cover_vrt,
  from = from,
  to = to,
  names = cols
)

terra::plotRGB(land_cover_bosnia)

3. 验证修复效果 #

运行修改后的代码后，terra::plotRGB(land_cover_bosnia)应显示多分类颜色的土地覆盖图，后续rayshader渲染会生成符合预期的3D效果。

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