合并多个RasterBrick时出现'data'类型错误的技术求助

问题背景 #

我有6个基于野火排放数据的RasterBrick对象，想要合并为一个完整的RasterBrick。使用stack()堆叠所有对象后调用brick()时，出现如下错误：

Error in matrix(unlist(ini), ncol = 2, byrow = TRUE) :  'data' must be of a vector type, was 'NULL' In addition: There were 50 or more warnings (use warnings() to see the first 50)

仅堆叠前两个brick时brick()能正常运行，但加入第三个就触发错误，推测部分RasterBrick的格式或属性存在问题。

以下是创建和堆叠这些RasterBrick的代码：

lng <- seq(from=-130, to=-110, by=res_xy)
lat <- seq(from=30, to=45, by=res_xy)

lng2 <- rep.int(lng, length(lat))
lat2 <- rep.int(lat, length(lng))
lat2 <- sort(lat2)
grd <- as.data.frame(cbind(lng2, lat2))
grd$PM25 <- 0 #no emissions by default

raster_cat <- list()
rdsChunk <- 60 #roughly two months

start_time <- Sys.time()
for (i in 1:365) {
  print(paste('Processing days (total 365):', i))
  pm25_tmp <- subset(finn24_2018_pm25_cal2, DAY == i)[, c("LONGI", "LATI", "PM25")]
  
  if (nrow(pm25_tmp) == 0){#no emissions data for that day
    r <- rasterFromXYZ(grd, res = c(res_xy, res_xy), crs = "+proj=longlat +datum=WGS84")    
  } else{
    # for each point in pm25_tmp, find the nearest neighbor from grd
    closest <- RANN::nn2(data = grd[, 1:2], query = pm25_tmp[, 1:2], k = 1)
    
    # Get coordinates of nearest neighbor
    pm25_tmp$lng2 <- grd[closest$nn.idx, "lng2"]
    pm25_tmp$lat2 <- grd[closest$nn.idx, "lat2"]
    
    # it might happen that multiple pm25_tmp records have a single grd record match
    pm25 <- ddply(pm25_tmp, .(lng2, lat2), summarise,
                  PM25=mean(PM25, na.rm=TRUE))
    grd_sub1 <- merge(pm25, grd[, 1:2], by=c("lng2", "lat2"), all.x=TRUE, all.y=FALSE)
    # use anti_join to find those ones that are in grd but not grid_sub1
    grd_sub2 <- anti_join(grd[,1:2], grd_sub1[, 1:2], by = NULL) # NULL indicates using all the columns
    grd_sub2$PM25 <- 0
    grd2 <- rbind(grd_sub1, grd_sub2)
    
    r <- rasterFromXYZ(grd2, res = c(res_xy, res_xy), crs = "+proj=longlat +datum=WGS84")
  } 
  
  raster_cat[[i]] <- r
  if(i==365){ # the last chunk of data with total records more than rdsChunk 
    file_name <- paste("PM25_2018Raster_Cal", i, ".grd", sep = "")
    raster_stack <- stack(lapply(((i%/%rdsChunk - 1) * rdsChunk + 1):i, function(i)  raster_cat[[i]])) #301-365, %/% and %% for quotient and remainder
    writeRaster(raster_stack, file = file.path(stilt_wd, "calsacto/data", file_name), format="raster", overwrite=TRUE)
  } else if ((i%%rdsChunk == 0) & ((365 - i )> rdsChunk)) { #i.e., not the last chunk yet
    file_name <- paste("PM25_2018Raster_Cal", i, ".grd", sep = "") #default raster format is grd
    raster_stack <- stack(lapply((i-rdsChunk+1):i, function(i)  raster_cat[[i]]))
    writeRaster(raster_stack, file = file.path(stilt_wd, "calsacto/data", file_name), format="raster", overwrite=TRUE)
    raster_cat <- list() #reset raster list
  }
}
end_time <- Sys.time()
print(paste("Total time spent: ", end_time-start_time, sep = ""))

#read one raster stack and convert it to raster brick
stack_1 <- file.path(stilt_wd, "calsacto/data", "PM25_2018Raster_Cal60.grd")
brick_1 <- brick(stack_1)
nlyrs_1 <- nlayers(brick_1)
nlyrs_1

stack_2 <- file.path(stilt_wd, "calsacto/data", "PM25_2018Raster_Cal120.grd")
brick_2 <- brick(stack_2)
nlyrs_2 <- nlayers(brick_2)
nlyrs_2

stack_3 <- file.path(stilt_wd, "calsacto/data", "PM25_2018Raster_Cal180.grd")
brick_3 <- brick(stack_3)
nlyrs_3 <- nlayers(brick_3)
nlyrs_3

