如何在两个SpatialPolygonsDataFrame对象之间执行矢量叠加？

Question

我有两个GIS图层 -- 分别称为Soils和Parcels -- 存储为SpatialPolygonsDataFrame（SPDF），我想要"叠加"它们，如此处所述。

叠加操作的结果应该是一个新的SPDF，其中:

1. SpatialPolygons组件包含由两个图层交集形成的多边形。(想象一下在投影仪上重叠两个mylar时形成的所有原子多边形)。

2. data.frame组件记录每个原子多边形所在的Soils和Parcels多边形的属性。

我的问题：是否存在一个现有的R函数可以做到这一点？（我甚至会很高兴地了解一个只能正确获取SpatialPolygons组件的函数，计算从两个图层相交形成的原子多边形。）我觉得rgeos应该至少有一个可以实现（1）的函数，但似乎并没有... 这里是一张图，可能有助于更清楚地理解我想要的内容，接下来是创建在图中显示的Soils和Parcels图层的代码。

library(rgeos)

## Just a utility function to construct the example layers.
flattenSpatialPolygons <- function(SP) {
    nm <- deparse(substitute(SP))
    AA <- unlist(lapply(SP@polygons, function(X) X@Polygons))
    SpatialPolygons(lapply(seq_along(AA),
                           function(X) Polygons(AA[X], ID=paste0(nm, X))))
}

## Example Soils layer
Soils <-
local({
    A <- readWKT("MULTIPOLYGON(((3 .5,7 1,7 2,3 1.5,3 0.5), (3 1.5,7 2,7 3,3 2.5,3 1.5)))")
    AA <- flattenSpatialPolygons(A)
    SpatialPolygonsDataFrame(AA,
           data.frame(soilType = paste0("Soil_", LETTERS[seq_along(AA)]),
                      row.names = getSpPPolygonsIDSlots(AA)))
})

## Example Parcels layer
Parcels <-
local({
    B <- readWKT("MULTIPOLYGON(((0 0,2 0,2 3,0 3,0 0),(2 0,4 0,4 3,2 3,2 0)),((4 0,6 0,6 3,4 3,4 0)))")
    BB <- flattenSpatialPolygons(B)
    SpatialPolygonsDataFrame(BB,
           data.frame(soilType = paste0("Parcel_", seq_along(BB)),
                      row.names = getSpPPolygonsIDSlots(BB)))
})

- Josh O'Brien

---

%over% 还是 overlay？请参考 gis.SE 和 @Spacedman 的备忘单：http://www.maths.lancs.ac.uk/~rowlings/Teaching/UseR2012/cheatsheet.html - Ari B. Friedman

@AriB.Friedman -- 谢谢，但不是这个-- over() 不是我想要的。就像我说的，我需要一个返回单个多边形的函数，而不仅仅是它们是否重叠的指示器。（此外，即使我可以构造原子多边形，over 等也没有帮助，因为它们将共享边界的多边形计算为彼此“覆盖”。rgeos::gRelate() 在这方面更好。） - Josh O'Brien

gContains和gOverlaps，然后是gIntersection？ - Ari B. Friedman

@AriB.Friedman -- 能给我展示一下你的意思吗？在放弃并发布之前，gIntersection和gSymdifference的组合是我得到的最接近的结果，但它没有获取原子多边形（即底部图中不包含任何线条的10个多边形）。这似乎是一个非常基本的GIS操作，我将很高兴向展示如何在R中完成它（或通过调用其他易于链接的GIS）的人提供赏金... - Josh O'Brien

这周非常忙，但希望下面的答案足以让你开始。如果不行，请告诉我，我会尽力做更多。 - Ari B. Friedman

Answer 1

自2014年1月，raster库包括了一个union()函数，使得这个过程变得轻而易举：

library(raster)
Intersects <- Soils + Parcels  ## Shorthand for union(Soils, Parcels)

## Check that it work
data.frame(Intersects)
## soilType.1 soilType.2
## 1     Soil_A       <NA>
## 2     Soil_B       <NA>
## 3       <NA>   Parcel_1
## 4       <NA>   Parcel_2
## 5       <NA>   Parcel_3
## 6     Soil_A   Parcel_2
## 7     Soil_A   Parcel_3
## 8     Soil_B   Parcel_2
## 9     Soil_B   Parcel_3
plot(Intersects, col = blues9)

- Josh O'Brien

Answer 2

自2014年1月起，此处发布的解决方案已被完全取代，现可使用 raster::union() 函数，该函数在上面正式接受的答案中进行了演示。

---

此代码获取通过叠加两个不同的 SpatialPolygons 图层形成的“原子”多边形，解决了我问题的第一部分。

library(rgeos)

gFragment <- function(X, Y) {
    aa <- gIntersection(X, Y, byid = TRUE)
    bb <- gDifference(X, gUnionCascaded(Y), byid = T)
    cc <- gDifference(Y, gUnionCascaded(X), byid = T)
    ## Note: testing for NULL is needed in case any of aa, bb, or cc is empty,
    ## as when X & Y don't intersect, or when one is fully contained in the other
    SpatialPolygons(c(if(is.null(aa)) NULL else aa@polygons,
                      if(is.null(bb)) NULL else bb@polygons,
                      if(is.null(cc)) NULL else cc@polygons)) 
}

## Try it out
Fragments <- gFragment(Parcels, Soils)
plot(Fragments, col=blues9)

这会提取每个由gFragment()输出的“原子”多边形所覆盖的两个输入图层中多边形的ID（如果有），解决了我问题的第二部分。

