In R, I wish to convert the alpha shape polygon surrounding a bunch of points into one single spatial polygon object.
library(sf)
library(alphahull)
To start out, I create the point random points distribution
dat <- matrix(c(1,2,3,4,5, 3,3,5,6,9), ncol = 2)
I find the alpha shape covering the points (i.e. a polygon encompassing all points). I am particularly interested in this function as it has the feature to find a more or less tight polygon shape according to the given alpha
dat.ashape<- ashape(dat, alpha= 7)
I take the coordinates of the extreme
coords<- dat.ashape$x[dat.ashape$alpha.extreme,]
I make the last point same as the first (to have a closed shape)
coords<- rbind(coords, coords[1,])
To make things to work I need to order the point in sequence
coords<- cbind(coords, NA)
coords[,3]<- c(1, 5, 3, 2, 4, 6)
coords<- coords[order(coords[,3]),]
I create the simple spatial point feature from the coordinate matrix
dat.sf <- st_multipoint(coords, dim = "XYZ")
... and create the polygon
tst<- dat.sf %>% #
st_cast('POLYGON')
Finally, comparing the point and shape distribution and the polygon, I was able to build the polygon correctly, but this is rather easy with six points! (Because I made myself manually the right order)
plot(dat.ashape)
plot(tst, add=T, col=adjustcolor('red', alpha.f=.3), border=2)
In a more sophisticated example with say 100 points, I get stuck in the part where I should get the sequence of points right, before st_cast into polygon.
set.seed(1)
dat <- matrix(stats::rnorm(100), ncol = 2)
dat.ashape<- ashape(dat, alpha=7)
coords<- dat.ashape$x[dat.ashape$alpha.extreme,]
coords<- rbind(coords, coords[1,])
dat.sf <- st_multipoint(coords, dim = "XY")
tst <- dat.sf %>%
st_cast('POLYGON')
plot(dat.ashape)
plot(tst, add=T, col=adjustcolor('red', alpha.f=.3), col.line='red', border=2)
.... and I obviously do not get trick done.
I am grateful for any help!
OK, I was not happy with the concaveman. I really wanted the Delaunay triangulation as basis of my hull computation as I like alphahull a lot. Also, after reading this I wanted to find a (or my) viable way for converting the hull retrieved from alphahull package to a spatial polygon, which I could further use for my broader spatial analysis. Therefore I wrote the following function to do the job:
hull2poly <- function(my.ashape){
require(sf)
if(class(my.ashape) != "ashape") {stop('error, your input must be
ashape class')} else
my.edge<- data.frame(my.ashape$edges)[,c( 'x1', 'y1', 'x2', 'y2')]
x<- my.edge[,1:2]
y<- my.edge[,3:4]
my.edge2<- matrix(t(cbind(x,y)), byrow=T,ncol=2)
my.edge2<- as.data.frame(my.edge2)
names(my.edge2)<- c('x','y')
my.edge2$id <- unlist(lapply((1: (nrow(my.edge2)/2)),
FUN=function(x){c(rep(x,2))}))
start.edge<- 1
new.id<- start.edge
new.edges<- my.edge2[which(my.edge2$id== start.edge ),]
while(length(new.id)<= length(unique(my.edge2$id))-1){
internal.id<- new.id[length(new.id)]
edge <- my.edge2[which(my.edge2$id== internal.id ),]
where.to.search <- my.edge2[which(my.edge2$id %in% new.id ==F ),]
index1<- apply(where.to.search[,1:2], 1, function(x){x == edge[1,1:2]})
index1<- as.numeric(names(which(apply(index1,2, sum)>0)))[1]
index2<- apply(where.to.search[,1:2], 1, function(x){x == edge[2,1:2]})
index2<- as.numeric(names(which(apply(index2,2, sum)>0)))[1]
main.index<- c(index1, index2)
ifelse(all(!is.na(main.index)),
# yes
{flag<- c(T,T)
main.index<- main.index[2]
point.coord<- my.edge2[main.index,]
segment<- my.edge2[my.edge2$id==my.edge2[main.index,'id'],]
new.id<- c( new.id, my.edge2[main.index,]$id) },
# no
ifelse(which(!is.na(main.index))==1,
# yes
{flag<- c(T,F)
main.index<- main.index[flag]
point.coord<- my.edge2[main.index,]
segment<-
my.edge2[my.edge2$id==my.edge2[main.index,'id'],]
new.id<- c( new.id, my.edge2[main.index,]$id)},
# no
{flag<- c(F,T)
main.index<- main.index[flag]
point.coord<- my.edge2[main.index,]
segment<- my.edge2[my.edge2$id==my.edge2[main.index,'id'],]
new.id<- c( new.id, my.edge2[main.index,]$id)} ) )
index3<- t(apply(segment, 1, function(x){x ==point.coord}))
new.edges<- rbind(new.edges, rbind(point.coord, segment[which(apply(index3,1, sum)<3),]))
}
tst <- st_multipoint(as.matrix(new.edges), dim = "XYZ")
poly<- tst %>% #
st_cast('POLYGON')
return(poly)}
So, if you wish to give a try with a cloud of 1000 points:
library(alphahull)
set.seed(1)
dat <- matrix(stats::rnorm(1000), ncol = 2)
dat <- as.data.frame(dat)
dat.ashape<- ashape(dat, alpha= 2)
tmp<- hull2poly(dat.ashape)
plot(tmp)
I hope, it comes useful for someone.