Calculate mean line

This document will show how to calculate mean line from a set of polygons.

The approach is straightforward: first, a central point is found. Then, lines from this central point are drawn in a number of directions. In every direction, distances of intersection with boundaries of polygon are calculated. Mean line then lies in mean distance from the central point.

library(maptools)

## Loading required package: sp
## Checking rgeos availability: TRUE

library(rgeos)

## rgeos version: 0.3-4, (SVN revision 438)
##  GEOS runtime version: 3.3.8-CAPI-1.7.8 
##  Polygon checking: TRUE

First, a set of ellipses is generated. Centres, lengths of minor and major semi-axis, and angle of major axis are randomized.

# function to generate ellipse
ellipse <- function(xc, yc, a, b, alpha) {
    alpha <- alpha * (pi/180)
    t <- seq(0, 2 * pi, length.out = 100)
    x <- xc + a * cos(t) * cos(alpha) - b * sin(t) * sin(alpha)
    y <- yc + a * cos(t) * sin(alpha) + b * sin(t) * cos(alpha)
    return(data.frame(x, y))
}

# randomly generate 100 ellipses, transform them to SpatialPolygon class
lst <- c()
for (i in 1:100) {
    center <- rnorm(n = 2, mean = 0, sd = 1)
    ab <- c(rnorm(n = 1, mean = 10, sd = 0.5), rnorm(n = 1, mean = 3, sd = 0.2))
    alpha <- rnorm(n = 1, mean = 45, sd = 7)
    e <- ellipse(xc = center[1], yc = center[2], a = ab[1], b = ab[2], alpha = alpha)
    e[nrow(e), ] <- e[1, ]
    polygon <- SpatialPolygons(list(Polygons(list(Polygon(e)), ID = as.character(i))))
    lst <- c(lst, polygon@polygons[[1]])
    lst[[i]]@ID <- as.character(i)
}

# plot
spp <- SpatialPolygons(Srl = lst)
plot(spp)
plot(gCentroid(spp), add = T, pch = 16, col = "red")

Now, calculate intersection of polygons with lines drawn in 180 different direction from the center (distributed each 2 degrees along the circle)

#### calculate intersections
alfa <- seq(0, 360 * pi/180, length.out = 180)
len <- 100  # length of line
cross <- which(gIntersects(spgeom1 = gCentroid(spp), spp, byid = T))
selection <- spp[cross, ]
centroid <- c(0, 0)
out <- c()
# for each direction
for (i in alfa) {
    line <- rbind(centroid, c(centroid[1] + cos(i) * len, centroid[2] + sin(i) * 
        len))
    sl <- SpatialLines(list(Lines(slinelist = list(Line(coords = line)), ID = "1")))
    lens <- c()
    # for each polygon
    for (j in 1:length(selection)) {
        int <- gIntersection(sl, selection[j, ])
        lens <- c(lens, gLength(int))
    }
    # keep only summary statistics for every direction
    out <- rbind(out, c(i/pi * 180, quantile(lens, c(0.05, 0.25, 0.5, 0.75, 
        0.95)), mean(lens), sd(lens)))
}

# names of columns of data frame
colnames(out) <- c("angle", "q05", "q25", "q50", "q75", "q95", "mean", "sd")
out <- data.frame(out)
head(out)

##    angle   q05   q25   q50   q75   q95  mean    sd
## 1  0.000 1.744 3.385 4.272 5.027 6.467 4.244 1.444
## 2  2.011 1.824 3.500 4.382 5.228 6.667 4.374 1.489
## 3  4.022 1.890 3.628 4.503 5.438 6.882 4.516 1.537
## 4  6.034 1.952 3.765 4.653 5.627 7.117 4.670 1.589
## 5  8.045 2.021 3.897 4.818 5.904 7.370 4.838 1.644
## 6 10.056 2.096 4.021 5.007 6.158 7.644 5.020 1.703

Finaly, create new spatial objects from calculated statistics.

# calculate position of points, based on the angle and distance from the
# centroid, as stored in variable out
m <- c()  # mean
lower <- c()
upper <- c()
for (i in 1:nrow(out)) {
    len <- out$mean[i]
    sdev <- out$sd[i]
    angle <- out$angle[i] * pi/180
    point.m <- c(centroid[1] + cos(angle) * len, centroid[2] + sin(angle) * 
        len)
    point.lower <- c(centroid[1] + cos(angle) * (len - 1.96 * sdev), centroid[2] + 
        sin(angle) * (len - 1.96 * sdev))
    point.upper <- c(centroid[1] + cos(angle) * (len + 1.96 * sdev), centroid[2] + 
        sin(angle) * (len + 1.96 * sdev))
    m <- rbind(m, point.m)
    lower <- rbind(lower, point.lower)
    upper <- rbind(upper, point.upper)
}

sl.m <- SpatialLines(list(Lines(slinelist = list(Line(coords = m)), ID = "1")))
sl.lower <- SpatialLines(list(Lines(slinelist = list(Line(coords = lower)), 
    ID = "1")))
sl.upper <- SpatialLines(list(Lines(slinelist = list(Line(coords = upper)), 
    ID = "1")))


plot(spp, border = "grey80")
plot(SpatialPoints(matrix(c(0, 0), ncol = 2)), pch = 19, col = "black", add = T, 
    cex = 1.4)
plot(sl.m, add = T, col = "black", lwd = 2)
plot(sl.lower, add = T, col = "grey40", lwd = 2, lty = 2)
plot(sl.upper, add = T, col = "grey40", lwd = 2, lty = 2)