Znajdź najlepsze dopasowanie nakładających się wielokątów w R

Mam dwa shapefiles, które zostały odczytane do R przy użyciu readOGR() jako obiekty SpatialPolygonsDataFrame. Obie są mapami Nowej Zelandii z różnymi granicami wewnętrznymi. Jeden ma około 70 wielokątów reprezentujących granice terytorialne; drugi ma około 1900 reprezentujących jednostki powierzchni.Znajdź najlepsze dopasowanie nakładających się wielokątów w R

Moim celem - irytująco podstawowym elementem większego projektu - jest użycie tych map do stworzenia tabeli referencyjnej, która może wyszukać jednostkę powierzchni i zwrócić, w którym jest ona w większości terytorialna. Mogę użyć funkcji over() znaleźć, które wielokąty nakładają się, ale w wielu przypadkach jednostki terytorialne wydają się być, przynajmniej w małej części, w obrębie wielu władz terytorialnych - nawet jeśli rozpatrywanie indywidualnych przypadków sugeruje, że normalnie 90% + jednostki powierzchniowej znajduje się w jednym organie terytorialnym.

Czy jest gotowy do użycia oznacza, że robi to, co ponad(), ale które mogą zidentyfikować nie tylko (lub nawet nie) wszystkie nakładające się wielokąty, ale który z kilku nakładających się wielokątów jest najbardziej pokrywające się w każdym przypadku?

Źródło

2013-01-08 Peter Ellis

myślę co chcesz została pokryta na geograficznych systemów informacyjnych SE:

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

W szczególności jeśli wielokąty terytorium są T1, T2, T3 itp a wielokąt próbujesz klasyfikować jest A, prawdopodobnie chcesz użyć gArea na gIntersection z A i T1, następnie A i T2, następnie A i T3 itd., A następnie wybierz, który ma maksymalny obszar. (Czy trzeba pakiet rgeos.)

Źródło

2013-01-08 04:32:23 Silverfish

Dzięki - gArea i gIntersection razem tworzą brakujące ogniwo dla mnie. Wygląda na to, że powinno działać. Jeśli tak, zaakceptuję to jako odpowiedź. –

Miło to słyszeć. Jeśli to zrobisz, byłoby świetnie, gdybyś mógł zostawić działającą próbkę kodu, ponieważ nie mogę uwierzyć, że jesteś jedyną osobą, która stara się wykonać tak oczywiste ważne zadanie i stara się znaleźć odpowiednie narzędzia! – Silverfish

Dzięki @Silverfish - dodałem odpowiedź, która spełnia tę funkcję i powinna być dostosowana do innych. –

Oto kod, który spełnił swoje zadanie, opierając się na @ odpowiedź Silverfish za

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([email protected], 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 <- [email protected][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, [email protected]$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()

enter image description here

Źródło

2013-01-09 05:58:13

Znajdź najlepsze dopasowanie nakładających się wielokątów w R

Odpowiedz

Powiązane problemy