检测栅格是否在 SpatialPolygons 对象内、不在或与 SpatialPolygons 对象相交

Question

我有很多栅格，我想检查它们是否完全包含在空间多边形内、完全没有空间多边形或与空间多边形相交（这可能意味着多边形完全在栅格内，或多边形和栅格重叠）。我正在做这项检查，以便在可能的情况下避免时间密集型屏蔽。

这是一个例子：

                                      # create 3 example rasters
r <- raster()
r[] <- rnorm(n = ncell(r))
e1 <- extent(c(45,55,45,50))
r1 <- crop(r,e1)
e2 <- extent(c(20,25,25,30))
r2 <- crop(r,e2)
e3 <- extent(c(38,55,57,65))
r3 <- crop(r,e3)


                                        #create SpatialPolygons

x <- c(40,60)
y <- c(40,60)
m <- expand.grid(x,y)
m <- m[c(1,2,4,3),]
p1 <- Polygon(m)
p1 <- Polygons(list(p1),1)

x <- c(10,15)
y <- c(10,15)
m <- expand.grid(x,y)
m <- m[c(1,2,4,3),]
p2 <- Polygon(m)
p2 <- Polygons(list(p2),2)

x <- c(30,45)
y <- c(70,80)
m <- expand.grid(x,y)
m <- m[c(1,2,4,3),]
p3 <- Polygon(m)
p3 <- Polygons(list(p3),3)


poly <- SpatialPolygons(list(p1,p2,p3))

绘制这些：

我将分别读取每个栅格并检查它是否在 SpatialPolygons 之内、之外或相交。

您认为在 R 中执行此操作的最有效方法是什么？我有数以千计的 4mb 光栅，我计划并行屏蔽这些光栅，希望此检查能够加快速度。

注意，还有这个问题：https://gis.stackexchange.com/questions/34535/detect-whether-there-is-a-spatial-polygon-in-a-spatial-extent

但是，我认为它没有提供我正在寻找的详细信息。例如，所有栅格都在空间多边形的 范围内，但并非所有栅格都在空间多边形内。

rgeos 中的函数（gIntersects、gContains）可能会很方便。我不确定这些是否最有效，或者我应该如何将栅格（或其范围）转换为 sp 对象。

谢谢！

Answer 1

这是我解决问题的方法：

library(rgeos)
rlist <- list(r1,r2,r3)

lapply(rlist, function(raster) {
  ei <- as(extent(raster), "SpatialPolygons")
  if (gContainsProperly(poly, ei)) {
    print ("fully within")
  } else if (gIntersects(poly, ei)) {
    print ("intersects")
  } else {
    print ("fully without")
  }
})

如果您知道更有效的解决方案，请告诉我。

Answer 2

您也可以为此使用 gRelate。它 returns 一个 DE-9IM 代码，描述了两个几何图形的内部、边界和外部组件之间的关系。

library(rgeos)
x <- sapply(rlist, function(x) 
  gsub('[^F]', 'T', gRelate(as(extent(x), 'SpatialPolygons'), poly)))

然后您可以将字符串与感兴趣的关系进行比较。例如，我们可以如下定义 within、disjoint 和 overlaps（但请注意，对于给定关系，其他一些交集是可选的 - "within" 是 defined by GEOS T*F**F***、"disjoint" 为 FF*FF****，"overlaps" 为 T*T***T**）：

pat <- c(TFFTFFTTT='within', FFTFFTTTT='disjoint', TTTTTTTTT='overlaps')
pat[x]

##  TFFTFFTTT  FFTFFTTTT  TTTTTTTTT 
##   "within" "disjoint"  "overlaps"

它似乎比 , but @Tedward's 稍微快一点，更易于理解，并且更符合 GEOS 定义（尽管创建特定关系定义的能力可能是可取的）。

DE-9IM字符串的元素依次表示：

几何体 A 的内部是否与几何体 B 的内部相交？
A 的边界是否与 B 的内部相交？
A 的外部是否与 B 的内部相交？
几何体A的内部是否与几何体B的边界相交？
A的边界是否与B的边界相交？
A的外部是否与B的边界相交？
几何体 A 的内部是否与几何体 B 的外部相交？
A 的边界是否与 B 的外部相交？
A 的外部是否与 B 的外部相交？

检测栅格是否在 SpatialPolygons 对象内、不在或与 SpatialPolygons 对象相交

Detect if raster is within, without, or intersecting a SpatialPolygons object

r

r-raster

sp