检测栅格是否在空间多边形对象内部、外部或相交。

Question

检测栅格是否在空间多边形对象内部、外部或相交。

4

我有很多光栅数据，我想检查它们是否完全包含在一个空间多边形内部、完全不在空间多边形内部或与空间多边形相交（这可能意味着多边形完全在光栅数据中，或者多边形和光栅数据重叠）。我进行这个检查是为了节省时间，避免不必要的蒙版操作。

以下是示例：

                                      # 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))

绘制以下内容：

我将分别读取每个栅格，并检查它是否在空间多边形内、外或相交。

您认为在R中最有效的方法是什么？我有数千个4MB的栅格图像，计划并行进行掩膜处理，并希望通过此检查加快速度。

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

但是，我认为它没有给出我需要的详细信息。例如，所有栅格都在空间多边形的范围内，但并不是所有栅格都在空间多边形内部。

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

谢谢！

- Tedward

2个回答

3

以下是我解决问题的方法：

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")
  }
})

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

- Tedward

1

通过删除“print”调用和“proj4string”赋值，您将节省一点时间。 - jbaums

谢谢，打印调用只是占位符。我将调用一些其他函数。proj4string赋值做得好。 - Tedward

我的意思是你可以只写 'fully within'，而不是 print('fully within')，但如果它们只是占位符，那就不用担心了。 - jbaums

知道了，很不错。非常感谢！ - Tedward

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- jbaums · Accepted Answer

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

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

您可以将这些字符串与感兴趣的关系进行比较。例如，我们可以将within，disjoint和overlaps定义如下（但请注意，对于给定的关系，其他一些交集是可选的 - “within”由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"

它似乎比gContainsProperly/gIntersects方法稍微快一点，但@Tedward的帖子更易理解，并且与GEOS定义更一致（尽管创建特定关系定义的能力可能是可取的）。

DE-9IM字符串的元素按顺序表示：

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