I have two GIS layers -- call them Soils
and Parcels
-- stored as SpatialPolygonsDataFrame
s (SPDF
s), and I would like to "overlay" them, in the sense described here.
The result of the overlay operation should be a new SPDF in which:
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).The
data.frame
component records the attributes of theSoils
andParcels
polygons within which each atomic polygon falls.
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...
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)))
})
%over%
oroverlay
? See also gis.SE and @Spacedman's cheat sheet: maths.lancs.ac.uk/~rowlings/Teaching/UseR2012/cheatsheet.html – Ari B. Friedmanover()
doesn't do what I'm after. Like I said, I need a function that returns the individual polygons, not just an indicator of whether they overlap. (Plus, even if I can construct the atomic polygonsover
et al. don't help, becauses they count polygons sharing a border as being "over" one another.rgeos::gRelate()
is better in that respect.) – Josh O'BriengIntersection
andgSymdifference
is the closest I got before giving up and posting, but it doesn't get the atomic polygons (i.e. the 10 polygons in the bottom figure that containing no lines in them). This seems a pretty fundamental GIS operation, and I'll be happy to give a bounty to whoever shows me how to do it in R (or with a call to some other easily linked GIS)... – Josh O'Brien