Alpha Shapes to Polygons

This document shows how to simplify a complex region using the alphahull package.

The resulting simplification will include all the original polygon.

With sincerest apologies to Slartibartfast, we'll get a bit of Norway:

library(alphahull)

## Loading required package: tripack

## Loading required package: sgeostat

## Loading required package: splancs

## Loading required package: sp

## Spatial Point Pattern Analysis Code in S-Plus
## 
## Version 2 - Spatial and Space-Time analysis

library(maps)

## Attaching package: 'maps'

## The following object(s) are masked from 'package:sgeostat':
## 
## in.polygon

norway = map("world", "Norway", exact = TRUE)

Now a bit of rearranging, and we'll remove the NA points, and remove duplicate points.

norway = data.frame(x = norway$x, y = norway$y)
norway = subset(norway, !is.na(x))
norway = norway[!duplicated(paste(norway$x, norway$y)), ]

Now we use the ashape function to simplify the coastline:

norwaySimple = ashape(norway, alpha = 2)
plot(norwaySimple)

The alpha parameter controls the level of simplification, use with care:

par(mfrow = c(2, 2))
plot(ashape(norway, alpha = 0.1))
plot(ashape(norway, alpha = 1))
plot(ashape(norway, alpha = 10))
plot(ashape(norway, alpha = 100))

Note that the alpha-shape always includes all the input points - in the cases with small alpha the set degenerates into disconnected points. If you have a set of data points within your polygon then you could add them to the polygon points and do the alpha shape calculation on that combined set of points.

The return structure of ashape doesn't make it easy to get the polygon out. The plotting function just draws segments in whatever order, so if we try and draw one continuous polygon by indexing via the extremes points, we get a mess:

n10 = ashape(norway, alpha = 10)
plot(n10)
lines(norway[n10$alpha.extremes, ], col = "red")

We can re-order the points using the igraph package. The ashape function returns an object that has the indexes of the start and end point of each segment. We can feed these in as an edge list for a graph object:

library(igraph)

## Attaching package: 'igraph'

## The following object(s) are masked from 'package:tripack':
## 
## convex.hull

n10g = graph.edgelist(cbind(as.character(n10$edges[, "ind1"]), as.character(n10$edges[, 
    "ind2"])), directed = FALSE)
plot(n10g)

This is a circular graph - if we've created a single enclosing polygon. We should check this before proceeding with a few tests:

if (!is.connected(n10g)) {
    stop("Graph not connected")
}
if (any(degree(n10g) != 2)) {
    stop("Graph not circular")
}
if (clusters(n10g)$no > 1) {
    stop("Graph composed of more than one circle")
}

In the next step, we delete one edge of the circle to create a chain, and find the path from one end of the chain to the other.

cutg = n10g - E(n10g)[1]
# find chain end points
ends = names(which(degree(cutg) == 1))
path = get.shortest.paths(cutg, ends[1], ends[2])[[1]]
# this is an index into the points
pathX = as.numeric(V(n10g)[path]$name)
# join the ends
pathX = c(pathX, pathX[1])
# now show the alpha shape plot with our poly on top
plot(n10, lwd = 10, col = "gray")
# get the points from the ashape object
lines(n10$x[pathX, ], lwd = 2)

It should then be a simple matter to turn this into a SpatialPolygons object.