Example of calculating and plotting ellipse overlap with phonR, sp, and rgeos

vowel-ellipse-overlap.png

vowel-ellipse-overlap.R

library(phonR)
data(indoVowels)

library(sp)     # SpatialPolygons, etc
library(rgeos)  # gIntersection

fem2 <- indo[indo$subj %in% "F02",]
vowels <- unique(as.character(fem2$vowel))
## important! ellipse.conf defaults to 0.6827 (2 s.d.) for plotVowels, but
## defaults internally to 0.95 for phonR:::ellipse.  You must specify the
## same value for both in order for calculated overlap areas to match plotted
## overlap areas.
conf <- 0.95

## get mean, covariance, and number of data points for each vowel
means <- by(fem2[c("f2", "f1")], fem2["vowel"], colMeans)
covars <- by(fem2[c("f2", "f1")], fem2["vowel"], cov)
n_pts <- by(fem2, fem2["vowel"], nrow)

## get points defining the edge of each vowel ellipse
ellipses <- list()
for (v in vowels) {
    ellipses[[v]] <- phonR:::ellipse(means[[v]], covars[[v]],
                                     n_pts[[v]], alpha=1-conf,
                                     draw=FALSE)
    # make sure the first and last points coincide exactly
    ellipses[[v]][1,] <- ellipses[[v]][nrow(ellipses[[v]]),]
}

## convert ellipse points into spatial polygons objects
pgon <- sapply(ellipses, Polygon, hole=FALSE)
pgons <- sapply(names(pgon), function(i) Polygons(pgon[i], ID=i))
spgons <- sapply(names(pgons), function(i) SpatialPolygons(pgons[i]))

## get overlap & total area for each pair of spatial polygons objects
## (use gUnion to avoid counting overlap area twice when calculating
## total area)
cmbns <- combn(vowels, 2)
overlap <- apply(cmbns, 2,
                 function(i) gIntersection(spgons[[i[1]]], spgons[[i[2]]]))
union <- apply(cmbns, 2, 
               function(i) gUnion(spgons[[i[1]]], spgons[[i[2]]]))
names(overlap) <- names(union) <- apply(cmbns, 2, paste, collapse="-")
ol_area <- sapply(overlap,
                  function(i) if (is.null(i)) {0} else {i@polygons[[1]]@area})
un_area <- sapply(union,
                  function(i) if (is.null(i)) {0} else {i@polygons[[1]]@area})
pct_ol_area <- apply(rbind(ol_area, un_area), 2,
                     function(i) 100 * i[1] / i[2])

## optional: plot the overlap
with(fem2, plotVowels(f1, f2, vowel, ellipse.line=TRUE, ellipse.conf=conf))
invisible(lapply(overlap, function(i) if (is.null(i)) {NULL}
                 else {polygon(i@polygons[[1]]@Polygons[[1]]@coords,
                               border=NA, col="#FF000066")}))

Here is the content of ol_area, un_area, and pct_ol_area

cbind(ol_area, un_area, pct_ol_area)
#       ol_area  un_area pct_ol_area
# a-e 168377.37 670452.3   25.113997
# a-i      0.00 709137.7    0.000000
# a-o  22135.82 694426.1    3.187642
# a-u      0.00 662454.8    0.000000
# e-i  35916.01 566342.2    6.341750
# e-o      0.00 609682.4    0.000000
# e-u      0.00 555575.2    0.000000
# i-o      0.00 479990.5    0.000000
# i-u      0.00 425883.3    0.000000
# o-u  85373.52 347934.0   24.537275