R: Различные результаты от intersect () и gIntersects () при проверке, пересекает ли полигон растр

Question

Основано на этом посте (https://gis.stackexchange.com/questions/255025/r-raster-masking-a-raster-by-polygon-also-remove-cells-partially-covered) Я настроил следующие растры и один SpatialPolygonsDataFrame, а затем использовал raster::intersect и rgeos::gIntersects, чтобы проверить, включают ли растры полигон:

library(raster)
library(sp)
library(rgeos)

# create rasters and store them in a list
r1 <- raster(xmn=1, xmx=5, ymn=1, ymx=5, nrows=4, ncols=4)
r1[] <- 1:length(r1)
r2 <- raster(xmn=5, xmx=9, ymn=1, ymx=5, nrows=4, ncols=4)
r2[] <- 10:(length(r2)+9)
r3 <- raster(xmn=1, xmx=5, ymn=5, ymx=9, nrows=4, ncols=4)
r3[] <- seq(0,1.5,0.1)
r_list <- list(r1,r2,r3)
# create SpatialPolygonsDataFrame
Sr1 = Polygon(cbind(c(2,4,4,1,2),c(2,3,5,4,2)))
Sr2 = Polygon(cbind(c(1,2,2),c(1,2,1)))
SpP = SpatialPolygons(list(Polygons(list(Sr1), "s1"), Polygons(list(Sr2), "s2")), 1:2)
dat = data.frame(ID = c("s1", "s2"), value = c("a", "b"))
row.names(dat) <- c("s1", "s2")
p <- SpatialPolygonsDataFrame(SpP, data = dat, 
                              match.ID = TRUE)

# check if rasters include the polygon
for (i in 1:length(r_list)) {

  inter1 <- raster::intersect(extent(r_list[[i]]), extent(p))
  print(paste0("p intersects r", i, ": raster::intersect ", isTRUE(inter1)) )
  inter2 <- gIntersects(as(extent(r_list[[i]]), 'SpatialPolygons'), p)
  print(paste0("p intersects r", i, ": rgeos::gIntersects ", inter2) )

}

Почему raster::intersect не распознает, что многоугольник покрывает часть r1 (и r3), но gIntersects делает это правильно?

Answer 1

Я получаю эти результаты

raster::intersect(extent(r1), extent(p))
#class      : Extent 
#xmin       : 1 
#xmax       : 4 
#ymin       : 1 
#ymax       : 5 
gIntersects(as(extent(r1), 'SpatialPolygons'), p)
#[1] TRUE

raster::intersect(extent(r2), extent(p))
#NULL
gIntersects(as(extent(r2), 'SpatialPolygons'), p)
#[1] FALSE

raster::intersect(extent(r3), extent(p))
#NULL
gIntersects(as(extent(r3), 'SpatialPolygons'), p)
#[1] TRUE

Так что только с r3 есть разница.Это потому, что вы можете спорить, есть ли пересечение или нет, так как многоугольники только соприкасаются друг с другом.Вы можете видеть, что gIntersection возвращает один SpatialPoint

gIntersection(as(extent(r3), 'SpatialPolygons'), p)
#class       : SpatialPoints 
#features    : 1 
#extent      : 4, 4, 5, 5  (xmin, xmax, ymin, ymax)
#crs         : NA 

Для raster пересечение между полигонами также должно быть многоугольником, поэтому пересечение здесь отсутствует.