Extract multiple polygons from alphahull in R

user_3321151 0
Viewed 321

I'm creating map boundaries with alphahull and result is sometimes discrete hulls (which is good) . . three nice clusters in the example below. I can get the number of discrete clusters using igraph, but I want to close the polygons and I don't see a simple way to assign points to the correct cluster. What am I missing? Eventually I want to pass the object to ggplot as polygons.

library(alphahull)
library(ggplot2)
library(igraph)
#> 
#> Attaching package: 'igraph'
#> The following objects are masked from 'package:stats':
#> 
#>     decompose, spectrum
#> The following object is masked from 'package:base':
#> 
#>     union
cluster_one.df <- data.frame("lon"=rnorm(20,135,0.2),"lat"=rnorm(20,35,0.2)) 
cluster_two.df <- data.frame("lon"=rnorm(20,130,0.2),"lat"=rnorm(20,30,0.2))
cluster_three.df <- data.frame("lon"=rnorm(20,125,0.2),"lat"=rnorm(20,25,0.2))
points.df <- rbind(cluster_one.df,cluster_two.df,cluster_three.df)
plot(points.df)



alpha_obj <- ashape(points.df, alpha=0.5)
find_no <- alpha_obj$edges
plot(alpha_obj)

Alphahull has three clear clusters . . .


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


clusters(network)$no
#> [1] 3
is.connected(network)
#> [1] FALSE

ggplot(points.df, aes(x = lon, y = lat)) +
  geom_point() +
  geom_segment(data=data.frame(alpha_obj$edges), aes(x = x1, y = y1, xend = x2, yend = y2), color="red")

igraph reports the right number of clusters, and the segments are correct in ggplot, but how to join the segments into three polygons?

Created on 2020-05-23 by the reprex package (v0.3.0)

Update: Continuing with igraph, I was able to subset the data by group

alpha.df <- data.frame(alpha_obj$edges)
V(network)$comp <- components(network, "strong")$membership
group_1 <- induced_subgraph(network, V(network)$comp==1)

And borrowing from https://rpubs.com/geospacedman/alphasimple, this get the order of the points

cutg <- group_1 - E(group_1)[1] ## break loop
ends <- names(which(degree(cutg) == 1))
path <-  get.shortest.paths(cutg, ends[1], ends[2])[[1]]

names(path[[1]]) # this gets the sequence for the first group

But, I still need to order the line segments for correct polygons...which seems to lead back to a loop.

The solution by Allan Cameron below works well, but throws this error on this real world data set when alpha is small, such as 0.1. Then two alpha shapes seem to share a line!

enter image description here

Error in m[, j] == m[i, j] : non-conformable arrays In addition: Warning messages: 1: In if (j == 1) xs[i] <- edge.df$x1[i] else xs[i] <- edge.df$x2[i] : the condition has length > 1 and only the first element will be used 2: In if (j == 1) ys[i] <- edge.df$y1[i] else ys[i] <- edge.df$y2[i] : the condition has length > 1 and only the first element will be used

