Come eseguire una sovrapposizione vettoriale di due oggetti SpatialPolygonsDataFrame?

Question

Ho due strati GIS - chiamiamoli Soils e Parcels - memorizzati come SpatialPolygonsDataFrame s (SPDF s), e mi piacerebbe "sovrapporli", in the sense described here.

Il risultato dell'operazione sovrapposizione dovrebbe essere una nuova SPDF in cui:

1. Il SpatialPolygons componente contiene poligoni formati dall'intersezione dei due strati. (Pensa a tutti i poligoni atomici formati sovrapponendo due mylar su una lavagna luminosa).

2. Il componente data.frame registra gli attributi dei poligoni Soils e Parcels all'interno dei quali cade ciascun poligono atomico.

La mia domanda (s): C'è una funzione di R esistente che fa questo? (Sarei anche felice di apprendere una funzione che ottiene il componente SpatialPolygons giusto, calcolando i poligoni atomici formati dall'intersezione dei due strati.) Mi sento come rgeos dovrebbe avere una funzione che fa (1) almeno , ma sembra non ...

Ecco una figura che può aiutare a rendere più chiaro ciò che cerco, seguito dal codice che crea i livelli Soils e Parcels mostrato nella figura.

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

Dal gennaio del 2014, il pacchetto raster include una funzione union() che rende questo facile-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

Ecco l'idea di base (fare questo come un ciclo nidificato su pacchi poi col suolo; io non sono sicuro che possa essere vettorializzare il modo in cui le g* funzioni sono scritte):

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) 

Questo dovrebbe essere sufficiente per iniziare. Il resto è solo in loop, tenendo traccia dei pezzi e unendoli tutti in un unico SPDF.

Un approccio alternativo che è più vettorizzati è:

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

Questo vi dà più pezzi che in realtà esiste, per qualche ragione. Quindi dovrai tornare indietro e raggrupparli e selezionare i pezzi extra che sono gEquals.

Answer 3

Ecco il mio crack, questo dà solo un elenco dei pezzi con i dati Parcel (Parcels-> Soils). È ancora necessario ottenere gli attributi dall'oggetto Suoli, quindi eseguire un lavoro simile con le "differenze" e viceversa (Suoli-> Pacchi) in modo da poter avere qualsiasi tipo di relazioni di sovrapposizione.

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

Dal gennaio 2014, la soluzione postato qui è stato completamente sostituito dalla funzione raster::union(), dimostrata nella risposta ufficialmente accettato sopra.

---

Questo ottiene i poligoni "Atomic" formate dalla sovrapposizione di due diversi strati SpatialPolygons, risolvendo parte 1 della mia domanda.

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) 

E questo estrae i ids (eventuali) dei poligoni nei due strati di input overlain da ciascun poligono "atomica" emessi da gFragment() sopra, risolvendo parte 2 della mia domanda.

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

Edit:

Si noti che questo codice potrebbe non riuscire se uno qualsiasi dei vostri poligoni di ingresso sono di id chiamato "1". In tal caso, una soluzione alternativa è rinominare l'id, magari facendo qualcosa come Parcels <- spChFIDs(Parcels, paste0("pp", row.names(Parcels))).