stack_4 <- file.path(stilt_wd, "calsacto/data", "PM25_2018Raster_Cal240.grd")
brick_4 <- brick(stack_4)
nlyrs_4 <- nlayers(brick_4)
nlyrs_4

stack_5 <- file.path(stilt_wd, "calsacto/data", "PM25_2018Raster_Cal300.grd")
brick_5 <- brick(stack_5)
nlyrs_5 <- nlayers(brick_5)
nlyrs_5

stack_6 <- file.path(stilt_wd, "calsacto/data", "PM25_2018Raster_Cal365.grd")
brick_6 <- brick(stack_6)
nlyrs_6 <- nlayers(brick_6)
nlyrs_6

nlyrs <- nlyrs_1 + nlyrs_2 + nlyrs_3 + nlyrs_4 + nlyrs_5 + nlyrs_6
nlyrs

emi_stack_2018 <- stack(brick_1, brick_2, brick_3, brick_4, brick_5, brick_6)

emi_brick_2018 <- brick(emi_stack_2018)

问题根源与修复方案 #

1. 先验证所有brick的空间属性一致性 #

报错的核心原因是部分RasterBrick的空间属性不匹配，比如范围、分辨率或CRS不一致。先运行以下代码检查：

# 定义检查函数
check_brick_attr <- function(brick_obj) {
  list(
    extent = as.character(extent(brick_obj)),
    resolution = paste(res(brick_obj), collapse = ","),
    crs = as.character(crs(brick_obj)),
    layer_count = nlayers(brick_obj)
  )
}

# 批量检查所有brick
lapply(list(brick_1, brick_2, brick_3, brick_4, brick_5, brick_6), check_brick_attr)

如果输出结果中有不一致的属性，说明对应chunk的栅格生成有问题，需要重新生成该部分数据。

2. 修复chunk生成的逻辑漏洞 #

原代码中最后一个chunk的写入逻辑没问题，但建议直接用brick()生成并写入，避免stack()转brick()的潜在问题：

# 修改循环中写入栅格的代码
if(i==365){ 
  file_name <- paste("PM25_2018Raster_Cal", i, ".grd", sep = "")
  # 直接生成brick而非stack
  raster_brick <- brick(lapply(((i%/%rdsChunk - 1) * rdsChunk + 1):i, function(j)  raster_cat[[j]]))
  writeRaster(raster_brick, file = file.path(stilt_wd, "calsacto/data", file_name), format="raster", overwrite=TRUE)
} else if ((i%%rdsChunk == 0) & ((365 - i )> rdsChunk)) { 
  file_name <- paste("PM25_2018Raster_Cal", i, ".grd", sep = "")
  raster_brick <- brick(lapply((i-rdsChunk+1):i, function(j)  raster_cat[[j]]))
  writeRaster(raster_brick, file = file.path(stilt_wd, "calsacto/data", file_name), format="raster", overwrite=TRUE)
  raster_cat <- list() #reset raster list
}

3. 简化合并流程，跳过多余的stack步骤 #

不需要先stack()再转brick()，直接用brick()合并多个RasterBrick对象：

# 直接合并所有brick
emi_brick_2018 <- brick(brick_1, brick_2, brick_3, brick_4, brick_5, brick_6)

如果仍然报错，可逐步合并排查问题brick：

# 逐步验证，定位问题brick
test_merge <- brick(brick_1, brick_2)
test_merge <- brick(test_merge, brick_3) # 此处报错则brick_3有问题
test_merge <- brick(test_merge, brick_4)
# ...以此类推

4. 提前验证空数据日的栅格模板 #

当某天无排放数据时，grd的坐标有效性直接影响栅格生成，提前检查：

# 检查grd是否存在重复坐标
anyDuplicated(grd[,1:2]) # 正常应返回0

最终优化后的合并代码 #

# 先验证所有brick的空间属性完全一致
all_extents <- sapply(list(brick_1, brick_2, brick_3, brick_4, brick_5, brick_6), function(x) as.character(extent(x)))
all_resolutions <- sapply(list(brick_1, brick_2, brick_3, brick_4, brick_5, brick_6), function(x) paste(res(x), collapse=","))
all_crs <- sapply(list(brick_1, brick_2, brick_3, brick_4, brick_5, brick_6), function(x) as.character(crs(x)))

if(length(unique(all_extents)) == 1 && length(unique(all_resolutions)) == 1 && length(unique(all_crs)) == 1){
  # 属性一致，执行合并
  emi_brick_2018 <- brick(brick_1, brick_2, brick_3, brick_4, brick_5, brick_6)
  # 写入合并后的完整brick
  writeRaster(emi_brick_2018, file = file.path(stilt_wd, "calsacto/data", "PM25_2018Raster_Cal_Full.grd"), format="raster", overwrite=TRUE)
}else{
  stop("部分栅格的空间属性不一致，请检查对应chunk的生成逻辑")
}

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