points.df <- 
structure(list(lon = c(132.6654637, 132.628814738644, 132.539254, 
132.5553958, 132.539302181878, 132.7950567, 132.641769, 132.7684343, 
132.7888039, 132.6647853, 132.8347288, 132.8278186, 132.8982418, 
132.909312197459, 132.8625508, 132.82386, 132.809095, 132.7759964, 
132.7756813, 132.7870982, 132.7851847, 132.7681695, 132.8032465, 
132.8160503, 132.6358063, 132.635823981036, 132.5445085, 132.6403534, 
132.5769612, 132.609291, 132.6340415, 132.6381995, 132.6648159, 
132.7101352, 132.7194943, 132.7876312, 132.7825205, 132.7585091, 
132.7559927, 132.7685386, 132.7875584, 132.7872259, 132.7874004, 
132.7902362, 132.8010577, 132.7922278, 132.7861715, 132.768931, 
132.7619977, 132.7513127, 132.7391948, 132.790737970094, 132.699771, 
132.6791708, 132.6854316, 132.6888423, 132.6958107, 132.7075791, 
132.7359368, 132.7008823, 132.7094064, 132.7131504, 132.708409, 
132.7251742, 132.7357202, 132.742711, 132.7600732, 132.7732886, 
132.727009, 132.5486109, 132.5470441, 132.5480779, 132.5586113, 
132.541613, 132.516896, 132.5240977, 132.529158, 132.576234, 
132.5889441, 132.6164585, 132.6411233, 132.6148395, 132.7246135, 
132.6422625, 132.7226471, 132.7363183, 132.7421675, 132.691419, 
132.6641904, 132.665175, 132.661818, 132.6559487, 132.6567813, 
132.6579793, 132.7164382, 132.6956271, 132.6983378, 132.6573281, 
132.6543369, 132.6392342, 132.648159, 132.6758976, 132.6813358, 
132.6640595, 132.7096546, 132.587736, 132.5992522, 132.5865522, 
132.5473471, 132.5728439, 132.5487525, 132.5445001, 132.547059761893, 
132.5538985, 132.556837, 132.5422245, 132.5413181, 132.5451487, 
132.5387885, 132.5249464, 132.513415, 132.4995801, 132.4930108, 
132.5204733, 132.4952279, 132.4331523, 132.4692833, 132.4746544, 
132.5247713, 132.5374035, 132.5444085, 132.5947945, 132.6278865, 
132.7003073, 132.724446383125, 132.7509583, 132.5697286, 132.5031789, 
132.5261544, 132.433172666015, 132.7995, 132.6734, 132.9806, 
132.5196, 132.5454, 132.7716, 132.7272, 132.6548, 132.654852573297, 
132.7008, 132.738540411042), lat = c(33.9564431, 33.9754360356657, 
33.933059, 33.9809751, 33.9329099108545, 33.5110363, 33.378098, 
33.5268887, 33.5509639, 33.5365971, 33.575854, 33.5742345, 33.6157401, 
33.6432153491279, 33.6330822, 33.622899, 33.6338496, 33.7256174, 
33.6549261, 33.6135064, 33.5958084, 33.5802348, 33.5830577, 33.7204944, 
33.6225688, 33.6224810289754, 33.5061936, 33.6027543, 33.651187, 
33.667152, 33.6850872, 33.6879397, 33.70622, 33.6404735, 33.6829081, 
33.7565837, 33.750887, 33.7458008, 33.7331328, 33.7286281, 33.7365333, 
33.7071535, 33.7465765, 33.7467099, 33.7060568, 33.7492625, 33.6732325, 
33.692618, 33.7126021, 33.7205759, 33.7217591, 33.7759301629704, 
33.774027, 33.7372899, 33.7263687, 33.7433889, 33.7534323, 33.7168509, 
33.7228451, 33.754549, 33.7582016, 33.7683532, 33.7761549, 33.7663288, 
33.7690062, 33.7768026, 33.7815469, 33.7863218, 33.801855, 33.5033857, 
33.5442034, 33.4938868, 33.4587104, 33.4522772, 33.4445049, 33.4605537, 
33.4735192, 33.4983482, 33.5126441, 33.5051541, 33.484071, 33.4607175, 
33.5217105, 33.461865, 33.4978878, 33.4847207, 33.5018607, 33.502981, 
33.4882015, 33.521964, 33.511535, 33.5266055, 33.5497139, 33.5328971, 
33.5744565, 33.5492812, 33.5675655, 33.561057, 33.5538189, 33.5993602, 
33.5783269, 33.5976171, 33.6069939, 33.6064045, 33.6441148, 33.6010336, 
33.5800368, 33.523926, 33.516037, 33.533482, 33.5201949, 33.5251998, 
33.5440333757024, 33.5622518, 33.5667931, 33.5634791, 33.593845, 
33.5816383, 33.5713532, 33.579981, 33.5975215, 33.5963463, 33.6126805, 
33.5777769, 33.5926493, 33.5579864, 33.5623084, 33.5889524, 33.5128369, 
33.5181171, 33.5069546, 33.573771, 33.5544736, 33.5249674, 33.5215648578441, 
33.5097935, 33.5181873, 33.5354703, 33.6342814, 33.5579579304105, 
33.64085, 33.53592, 33.62881, 33.44105, 33.46809, 33.49512, 33.57256, 
33.53509, 33.5351158000868, 33.52385, 33.48571)), row.names = c("36855", 
"36856", "36857", "36858", "36859", "36996", "36997", "36998", 
"36999", "37000", "37001", "37002", "37003", "37004", "37005", 
"37006", "37007", "37008", "37009", "37010", "37011", "37012", 
"37013", "37014", "37015", "37016", "37017", "37018", "37019", 
"37020", "37021", "37022", "37023", "37024", "37025", "37026", 
"37027", "37028", "37029", "37030", "37031", "37032", "37033", 
"37034", "37035", "37036", "37037", "37038", "37039", "37040", 
"37041", "37042", "37043", "37044", "37045", "37046", "37047", 
"37048", "37049", "37050", "37051", "37052", "37053", "37054", 
"37055", "37056", "37057", "37058", "37059", "37077", "37078", 
"37079", "37080", "37081", "37082", "37083", "37084", "37085", 
"37086", "37087", "37088", "37089", "37090", "37091", "37092", 
"37093", "37094", "37095", "37096", "37097", "37098", "37099", 
"37100", "37101", "37102", "37103", "37104", "37105", "37106", 
"37107", "37108", "37109", "37110", "37111", "37112", "37113", 
"37114", "37115", "37116", "37117", "37118", "37119", "37120", 
"37121", "37122", "37123", "37124", "37125", "37126", "37127", 
"37128", "37129", "37130", "37131", "37132", "37133", "37134", 
"37135", "37136", "37137", "37138", "37142", "37143", "37144", 
"37145", "37146", "37147", "37148", "37149", "37150", "88530", 
"88537", "94189", "94192", "94193", "94194", "94195", "94196", 
"94197", "94198", "94199"), class = "data.frame")
igraph ggplot2 geospatial r
2 Answers

