r - 如何执行两个 SpatialPolygonsDataFrame 对象的矢量叠加？

Question

我有两个 GIS 层——调用它们Soils并Parcels存储为SpatialPolygonsDataFrames( SPDFs)，并且我想“覆盖”它们，在此处描述的意义上。

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

该SpatialPolygons组件包含由两层相交形成的多边形。（想想在投影仪上覆盖两个聚酯薄膜形成的所有原子多边形）。
该组件记录每个原子多边形所属的多边形和多边形data.frame的属性。SoilsParcels

我的问题：是否存在执行此操作的现有 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)))
})

score 7 · Accepted Answer

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

在此处输入图像描述

score 4 · Accepted Answer

自 2014 年 1 月以来，此处发布的解决方案已完全被该raster::union()功能取代，如上面官方接受的答案所示。

这得到了通过覆盖两个不同SpatialPolygons层形成的“原子”多边形，解决了我问题的第 1 部分。

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（如果有），解决了我的问题的第 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))

编辑：

请注意，如果您的任何输入多边形具有 id 的 named ，则此代码可能会失败"1"。在这种情况下，一种解决方法是重命名 id，也许通过执行类似Parcels <- spChFIDs(Parcels, paste0("pp", row.names(Parcels))).

score 2 · Accepted Answer

这是我的破解，这只是给出了带有包裹（包裹->土壤）数据的碎片列表。您仍然需要从 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"])

score 0 · Accepted Answer

这是基本思想（将其作为嵌套循环在地块上然后在土壤上进行；我不确定它是否可以按照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.

r - 如何执行两个 SpatialPolygonsDataFrame 对象的矢量叠加？

4 回答 4

Related

Reference