R语言绘制栅格与Shapefile箱线图遇depth方法NULL类错误求助

正常运行的示例代码 #

以下示例代码可正常运行：

library(raster)
library(sp)
library(rgdal)
library(ggplot2)

# 创建示例栅格数据
raster_file <- raster(ncol=36, nrow=18)
raster_file[] <- 1:ncell(raster_file)
plot(raster_file)

# 创建示例多边形
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0)) 
shape_file <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1), 
                          Polygons(list(Polygon(cds2)), 2)))
plot(shape_file)

# 提取多边形内的栅格值
extracted_values <- extract(raster_file, shape_file)

# 为每个多边形构建数据框并合并
data_list <- lapply(1:length(extracted_values), function(i) {
  data.frame(value = extracted_values[[i]], polygon = i)
 })
data <- do.call(rbind, data_list)

# 绘制箱线图
bp <- ggplot(data, aes(x = factor(polygon), y = value)) +
  geom_boxplot() +
  labs(x = "Polygon", y = "Raster Values") +
  theme_minimal()
bp

自有数据集的操作及报错 #

使用自有数据集时，读取多边形后绘制箱线图出现报错，操作过程及报错如下：

# 加载shapefile并查看字段
shape_file <- shapefile("C:/Users/.....BiogeoRegion.shp")
Warning message:
[vect] Z coordinates ignored 

> names(shape_file)
 [1] "RegionNumm" "RegionName" "Unterregio" "Unterreg_1" "ObjNummer"
[6] "Version"    "Shape_Leng" "Shape_Area" "DERegionNa" "FRRegionNa"
[11] "ITRegionNa" "DEBioBedeu" "FRBioBedeu" "ITBioBedeu"

> str(shape_file)
Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 12 obs. of  14 variables:
  .. ..$ RegionNumm: int [1:12] 1 2 2 2 2 3 3 4 5 5 ...
  .. ..$ RegionName: chr [1:12] "R1" "R2" "R2" "R2" ...
  .. ..$ Unterregio: int [1:12] 11 21 22 23 24 31 32 41 51 52 ...
  .. ..$ Unterreg_1: chr [1:12] "U11" "U21" "U22" "U23" ...
  .. ..$ ObjNummer : chr [1:12] "1" "2" "3" "4" ...
  .. ..$ Version   : chr [1:12] "2020/05/08" "2020/05/08" "2020/05/08" "2020/05/08" ...

# 从列表中选择栅格文件
raster_file <- allrasters_pres[[1]]
> names(raster_file)
[1] "Andrena.barbilabris_glo_ensemble"

# 提取多边形内的栅格值
extracted_values <- extract(raster_file, shape_file)

# 构建箱线图数据
data_list <- lapply(1:length(extracted_values), function(i) {
    data.frame(value = extracted_values[[i]], polygon = i)
})
data <- do.call(rbind, data_list)

> names(data)
[1] "value"   "polygon"

# 绘制箱线图
bp <- ggplot(data, aes(x = factor(polygon), y = value)) +
    geom_boxplot() +
    labs(x = "Polygon", y = "Raster Values") +
    theme_minimal()
> bp
Error in UseMethod("depth") : 
  no applicable method for 'depth' applied to an object of class "NULL"
In addition: Warning message:
Removed 452451 rows containing non-finite outside the scale range (`stat_boxplot()`). 

解决技术方案 #

错误核心原因 #

1. 提取的栅格值中存在大量非有限值（NA/Inf/-Inf），stat_boxplot()过滤后部分多边形分组没有有效数据，无法计算箱线图统计量，触发depth方法的NULL错误。
2. 可能存在栅格与多边形坐标系不匹配，或部分多边形与栅格无空间重叠，导致提取结果全为NA。

分步解决 #

1. 排查提取结果的有效性

先确认每个多边形的提取值情况，找出无有效数据的分组：

# 查看每个多边形的提取值总数、NA数量及有效占比
lapply(extracted_values, function(x) {
  total <- length(x)
  na_count <- sum(is.na(x))
  list(total_count = total, na_count = na_count, valid_ratio = (total - na_count)/total)
})

2. 清理无效数据

构建数据时直接过滤NA值，并跳过无有效数据的多边形：

# 生成数据列表时仅保留有有效值的分组
data_list <- lapply(1:length(extracted_values), function(i) {
  vals <- extracted_values[[i]]
  # 筛选非NA且为有限值的记录
  valid_vals <- vals[!is.na(vals) & is.finite(vals)]
  # 只有有效值数量>0时才生成数据框
  if (length(valid_vals) > 0) {
    data.frame(value = valid_vals, polygon = i)
  }
})
# 移除列表中的空元素
data_list <- Filter(Negate(is.null), data_list)
# 合并数据
data <- do.call(rbind, data_list)

3. 检查并统一坐标系

确认栅格和多边形的坐标系是否一致，不一致则转换：

# 查看两者的坐标系
print(crs(raster_file))
print(crs(shape_file))

# 如果坐标系不匹配，将多边形转换为栅格的坐标系
if (crs(raster_file) != crs(shape_file)) {
  shape_file <- spTransform(shape_file, crs(raster_file))
  # 重新提取栅格值
  extracted_values <- extract(raster_file, shape_file)
}

4. 绘制箱线图

使用清理后的数据绘图，加上na.rm = TRUE确保忽略剩余NA：

bp <- ggplot(data, aes(x = factor(polygon), y = value)) +
  geom_boxplot(na.rm = TRUE) +
  labs(x = "Polygon", y = "Raster Values") +
  theme_minimal()
print(bp)

内容的提问来源于stack exchange，提问作者Sibylle Stöckli