If there's an easy way to do this, I didn't find it. There is a hard way though: each row of the find_no data frame contains the start and end co-ordinates of an edge. You therefore start at the first row and find which other row shares a co-ordinate with it. You label the first row as being the first point of polygon 1, then move to the row that shares a point with it. You repeat this process until you arrive back at the edge you started, and you have labelled your first polygon. Now repeat starting at the first available unlabelled row. Continue this process until no unlabelled rows are left.

Here's a function that does all that:

extract_polygons <- function(alpha_obj)
{
  if(class(alpha_obj) != "ashape") stop("extract_polygons requires an ashape")
  edge.df <- as.data.frame(alpha_obj$edges)
  groups <- ns <- xs <- ys <- numeric(nrow(edge.df))
  m <- cbind(edge.df[[1]], edge.df[[2]])
  group <- 1
  repeat {
    i <- which(groups == 0)[1]
    if (length(i) == 0 | is.na(i)) break()
    j <- n <- 1
    repeat {
      groups[i] <- group
      ns[i] <- n
      if(j == 1) xs[i] <- edge.df$x1[i] else xs[i] <- edge.df$x2[i]
      if(j == 1) ys[i] <- edge.df$y1[i] else ys[i] <- edge.df$y2[i]
      next_ind <- which((m[, j] == m[i, j] | m[, (j %% 2 + 1)] == m[i, j]) & groups == 0)
      if (length(next_ind) == 0) break()
      j <- which(m[next_ind,] == m[i, j]) %% 2 + 1
      i <- next_ind
      n <- n + 1
    }
    group <- group + 1
  }
  data.frame(x = xs, y = ys, group = as.factor(groups))[order(groups, ns), ]
}

So we can extract the polygons right out of your "ashape" object and use them in a plot.

polygon.df <- extract_polygons(alpha_obj)

ggplot(points.df, aes(lon, lat)) + 
  geom_point() + 
  geom_polygon(data = polygon.df, aes(x, y, fill = group), colour = "black", alpha = 0.5)

enter image description here

Update

On the data set given by the OP, I get the following result with exactly the same code, with no errors:

enter image description here

user_125003150

Here's an alternative solution, wrapped as a function. It has a major dependency (igraph) and also uses rbindlist from data.table, but it works even for tiny alpha settings

alpha_map <- function(lon_lat, alpha=0.5){
  alpha_obj <- ashape(lon_lat, alpha = alpha)

  network <-  graph.edgelist(cbind(as.character(alpha_obj$edges[, "ind1"]), as.character(alpha_obj$edges[, "ind2"])), directed = FALSE)
  group.df <- data.frame("group"=components(network, "strong")$membership)
  group.df$IDs <- names(components(network, "strong")$membership)


  alpha.df <- data.frame(alpha_obj$edges)
  alpha_temp.df <- data.frame("IDs"= c(alpha.df$ind1,alpha.df$ind2), "x"=c(alpha.df$x1,alpha.df$x2),"y"=c(alpha.df$y1,alpha.df$y2))


  ## all three parts as a list
  segments.ls<- lapply(1:clusters(network)$no, function (x) atan_sort(alpha_temp.df[which(alpha_temp.df$IDs%in%group.df$IDs[which(group.df$group==x)]),]))
  # segments.ls<- lapply(1:1, function (x) atan_sort(alpha_temp.df[which(alpha_temp.df$IDs%in%group.df$IDs[which(group.df$group==x)]),]))

  segments.df <- rbindlist(segments.ls)
  ## now assign groups and unlist to data frame

  group_vector <- NULL
  for (i in 1:clusters(network)$no)
  {
    temp <- rep(i,nrow(segments.ls[[i]]))
    group_vector <- c(group_vector,temp)
  }

  segments.df$group <- group_vector
  return(segments.df)
}
user_33211510
Related