Trova i migliori poligoni sovrapposti corrispondenti in R

Question

Ho due shapefile che ho letto in R usando readOGR() come oggetti SpatialPolygonsDataFrame. Entrambe sono mappe della Nuova Zelanda con diversi confini interni. Uno ha circa 70 poligoni che rappresentano i confini delle autorità territoriali; l'altro ha circa 1900 unità rappresentative dell'area.

Il mio obiettivo - una parte fastidiosamente fondamentale di un progetto più grande - è utilizzare queste mappe per produrre una tabella di riferimento che può cercare un'area e restituire quale autorità territoriale è principalmente in. Posso usare over() per trova i poligoni che si sovrappongono, ma in molti casi le unità di area sembrano essere, almeno in piccola parte, all'interno di più autorità territoriali - anche se i singoli casi suggeriscono che normalmente il 90% + di un'unità di area si trova in un'unica autorità territoriale.

C'è un mezzo pronto per l'utente che fa ciò che over() fa ma che può identificare non solo (o nemmeno) tutti i poligoni sovrapposti, ma quale dei vari poligoni sovrapposti è il più sovrapposto in ogni caso?

Answer 1

Penso che ciò che si desidera è stato coperto sui sistemi di informazione geografica SE:

https://gis.stackexchange.com/questions/40517/using-r-to-calculate-the-area-of-multiple-polygons-on-a-map-that-intersect-with?rq=1

In particolare se i vostri poligoni territorio sono T1, T2, T3, ecc e poligono si sta tentando di classificare è A, probabilmente si desidera utilizzare gArea su gIntersection di A e T1, quindi A e T2, quindi A e T3, ecc., Quindi selezionare quale ha l'area massima. (Si avrebbe bisogno il pacchetto rgeos.)

Answer 2

Qui è il codice che ha fatto il lavoro, attingendo @ risposta di Silverfish

library(sp) 
library(rgeos) 
library(rgdal) 

### 
# Read in Area Unit (AU) boundaries 
au <- readOGR("C:/Users/Peter Ellis/Documents/NZ", layer="AU12") 

# Read in Territorial Authority (TA) boundaries 
ta <- readOGR("C:/Users/Peter Ellis/Documents/NZ", layer="TA12") 

### 
# First cut - works ok when only one TA per area unit 
x1 <- over(au, ta) 
au_to_ta <- data.frame(au@data, TAid = x1) 

### 
# Second cut - find those with multiple intersections 
# and replace TAid with that with the greatest area. 

x2 <- over(au, ta, returnList=TRUE) 

# This next loop takes around 10 minutes to run: 
for (i in 1:nrow(au_to_ta)){ 
    tmp <- length(x2[[i]]) 
    if (tmp>1){ 
     areas <- numeric(tmp) 
     for (j in 1:tmp){ 
      areas[j] <- gArea(gIntersection(au[i,], ta[x2[[i]][j],])) 
      } 
#  Next line adds some tiny random jittering because 
#  there is a case (Kawerau) which is an exact tie 
#  in intersection area with two TAs - Rotorua and Whakatane 

     areas <- areas * rnorm(tmp,1,0.0001) 

     au_to_ta[i, "TAid"] <- x2[[i]][which(areas==max(areas))] 
    } 

} 


# Add names of TAs 
au_to_ta$TA <- ta@data[au_to_ta$TAid, "NAME"] 

#### 
# Draw map to check came out ok 
png("check NZ maps for TAs.png", 900, 600, res=80) 
par(mfrow=c(1,2), fg="grey") 
plot(ta, col=ta@data$NAME) 

title(main="Original TA boundaries") 
par(fg=NA) 
plot(au, col=au_to_ta$TAid) 
title(main="TA boundaries from aggregated\nArea Unit boundaries") 
dev.off()