如何检查栅格和多边形是否重叠？

Question

我有一个栅格，我想检查具有值的单元格是否与我感兴趣的区域重叠，这是一个空间多边形文件。

有没有快速的方法来做到这一点？

我的光栅看起来像这样：

我尝试了什么：

r <- raster("C:/user/test.jp2")
studyarea <- readOGR("C:/user/boundary.shp")
dummi <- as(extent(r), "SpatialPolygons")
if (gContainsProperly(studyarea, dummi)) {
  print("Raster extent is fully within the Area of interest")
} 
else if (gIntersects(studyarea, dummi)) {
  print("Raster extent is not fully within the Area of interest")
} else {
  print("Raster extent is fully outside the Area of interest")
}

as(extent(r), "SpatialPolygons") 在值和 NA 内创建一个多边形。但我只想检查值部分和研究区域是否重叠

栅格和研究区域的数据示例：

matrix <- matrix(c(NA,NA,NA,1,NA,NA,2,1,NA,1,3,4,2,3,4,1), nrow=4)
r <- raster(matrix)
extent(r) <- c(399960, 509760, 5290200, 5400000)
crs(r) <- "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs" 
writeRaster(r, paste(getwd(),"filename.jp2"))
                 
matrix <- matrix(c(1,1,NA,NA,1,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA), nrow=4)
rpoly<- raster(matrix)
extent(rpoly) <- c(399960, 509760, 5290200, 5400000)
crs(rpoly) <- "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs" 
studyarea <- rasterToPolygons(rpoly)

谢谢

Answer 1

所以这里是一个使用数据示例的方法。

想法是将具有值的栅格单元转换为单个多边形，并将该多边形与您的研究区域多边形进行比较。

sf 提供多种谓词来比较几何。

matrix <- matrix(c(NA,NA,NA,1,NA,NA,2,1,NA,1,3,4,2,3,4,1), nrow=4)
r <- raster::raster(matrix)
raster::extent(r) <- c(399960, 509760, 5290200, 5400000)
raster::crs(r) <- "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs" 
r <- sf::st_as_sf(raster::rasterToPolygons(r)) # This will create a sfc polygon dataframe from your cells with values
r <- sf::st_union(r) # This will combine all cells with values to a single polygon

matrix <- matrix(c(1,1,NA,NA,1,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA), nrow=4)
rpoly<- raster::raster(matrix)
raster::extent(rpoly) <- c(399960, 509760, 5290200, 5400000)
raster::crs(rpoly) <- "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs" 
studyarea <- sf::st_as_sf(raster::rasterToPolygons(rpoly)) # This will create a sfc polygon dataframe from your cells with values 
studyarea <- sf::st_union(studyarea) # This will combine all cells with values to a single polygon

sf::st_contains(r, studyarea, sparse = FALSE) #  This checks if studyarea is within r, see ?sf::st_contains for other predicates you might want (touches, crosses, etc.)

Answer 2

我找到了一个更快的解决方案。可以通过多边形裁剪栅格，然后用多边形掩盖栅格。如果多边形区域中没有栅格值。栅格将具有最小值 = 0 和最大值 = 0。

创建数据

matrix <- matrix(c(NA,NA,NA,1,NA,NA,2,1,NA,1,3,4,2,3,4,1), nrow=4)
r <- raster(matrix)
extent(r) <- c(399960, 509760, 5290200, 5400000)
crs(r) <- "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs" 
writeRaster(r, paste(getwd(),"filename.jp2"))
                 
matrix <- matrix(c(1,1,NA,NA,1,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA,NA), nrow=4)
rpoly<- raster(matrix)
extent(rpoly) <- c(399960, 509760, 5290200, 5400000)
crs(rpoly) <- "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs" 
studyarea <- rasterToPolygons(rpoly)

解决方案

# Reproject. Make sure both have the same projection
studyarea <-  spTransform(studyarea, projection(r))
  
# Crop and clip raster with study area polygon
r_crop <- crop(r, studyarea)
r_mask <- mask(r_crop, studyarea, updatevalue = NA, updateNA= TRUE)

if (r_mask@data@min == 0 && r_mask@data@max == 0){
  print("raster and studyarea do not overlap))
  }
else {
  print("Raster and studyarea do overlap")
}