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

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

我有一个栅格，我想检查具有值的像元是否与我的关注区域（一个空间多边形文件）重叠。

我的栅格看起来像这样：

我尝试了什么：

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,1,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,NA),nrow=4)
rpoly<- raster(matrix)
extent(rpoly) <- c(399960,5400000)
crs(rpoly) <- "+proj=utm +zone=32 +datum=wgs84 +units=m +no_defs" 
studyarea <- rasterTopolygons(rpoly)

谢谢

解决方法

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

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

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

matrix <- matrix(c(NA,NA,1,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,NA),nrow=4)
rpoly<- raster::raster(matrix)
raster::extent(rpoly) <- c(399960,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.)

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

创建数据

matrix <- matrix(c(NA,nrow=4)
r <- raster(matrix)
extent(r) <- c(399960,5400000)
crs(r) <- "+proj=utm +zone=32 +datum=WGS84 +units=m +no_defs" 
writeRaster(r,paste(getwd(),"filename.jp2"))
                 
matrix <- matrix(c(1,nrow=4)
rpoly<- raster(matrix)
extent(rpoly) <- c(399960,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,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")
}