如何提取栅格像素覆盖达指定比例的多边形？R代码问题排查

我已多次尝试解决该问题，但始终未找到理想方案。需编写R代码生成满足以下规则的矢量图层：每个多边形内值为1的栅格像素占其表面积的30%及以上。现有两个图层：

• 矢量图层BDPPAD2014_crop：代表某区域农业用地
• 栅格图层degraded（tif格式）：像素值标识农业用地退化区域

目标是提取退化区域占比≥30%的农田，但生成的矢量图层在QGIS中查看时，部分选中的多边形完全没有退化像素（无值为1的像素），不符合需求，恳请排查代码问题。

加载所需图层 #

BDPPAD2014_crop<-st_read(dsn=paste0(path_name, "Data/Vector/BDPPAD2014_crop.shp"))
degraded <- raster(paste0(path_name, "Data/Raster/negativeSignificant_trend.tif"))
plot(degraded)

检查CRS兼容性 #

# 检查CRS兼容性
if (!identical(crs(BDPPAD2014_crop), crs(degraded))) {
  print("CRS不匹配。正在将矢量图层投影到栅格图层的CRS。")
  BDPPAD2014_crop <- st_transform(BDPPAD2014_crop, crs(degraded))
}
# 检查空间范围兼容性
extent_vector <- st_bbox(BDPPAD2014_crop)
extent_raster <- extent(degraded)

if (!identical(extent_vector, extent_raster)) {
  print("空间范围不匹配。正在调整栅格范围以匹配矢量图层。")
  
  # 调整栅格范围匹配矢量
  extent_raster <- extent_vector
  extent(degraded) <- extent_raster
}
# 打印范围
cat("矢量图层范围:
", st_bbox(BDPPAD2014_crop), "
")
# 提取栅格范围
raster_extent <- extent(degraded)

# 打印栅格范围
print(raster_extent)

重分类栅格与筛选多边形 #

# 将所有退化像素设为1
m = c(-Inf, 0, 1)
rcl = matrix(m, ncol = 3, byrow = TRUE)
degraded.rcl <- reclassify(degraded, rcl)
# 使用exact_extract提取像素值
pixel_values <- exact_extract(degraded.rcl, BDPPAD2014_crop)

# 计算每个多边形内1值像素的占比
percentage_1s <- sapply(pixel_values, function(pixel_values) {
  if (length(pixel_values) == 0) return(NA)
  mean(pixel_values == 1, na.rm = TRUE) * 100
})

# 设置阈值百分比
threshold_percentage <- 30

# 根据阈值筛选多边形（同时使用存储面积的"SUP"列）
selected_polygons <- BDPPAD2014_crop[percentage_1s >= threshold_percentage & BDPPAD2014_crop$SUP >= threshold_percentage, ]

# 检查并打印不符合阈值的多边形
non_compliant_polygons <- BDPPAD2014_crop[percentage_1s < threshold_percentage | BDPPAD2014_crop$SUP < threshold_percentage, ]
print(non_compliant_polygons)

# 绘制选中多边形与原始图层对比
ggplot() +
  geom_sf(data = BDPPAD2014_crop, fill = "lightgray") +
  geom_sf(data = selected_polygons, fill = "blue")

st_write(selected_polygons,dsn=paste0(path_name, "Data/Vector/20230919_selected_degradedfield.shp"))

1. 栅格重分类逻辑反转：原代码m = c(-Inf, 0, 1)将所有≤0的像素设为1，与“退化区域值为1”的需求完全相反，导致占比计算完全错误。
2. 筛选条件逻辑混淆：BDPPAD2014_crop$SUP >= threshold_percentage将多边形面积（数值）与30%的占比阈值直接对比，逻辑无效，会错误过滤或保留多边形。
3. NA值处理漏洞：无像素覆盖的多边形返回NA，R子集筛选时会保留NA对应的多边形，导致无退化像素的多边形被错误选中。

# 加载所需包
library(sf)
library(raster)
library(exactextractr)
library(ggplot2)

# 加载图层
BDPPAD2014_crop <- st_read(dsn = paste0(path_name, "Data/Vector/BDPPAD2014_crop.shp"))
degraded <- raster(paste0(path_name, "Data/Raster/negativeSignificant_trend.tif"))
plot(degraded)

# 检查并统一CRS
if (!st_crs(BDPPAD2014_crop) == st_crs(degraded)) {
  message("CRS不匹配，正在投影矢量图层至栅格CRS")
  BDPPAD2014_crop <- st_transform(BDPPAD2014_crop, crs = st_crs(degraded))
}

# 裁剪栅格以匹配矢量范围（替代直接修改extent的错误操作）
degraded_cropped <- crop(degraded, extent(BDPPAD2014_crop))

# 重分类：仅将退化区域（原始值为1）保留为1，其余设为0
# 若原始退化栅格的退化值为其他数值，替换下面的1即可
m <- c(-Inf, 1, 0,
       1, 1, 1,
       1, Inf, 0)
rcl <- matrix(m, ncol = 3, byrow = TRUE)
degraded.rcl <- reclassify(degraded_cropped, rcl)

# 使用exact_extract计算加权占比（考虑像素被多边形覆盖的比例）
pixel_stats <- exact_extract(degraded.rcl, BDPPAD2014_crop, 
                             fun = function(values, coverage_fraction) {
                               total_weight <- sum(coverage_fraction, na.rm = TRUE)
                               if (total_weight == 0) return(0)
                               degraded_weight <- sum(coverage_fraction[values == 1], na.rm = TRUE)
                               (degraded_weight / total_weight) * 100
                             })

# 设置占比阈值
threshold_percentage <- 30

# 筛选符合条件的多边形：占比≥30%且存在退化像素
selected_polygons <- BDPPAD2014_crop[pixel_stats >= threshold_percentage & pixel_stats > 0, ]

# 对比绘图
ggplot() +
  geom_sf(data = BDPPAD2014_crop, fill = "lightgray", alpha = 0.5) +
  geom_sf(data = selected_polygons, fill = "blue", alpha = 0.7) +
  labs(title = "退化占比≥30%的农田")

# 保存结果
st_write(selected_polygons, dsn = paste0(path_name, "Data/Vector/20230919_selected_degradedfield_corrected.shp"), delete_layer = TRUE)

• 修正重分类规则：根据原始栅格的退化值定义，正确将退化区域设为1，其余为0。
• 裁剪栅格而非修改extent：使用crop()函数裁剪栅格到矢量范围，避免数据错位。
• 加权计算占比：利用coverage_fraction参数计算像素覆盖比例，得到更精确的面积占比。
• 修复筛选逻辑：移除面积与百分比的错误对比，同时排除占比为0的多边形。
• 处理无覆盖情况：无栅格覆盖的多边形返回占比0，避免NA导致的错误筛选。

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