Converting alpha hulls to spatial objects

Alpha shapes to polygons

In order to use α-shapes with spatial data in R, I wanted to convert these shapes to polygons. To accomplish this, I modified some of the code from an RPubs by Barry Rowlingson to create the following function:

library(igraph)
## 
## Attaching package: 'igraph'
## The following objects are masked from 'package:stats':
## 
##     decompose, spectrum
## The following object is masked from 'package:base':
## 
##     union

ashape2poly <- function(ashape){
  # Convert node numbers into characters
  ashape$edges[,1] <- as.character(ashape$edges[,1])
  ashape_graph <- graph_from_edgelist(ashape$edges[,1:2], directed = FALSE)
  if (!is.connected(ashape_graph)) {
    stop("Graph not connected")
  }
  if (any(degree(ashape_graph) != 2)) {
    stop("Graph not circular")
  }
  if (clusters(ashape_graph)$no > 1) {
    stop("Graph composed of more than one circle")
  }
  # Delete one edge to create a chain
  cut_graph <- ashape_graph - E(ashape_graph)[1]
  # Find chain end points
  ends = names(which(degree(cut_graph) == 1))
  path = get.shortest.paths(cut_graph, ends[1], ends[2])$vpath[[1]]
  # this is an index into the points
  pathX = as.numeric(V(ashape_graph)[path]$name)
  # join the ends
  pathX = c(pathX, pathX[1])
  return(pathX)
}

For the reasoning behind the function, check out the RPubs I referred to for guidance. For a sanity check, we can compare the resulting shape to the original α-shape we were trying to replicate:

alphapoly_1 <- ashape2poly(alphashape_1)
## Warning: `is.connected()` was deprecated in igraph 2.0.0.
## ℹ Please use `is_connected()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: `clusters()` was deprecated in igraph 2.0.0.
## ℹ Please use `components()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: `get.shortest.paths()` was deprecated in igraph 2.0.0.
## ℹ Please use `shortest_paths()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

plot(iris_sepals, pch = 19, col = "darkseagreen")
# show the original alpha shape
plot(alphashape_1, lwd = 5, col = "gray", add = TRUE)
# plot the new polygon
lines(iris_sepals[alphapoly_1, ], col = "magenta")

Arcs to lines

To deal with this curvature, I wrote the following function to convert the arcs between points in the hull to a series of very short line segments in order to approximate the curve.

# function to convert an arc into line segments
# given the center of the arc, the radius, the vector, and the angle (radians)
arc2line <- function(center, r, vector, theta, npoints = 100) {
  # Get the angles at the extremes of the arcs
  angles <- anglesArc(vector, theta)
  # Generate sequence of angles along the arc to determine the points
  seqang <- seq(angles[1], angles[2], length = npoints)
  # Generate x coordinates for points along the arc
  x <- center[1] + r * cos(seqang)
  # Generate y coordinates for points along the arc
  y <- center[2] + r * sin(seqang)
  coords.xy <- cbind(x,y)
  line <- Line(coords = coords.xy)
  return(line)
}

Hulls to lines

Using the previous function, I wrote another function to take an α-hull and convert it into a set of SpatialLines objects. The function uses the arc2line function from above to convert each arc in the α-hull into a series of lines, before adding each of these sets of lines together. (This function was updated on 6 January 2021 based on my response to kostas_k84 – thank you to everyone who brought this up to me!)

