Trouver les meilleurs polygones qui se chevauchent dans R

Question

J'ai deux shapefiles que j'ai lus dans R en utilisant readOGR() comme objets SpatialPolygonsDataFrame. Les deux sont des cartes de la Nouvelle-Zélande avec différentes limites internes. On a environ 70 polygones représentant les limites de l'autorité territoriale; l'autre a environ 1900 unités de surface. Mon but - une partie de base énervante d'un plus grand projet - est d'utiliser ces cartes pour produire une table de référence qui peut rechercher une unité de surface et retourner à quelle autorité territoriale elle appartient le plus. Je peux utiliser over() pour trouver les chevauchements des polygones, mais dans de nombreux cas, les unités territoriales semblent être, au moins en partie, au sein de plusieurs autorités territoriales - même si l'examen de cas individuels suggère que normalement 90% + d'une unité territoriale est dans une seule autorité territoriale.

Y a-t-il un prêt-à-faire qui fait quoi over() mais qui peut identifier non seulement tous les polygones qui se chevauchent, mais lequel des plusieurs polygones qui se chevauchent est le plus chevauchant dans chaque cas?

Answer 1

Je pense que ce que vous voulez a été couvert sur les systèmes d'information géographique SE:

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

En particulier si vos polygones de territoire sont T1, T2, T3 etc et polygone que vous essayez de classer est A, vous voulez probablement utiliser gArea sur le gIntersection de A et T1, puis A et T2, puis A et T3, etc, puis choisissez qui a la surface maximale. (Vous avez besoin du paquet rgeos.)

Answer 2

Voici le code qui a fait le travail, le dessin sur @ réponse de 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()