Find best matching overlapping polygons in R

Question

I have two shapefiles that I have read into R using readOGR() as SpatialPolygonsDataFrame objects. Both are maps of New Zealand with different internal boundaries. One has about 70 polygons representing territorial authority boundaries; the other has about 1900 representing area units.

My aim - an annoyingly basic part of a bigger project - is to use these maps to produce a reference table that can look up an area unit and return which territorial authority it is mostly in. I can use over() to find the which polygons overlap, but in many cases area units seem to be, at least in small part, within multiple territorial authorities - even though looking at individual cases suggests that normally 90%+ of an area unit is in a single territorial authority.

Is there a ready to hand means that does what over() does but which can identify not just (or not even) all the overlapping polygons, but which of the several overlapping polygons is the most overlapping in each case?

score 7 · Answer 1 · answered Jan 09 '13 at 05:58

Here is code that did the job, drawing on @Silverfish's answer

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

enter image description here

Seem like the core of the answer is in your `j` loop. By iterating through one geometry and calculating the area of each of its features intersection with a second geometry, you can then select features from the second geometry match a specific criterion. — Sean McKenzie, Sep 01 '21 at 06:35
Interesting to see this very old question bounced back up into my feed @SeanMcKenzie . These days this is handled much much simpler with `sf` - haven't checked on the original usecase but I think it would be just `st_join(..., largest = TRUE)`. Would run in seconds and be much easier to write. R really has become a geocomputational powerhouse since I first tackled this sort of question. — Peter Ellis, Sep 06 '21 at 00:50

score 5 · Accepted Answer · edited Apr 13 '17 at 12:33

5

I think what you want has been covered on the geographic information systems 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 particular if your territory polygons are T1, T2, T3 etc and your polygon you are trying to classify is A, probably want to use gArea on the gIntersection of A and T1, then A and T2, then A and T3, etc, then pick which has the maximum area. (You would need the rgeos package.)

edited Apr 13 '17 at 12:33

Community

1
1

answered Jan 08 '13 at 04:32

Silverfish

1,812
1
22
30

2

Thanks - gArea and gIntersection together formed the missing link for me. This looks like it should work. If it does I'll accept this as the answer. – Peter Ellis Jan 08 '13 at 07:10
1

Glad to hear it. If you get it working would be great if you could leave a working code sample as I can't believe you are the only person seeking to perform such an obviously important task and struggling to find the tools to do it! – Silverfish Jan 08 '13 at 13:44
Thanks @Silverfish - I have added an answer that does the job and should be adaptable for others. – Peter Ellis Jan 09 '13 at 05:59

Find best matching overlapping polygons in R

2 Answers2