将（alpha）外壳转换为空间多边形

Question

在 R 中，我希望将围绕一堆点的 alpha 形状多边形转换为一个空间多边形对象。

library(sf)
library(alphahull)

首先，我创建点随机点分布

dat <- matrix(c(1,2,3,4,5, 3,3,5,6,9), ncol = 2)

我找到了覆盖点的 alpha 形状（即包含所有点的多边形）。我对这个函数特别感兴趣，因为它具有根据给定的 alpha

找到或多或少紧密的多边形形状的功能

dat.ashape<- ashape(dat, alpha= 7)

我取极点坐标

coords<- dat.ashape$x[dat.ashape$alpha.extreme,]

我让最后一点和第一个点一样（有一个封闭的形状）

coords<- rbind(coords, coords[1,])

为了让事情正常进行，我需要按顺序对点进行排序

coords<- cbind(coords, NA) 
coords[,3]<- c(1, 5, 3, 2, 4, 6) 
coords<- coords[order(coords[,3]),]

我从坐标矩阵创建简单的空间点特征

dat.sf <- st_multipoint(coords, dim = "XYZ")

...并创建多边形

tst<- dat.sf %>% # 
  st_cast('POLYGON')

最后，比较点和形状分布以及多边形，我能够正确构建多边形，但是这对于六个点来说相当容易！（因为我自己手动设置了正确的顺序）

plot(dat.ashape) 
plot(tst, add=T, col=adjustcolor('red', alpha.f=.3), border=2)

在一个更复杂的例子中，比如说 100 个点，我被困在我应该正确排列点的部分，在 st_cast 进入多边形之前。

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)

.....我显然没有完成技巧。

非常感谢您的帮助！

Answer 1

寻找包 alphahull、rgeos、sf 的进一步替代方案来计算围绕一堆点的外壳，我终于找到了 concaveman 感谢 this post，它可以解决问题，与 sf 对象兼容。

library(concaveman)
library(sf)

set.seed(1)
dat <- matrix(stats::rnorm(100), ncol = 2)
dat <- as.data.frame(dat)
names(dat)<- c('x', 'y')
dat.sf<-st_as_sf(dat, coords=c("x","y"))

polygon <- concaveman(dat.sf)
plot(dat.sf, pch=16)
plot(polygon, add=T, col=adjustcolor('red', alpha.f=.3), col.line='red', border=2)

Answer 2

好吧，我对 concaveman 不满意。我真的很想要 Delaunay triangulation as basis of my hull computation as I like alphahull a lot. Also, after reading this 我想找到一种（或我的）可行方法将从 alphahull 包中检索到的船体转换为空间多边形，我可以进一步将其用于更广泛的空间分析。因此我编写了以下函数来完成这项工作：

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)}

所以，如果你想尝试一下 1000 点云：

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)

我希望它对某人有用。

将（alpha）外壳转换为空间多边形

Converting (alpha) hulls to spatial polygon

r

polygon

convex-hull

sf