I have an alphahull::ashape() result that I need to convert with spatstat::owin() for subsequent spatial analyses. Is there an elegant solution how to order the points from the concave hull and convert the shape to an owin object or any Spatial* object? With ordered points in a hypothetical matrix coords the polygon is defined as:
geo.owin=try(owin(poly=list(x=coords[,1],y=coords[,2])),silent=T)
if(class(geo.owin)=="try-error") geo.owin=owin(poly=list(x=rev(coords[,1]),y=rev(coords[,2])))
This is not possible. The ashape()$edges matrix contains coordinates of points that can be drawn to look like a polygon and the matrix is provided here. However, the points are not ordered as demonstrated with different plotting colors and directions of arrows.
# data
okraj = matrix(c(23.8808, 18.0106, 23.8808, 1.8265, 1.8266, 1.8265, 18.0106, 39.5352, 39.5352, 39.5352, 13.5519, 27.9675, 3.8102, 4.8269, 8.8236, 52.7248, 45.3385, 52.7248, 50.9600, 50.9600, 50.9600, 45.3385, 39.9042, 39.9042, 39.9042, 49.1204, 41.8421, 47.1450, 44.9423, 46.0246, 13.5519, 27.9675, 32.1116, 1.8266, 4.5644, 3.8102, 17.5557, 39.5840, 32.1116, 29.2158, 4.5644, 29.2158, 4.8269, 8.8236, 17.5557, 49.1204, 41.8421, 41.3710, 50.9600, 49.8638, 47.1450, 45.5063, 39.8987, 41.3710, 40.0750, 49.8638), ncol=4, dimnames=list(NULL,c("x1","y1","x2","y2")))
# draw polygon from arrows
farby=colorRampPalette(c("lightblue","black"))(nrow(okraj))
plot(0,type="n",xlim=range(okraj[,c(1,3)]),ylim=range(okraj[,c(2,4)]),xlab="",ylab="")
for(i in 1:nrow(okraj)){
arrows(okraj[i,"x1"],okraj[i,"y1"],okraj[i,"x2"],okraj[i,"y2"],col=farby[i],length=.1,lwd=2)
}
The search through other questions returned solutions for ordered points, convex not concave shapes, solutions in C, python and Lua languages. My current code searches for connected points across the matrix and uses if() conditions for every problem I ran into in interpreting the alpha shape edge coordinates (pardon the naive coding).
# function finding a row, which contains the connecting data
find.connection=function(bod, temp){
# test whether a connection exists in the first set of point coordinates
riadok1=ifelse(length(which(temp[,1]==bod[1]))==0,NA,which(temp[,1]==bod[1]))
riadok2=ifelse(length(which(temp[,2]==bod[2]))==0,NA,which(temp[,2]==bod[2]))
# test for a connection in the second set of coordinates
if(is.na(riadok1)){
riadok1=ifelse(length(which(temp[,3]==bod[1]))==0,NA,which(temp[,3]==bod[1]))
riadok2=ifelse(length(which(temp[,4]==bod[2]))==0,NA,which(temp[,4]==bod[2]))
}
# check multiple values in x or y coordinates and select a row where both occur
if(riadok1==riadok2){
riadok=riadok1
} else {
riadky1=ifelse(length(which(temp[,3]==bod[1]))==0,NA,which(temp[,3]==bod[1]))
riadky2=ifelse(length(which(temp[,4]==bod[2]))==0,NA,which(temp[,4]==bod[2]))
riadok=intersect(riadky1,riadky2)
}
return(riadok)
}
# setting up the variable with ordered points
coords=c(okraj[1,c("x2","y2")])
coords=rbind(coords,c(okraj[1,c("x1","y1")]))
# current working point and a matrix subset, in which to search for a connection
bod=okraj[1,1:2]
temp=okraj[-1,]
# consecutively search for connecting points
for(j in 1:nrow(okraj)){
if(any(is.na(match(bod,temp)))){
message(paste("Problem with: ", bod,", row: ",j ))
next
}
bod2=NA # next ordered point
smer=0 # controls the set of coordinates to use for a connection
for(x in 1:2){
riadok=find.connection(bod=bod,temp=temp)
# my solution to the current crash point can be implemented here
if(is.na(riadok)) bod2=NA
# select correct value for the connecting point from the two sets of coordinates
else bod2=temp[riadok,c(x+smer,x+smer+1)]
if(all(bod==bod2)) bod2=temp[riadok,c(x+2,x+3)]
if(any(is.na(bod2))){
smer=1
next
} else { break }
}
# store the connecting point and move to the next
coords=rbind(coords,bod2)
bod=bod2
# caveats for the last row in the matrix
if(!is.null(dim(temp))) temp = temp[-riadok,]
if(is.null(dim(temp))){
if(all(coords[1,]==coords[nrow(coords),])){
break
}
if(which(temp==bod[1])==1){
coords=rbind(coords,temp[3:4])
break
} else {
coords=rbind(coords,temp[1:2])
break
}
}
}
It works for some datasets. Here, the sample dataset contains an endpoint okraj[8,3:4] that does not connect to anything else in the matrix and the code crashes. My solution would be to skip the row and try again with previous point.
if(is.na(riadok)){
bod2=NA
temp=temp[-riadok,]
bod=coords[nrow(coords),]
next
}
However, the solution requires dropping a row of data, which is wrong. Your help with a systematic solution to converting an alpha shape to a mapped object would be very appreciated.
Edit:
A question about removing holes from polygons received an answer from @Spacedman using rgeos routines. How could it help to resolve my problem with unordered points from the example herein?
