如何执行两个SpatialPolygonsDataFrame对象的矢量叠加? [英] How to perform a vector overlay of two SpatialPolygonsDataFrame objects?

问题描述

我有两个GIS层-分别称为Soils和Parcels-存储为SpatialPolygonsDataFrame s(SPDF s)，我想对其进行覆盖"，

I have two GIS layers -- call them Soils and Parcels -- stored as SpatialPolygonsDataFrames (SPDFs), and I would like to "overlay" them, in the sense described here.

叠加操作的结果应为新的SPDF，其中:

The result of the overlay operation should be a new SPDF in which:

1. SpatialPolygons组件包含由两层的相交形成的多边形. (想一想在高射投影仪上重叠两个聚酯薄膜形成的所有原子多边形.)

1. The SpatialPolygons component contains polygons formed by the intersection of the two layers. (Think all of the atomic polygons formed by overlaying two mylars on an overhead projector).

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

The data.frame component records the attributes of the Soils and Parcels polygons within which each atomic polygon falls.

我的问题:是否存在执行此操作的现有R函数? (我什至乐于学习一个仅能正确获取SpatialPolygons组件的函数，计算由两层的交点形成的原子多边形.)我觉得 rgeos 应该具有一个函数至少可以做到(1)，但似乎没有...

My question(s): Is there an existing R function that does this? (I would even by happy to learn of a function that just gets the SpatialPolygons component right, calculating the atomic polygons formed from the intersection of the two layers.) I feel like rgeos should have a function that does (1) at least, but it seems not to...

下面的图可以帮助您更清楚地了解我所追求的内容，然后是创建图中所示的Soils和Parcels层的代码.

Here is a figure that may help make it clearer what I'm after, followed by code that creates the Soils and Parcels layers shown in the figure.

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

推荐答案

自2014年1月起， raster 软件包中包含一个union()函数，该函数使操作变得简单:

Since January of 2014, the raster package includes a union() function that makes this easy-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)

这篇关于如何执行两个SpatialPolygonsDataFrame对象的矢量叠加?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！