getAttributes <- function(Fragments, Layer1, Layer2, eps = 0) {
    ## Function to extract attributes from polygon layers
    ## overlain by fragments
    OVER <- function(AA, BB) {
        X <- gRelate(AA, BB, byid = TRUE, pattern="2********")
        ii <- sapply(seq_len(ncol(X)),
                    function(i) {
                        A <- which(X[,i])
                        if(!length(A)) NA else A
                    })
        rownames(X)[ii]
    }
    ## First need to (very slightly) trim Fragments b/c otherwise they
    ## tend to (very slightly) overlap adjacent ring(s)
    Frags <- gBuffer(Fragments, width = -eps, byid = TRUE)
    ## Construct data.frame of attributes
    df <- data.frame(OVER(Frags, Layer1), 
                     OVER(Frags, Layer2),
                     row.names = names(Fragments))
    names(df) <- c(deparse(substitute(Layer1)), deparse(substitute(Layer2)))
    ## Add data.frame to SpatialPolygons object
    SpatialPolygonsDataFrame(Fragments, data=df)
}

FragmentsDF <- getAttributes(Fragments = Fragments,
                             Layer1 = Parcels,
                             Layer2 = Soils)

## A few ways to examine the results
data.frame(FragmentsDF, row.names=NULL)
#   Parcels Soils
# 1      B2    A1
# 2      B2    A2
# 3      B3    A1
# 4      B3    A2
# 5      B1  <NA>
# 6      B2  <NA>
# 7      B3  <NA>
# 8    <NA>    A1
# 9    <NA>    A2
spplot(FragmentsDF, zcol="Soils", col.regions=blues9[3:4])
spplot(FragmentsDF, zcol="Parcels", col.regions=grey.colors(3))

编辑：

请注意，如果您的输入多边形中有名称为"1"的id，则此代码可能会失败。在这种情况下，一个解决方法是重新命名id，例如通过执行以下操作：Parcels <- spChFIDs(Parcels, paste0("pp", row.names(Parcels)))。

- Josh O'Brien

---

@mdsumner 谢谢。rgeos相当不错，是吧？ - Josh O'Brien

Answer 3

这里是我的建议，以下内容仅列出了捆绑数据中的一份清单（Parcels->Soils）。您还需要从Soils对象中获取属性，并执行类似的“差异”工作，反之亦然（Soils-> Parcels），以便可以拥有各种重叠关系。

intersects <- list()

## find all intersections (NULLs do nothing to the result)
for (i in 1:nrow(Soils)) {
  for (j in 1:nrow(Parcels)) {
    intersects[[sprintf("%sx%s", i, j)]] <- gIntersection(Soils[i,],
                                                          Parcels[j,]) 
  }
}

result <- list()
## let's try Parcels, transfer data attributes to new pieces
for (i in 1:nrow(Parcels)) {
  for (j in seq_along(intersects))
   if(gContains(Parcels[i,], intersects[[j]])) {
     result <- c(result, SpatialPolygonsDataFrame(intersects[[j]],     as.data.frame(Parcels[i,]), match.ID = FALSE))

   }
}


## plot
plot(Parcels, xlim = range(c(bbox(Parcels)[1,], bbox(Soils[1,]))),
     ylim = range(c(bbox(Parcels)[2,], bbox(Soils[2,]))))
plot(Soils, add = TRUE)

cols <- colorRampPalette(c("lightblue", "darkblue"))(length(result))
for (i in 1:length(result)) plot(result[[i]], col = cols[i], add = TRUE)
for (i in 1:length(result)) text(coordinates(result[[i]]), label =     as.data.frame(result[[i]])[,"soilType"])

- mdsumner

---

这里有很多有趣的想法。谢谢！我将来会同时使用您的sprintf()和多边形标签代码。 - Josh O'Brien

Answer 4

这里是基本的想法（将其作为嵌套循环在包裹和土壤上执行；我不确定它是否可以像g*函数那样进行向量化编写）：

i <- 2
j <- 2
pieces <- list()
pieces[["int"]] <- gIntersection(Parcels[i,],Soils[j,])
pieces[["diff1"]] <- gDifference(Parcels[i,],Soils[j,])
pieces[["diff2"]] <- gDifference(Soils[j,],Parcels[i,])
plot(Parcels)
plot(Soils,add=TRUE)
plot(pieces[["int"]],col="red",add=TRUE)
plot(pieces[["diff1"]],col="blue",add=TRUE)
plot(pieces[["diff2"]],col="green",add=TRUE)

这应该足以让你开始了。其余的只是在跟踪这些部分并将它们合并成一个大的SPDF时进行循环。

更矢量化的另一种方法是：

pieces <- list()
pieces[["int"]] <- gIntersection(Parcels,Soils,byid=TRUE)
pieces[["diff1"]] <- gDifference(Parcels,Soils,byid=TRUE)
pieces[["diff2"]] <- gDifference(Soils,Parcels,byid=TRUE)

出于某种原因，这会给你比实际存在的更多的部分。然后你需要返回并rbind它们并剔除那些gEquals的额外部分。

- Ari B. Friedman

---

嗨，Ari - 非常感谢你参与解决这个问题，但是这并不是我需要的。蓝色和绿色的多边形都不是我想要的原子多边形（即我底部图中的10个阴影区域）。 （顺便说一下，我得到了那个图，通过类似于您的建议的方法，使用gIntersection（Soils，Parcels，byid = TRUE）和gSymdifference（Soils，Parcels，byid = TRUE）的组合：byid参数避免了您提到的大量/全部循环。但根本问题是，该代码也给我带来了除我想要的之外的其他多边形。） - Josh O'Brien

我认为如果你以递归的方式解决，它仍然可以实现你想要的效果。例如，绿色的多边形与Soils[1，]相交，这会将其进一步分割成所需的部分。虽然不太优雅，但我认为它仍然有效，对吧？ - Ari B. Friedman