ahull2lines <- function(hull){
  arclist <- hull$arcs
  lines <- list()
  for (i in 1:nrow(arclist)) {
    # Extract the attributes of arc i
    center_i <- arclist[i, 1:2]
    radius_i <- arclist[i, 3]
    vector_i <- arclist[i, 4:5]
    theta_i <- arclist[i, 6]
    # Convert arc i into a Line object
    line_i <- arc2line(center = center_i, r = radius_i, vector = vector_i, theta = theta_i)
    list_length <- length(lines)
    if(list_length > 0){
      # If a line has already been added to the list of lines
      # Define last_line_coords as the coordinates of the last line added to the list before the ith line
      last_line_coords <- lines[[list_length]]@coords
    }
    if(i == 1){
      # Add the first line to the list of lines
      lines[[i]] <- line_i
    } else if(isTRUE(all.equal(line_i@coords[1,], last_line_coords[nrow(last_line_coords),]))){
      # If the first coordinate in the ith line is equal to the last coordinate in the previous line
      # then those lines should be connected
      # Row bind the coordinates for the ith line to the coordinates of the previous line in the list
      lines[[list_length]]@coords <- rbind(last_line_coords, line_i@coords[2:nrow(line_i@coords),])
    } else {
      # If the first coordinate in the ith line does not match the last coordinate in the previous line
      # then the ith line represents a new line
      # Add the ith line to the list as a new element
      lines[[length(lines) + 1]] <- line_i
    }
  }
  # Convert the list of lines to a Line object
  lines <- Lines(lines, ID = 'l')
  # Convert the Line object to a SpatialLines object
  sp_lines <- SpatialLines(list(lines))
  return(sp_lines)
}

The results look like this:

lines_1 <- ahull2lines(alphahull_1)
# the result is a SpatialLines object
class(lines_1)
## [1] "SpatialLines"
## attr(,"package")
## [1] "sp"

plot(iris_sepals, pch = 19, col = "darkseagreen")
# show the original alpha shape
plot(alphahull_1, lwd = 5, col = "gray", add = TRUE)
# plot the new polygon
plot(lines_1, col = "magenta", add = TRUE)

The resulting SpatialLines object is an almost spot-on approximation of the original α-hull (shown in gray).

SpatialLines to SpatialPolygon

Now, I needed a way to convert the SpatialLines object into a SpatialPolygon that would cover the same shape as the original α-hull. To accomplish this, I wrote a function that takes a SpatialLines object, checks which lines are part of polygons (i.e. form closed shapes), and converts those polygons to a SpatialPolygon.

spLines2poly <- function(sp_lines){
  # Extract the lines slot
  lines_slot <- sp_lines@lines[[1]]
  # Create a list of booleans indicating whether a given Line represents a polygon
  poly_bool <- sapply(lines_slot@Lines, function(x){
    coords <- lines_slot@Lines[[1]]@coords
    # Check if the first coordinate in the line is the same as the last
    all.equal(coords[1,], coords[nrow(coords),])
  })
  # Pull out the lines that form polygons
  poly_lines <- sp_lines[poly_bool]
  poly_lines_slot <- poly_lines@lines
  # Create SpatialPolygons
  sp_polys <- SpatialPolygons(list(Polygons(lapply(poly_lines_slot, function(x) {
    Polygon(slot(slot(x, "Lines")[[1]], "coords"))
  }), ID = "1")))
  return(sp_polys)
}

We can apply this function to the lines_1 object we made from the original alphahull_1:

SpPoly_1 <- spLines2poly(lines_1)
class(SpPoly_1)
## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"

plot(iris_sepals, pch = 19, col = "darkseagreen")
# show the original alpha shape
plot(alphahull_1, lwd = 5, col = "gray", add = TRUE)
# plot the new polygon
plot(SpPoly_1, border = "magenta", add = TRUE)

Alpha hulls to SpatialPolygons

Finally, we can string all of those functions together to create a single function that will convert an α-hull directly into a SpatialPolygon:

ahull2poly <- function(hull){
  # Convert the alpha hull to SpatialLines
  hull2SpatialLines <- ahull2lines(hull)
  # Convert SpatialLines to SpatialPolygon
  SpatialLines2SpatialPolygon <- spLines2poly(hull2SpatialLines)
  return(SpatialLines2SpatialPolygon)
}

As a final sanity check, we can see that the resulting shape is the same as the original shape produced by the alphahull package:

hullpoly_1 <- ahull2poly(alphahull_1)
class(hullpoly_1)
## [1] "SpatialPolygons"
## attr(,"package")
## [1] "sp"

plot(iris_sepals, pch = 19, col = "darkseagreen")
# show the original alpha shape
plot(alphahull_1, lwd = 5, col = "gray", add = TRUE)
# plot the new polygon
plot(hullpoly_1, border = "magenta", add = TRUE)