如何仅提取按设定百分比的光栅像素覆盖的多边形？

Question

我在很多地方寻求有关这个问题的帮助，但我似乎找不到我正在寻找的答案。我想创建一个 R 代码来获取包含遵守此规则的所有多边形的矢量图层：每个多边形需要有 30% 或更多的表面覆盖像素，值为 1。为了给您更多的背景信息，我有两个图层（矢量：BDPPAD2014_crop=代表一个地区的农田）（栅格：degraded=代表一个包含像素的tif文件，显示这些田地是否正在退化）。我的目标是能够使用 R 提取 30% 或更多多边形区域显示退化的所有字段。我能够使用下面的代码获得矢量图层，但是当我在 QGIS () 中查看它时，我可以清楚地看到所选的一些多边形上没有像素，这不是我想要的，因为像素这里的意思是退化。有人可以帮助我了解我的代码中的问题出在哪里吗？感谢您的宝贵时间和帮助:)

我尝试了很多不同的方法，但这就是几乎有效的方法：

加载所需图层

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兼容性

# Check the CRS compatibility
if (!identical(crs(BDPPAD2014_crop), crs(degraded))) {
  print("CRS mismatch. Reprojecting one of the layers to match the other's CRS.")
  # Reproject BDPPAD2014_crop to match degraded's CRS
  BDPPAD2014_crop <- st_transform(BDPPAD2014_crop, crs(degraded))
}
# Check spatial extent compatibility
extent_vector <- st_bbox(BDPPAD2014_crop)
extent_raster <- extent(degraded)

if (!identical(extent_vector, extent_raster)) {
  print("Spatial extent mismatch. Adjusting the extent of the raster to match the vector.")
  
  # Adjust the extent of the raster to match the vector
  extent_raster <- extent_vector
  extent(degraded) <- extent_raster
}
# Print the extents
cat("Extent of Vector Layer:\n", st_bbox(BDPPAD2014_crop), "\n")
# Extract the extent
raster_extent <- extent(degraded)

# Print the extent
print(raster_extent)

重新分类栅格

# All degraded pixels will have a value of 1
m = c(-Inf, 0, 1)
rcl = matrix(m, ncol = 3, byrow = TRUE)
degraded.rcl <- reclassify(degraded, rcl)
# Extract pixel values using exact_extract
pixel_values <- exact_extract(degraded.rcl, BDPPAD2014_crop)

# Calculate the percentage of 1s within each polygon
percentage_1s <- sapply(pixel_values, function(pixel_values) {
  if (length(pixel_values) == 0) return(NA)
  mean(pixel_values == 1, na.rm = TRUE) * 100
})

# Set the threshold percentage
threshold_percentage <- 30

# Filter polygons based on the threshold percentage using column "SUP" which contain the area of each polygons
selected_polygons <- BDPPAD2014_crop[percentage_1s >= threshold_percentage & BDPPAD2014_crop$SUP >= threshold_percentage, ]

# Check and print polygons that don't meet the threshold
non_compliant_polygons <- BDPPAD2014_crop[percentage_1s < threshold_percentage | BDPPAD2014_crop$SUP < threshold_percentage, ]
print(non_compliant_polygons)

# Plot the selected polygons vs the BDPPAD2014_crop 
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"))

Answer 1

追踪裁剪数据，组成一个退化的栅格（（单元格值设置为1，然后通过采样重置），摆弄让

crs

和

extent

如上排列，并从每个中取出一小部分裁剪部分，包裹和降级，以下方法可能有效

library(terra)
# some objects

# smaller_agri = draw(x = 'extent') # get small num of parcels from all BD2014
# smaller_agri
#SpatExtent : -73.0458375577738, -72.9682494683782, 45.8925658029653, 45.9361554490149 (xmin, xmax, ymin, ymax)
# smaller_agri2 = crop(agri_parcels, smaller_agri)
# small_agr2_rast = rast(smaller_agri2, extent = ext(smaller_agri2), nrows=1000, ncol=1000, crs = crs(smaller_agri2))
# degraded = rast(nrows = 1000, ncols = 1000, extent = ext(smaller_agri), crs = crs(smaller_agri))
# values(degraded) = 1
# degraded2 = crop(degraded, small_agri2, mask = TRUE) # also replaces your `classify`

# degraded2_1[not.na(degraded2_1)] = sample(c(0,1), length(which(values(degraded2_1) == 1)), replace = TRUE, prob = c(0.71, 0.29))
# degraded2_1_crs = project(degraded2_1, small_agr2_rast) #set crs
# degrade_mean = terra::extract(degraded2_1_crs, smaller_agri2, FUN = mean)
dim(degrade_mean)
[1] 827209      2 # there are lots of layers...
names(degrade_mean)
[1] "ID"    "lyr.1" 

# this should get us there and can `<- crop(` of course
plot(crop(degraded2_1_crs, smaller_agri2[degrade_mean[,2] >= 0.29], mask = TRUE)) # which layers
plot(smaller_agri2[degrade_mean[, 2] >= 0.29], add = TRUE)

认为这就是您正在寻找的。