如何执行两个SpatialPolygonsDataFrame对象的矢量叠加？(How to perf

我有两个GIS图层-给他们打电话Soils和Parcels -存储为SpatialPolygonsDataFrame S（ SPDF S），我想为“叠加”起来， 在这里所描述的意义 。

覆盖操作的结果应该是一个新的SPDF：

1. 所述SpatialPolygons组分包含由两个层的相交处形成的多边形。 （认为​​所有由高架投影机上叠置两个mylars形成的原子多边形）。

2. 所述data.frame部件记录的属性Soils和Parcels的多边形内，其中每个原子多边形下降。

我的问题（S）：有没有做这个工作的现有研发功能？ （我甚至高兴地得知，刚刚得到一个函数会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)))
})

Answer 1:

自2014年1月， 光栅包包括union()函数，使得这个易于peasy：

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)

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)

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

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， Parcels <- spChFIDs(Parcels, paste0("pp", row.names(Parcels)))

Answer 3:

这里是我的裂纹，这只是给出了与地块（Parcels->土壤）的数据碎片的列表。 你仍然需要得到从土壤对象的属性，然后做类似的工作与“差异”，然后反之亦然（Soils->包裹），所以你可以有任何种类的重叠关系。

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

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 。

文章来源: How to perform a vector overlay of two SpatialPolygonsDataFrame objects?