Using r and sf, create an concave envelope POLYGON around MULTILINESTRING

Question

I am interested in creating the minimum possible polygon around a set of nearly parallel cross sections. For example:

mls <- st_multilinestring(list(rbind(c(2,2),c(1,3)),
                    rbind(c(0,0),c(1,1),c(2,1)),
                    rbind(c(3,4),c(0,4))))
plot(mls)

I know there's the function st_convex_hull but that doesn't accomplish what I want because it doesn't generate the minimum possible polygon.

I came up with almost a solution, like this:

mls %>% st_cast(to = "MULTIPOINT") %>% 
st_cast(to = "POLYGON") %>% 
st_simplify() %>%
plot()

Unfortunately, this only produces the desired output for this simple case. But it quickly falls apart with more lines. Here's example bad output:

mls2 <- st_multilinestring(list(rbind(c(2,2),c(1,3)),
                               rbind(c(0,0),c(1,1),c(2,1)),
                               rbind(c(3,4),c(0,4)),
                               rbind(c(1,10),c(-1,12))))
plot(mls2)
mls2 %>% st_cast(to = "MULTIPOINT") %>% 
  st_cast(to = "POLYGON") %>% 
  st_simplify() %>%
  plot()

One idea I had was to pair adjacent lines and generate a polygon for each pair and then merge those polygons, but that feels cumbersome.

Answer 1

I think you can improve your wrapper with this:

1. Use concaveman to get coordinates of the shape
2. Convert it to LINESTRING
3. Snap it to the original lines
4. Convert this to POLYGON

You should adjust the tolerance parameter of st_snap to fit it to your needs:

library(sf)
library(concaveman)
library(dplyr)
library(ggplot2)

# mls-------
mls <- st_multilinestring(list(rbind(c(2,2),c(1,3)),
                               rbind(c(0,0),c(1,1),c(2,1)),
                               rbind(c(3,4),c(0,4))))

conc <- concaveman(st_coordinates(mls)[,1:2]) 
pol  <- st_linestring(conc)
pol <- st_snap(pol, mls, tolerance = 2) %>% st_polygonize()

ggplot() + geom_sf(data=mls,size=2,col="black") +
  geom_sf(data=pol,fill=NA,linetype=2,col="darkgreen",size=1.5) +
  theme_minimal()

# mls2------
mls2 <- st_multilinestring(list(rbind(c(2,2),c(1,3)),
                                rbind(c(0,0),c(1,1),c(2,1)),
                                rbind(c(3,4),c(0,4)),
                                rbind(c(1,10),c(-1,12))))

conc <- concaveman(st_coordinates(mls2)[,1:2]) 
pol  <- st_linestring(conc)
pol <- st_snap(pol, mls, tolerance = 2) %>% st_polygonize()



ggplot() + geom_sf(data=mls2,size=2,col="black") +
  geom_sf(data=pol,fill=NA,linetype=2,col="darkgreen",size=1.5) +
  theme_minimal()

^{Created on 2020-03-19 by the reprex package (v0.3.0)}

Answer 2

The suggestion posed by @TimSalabim in the comments was extremely helpful but even after tweaking the arguments to the concaveman function I couldn't get it to do quite what I wanted. So I combined their suggestion with my first idea in the question (with the help from some code from the github site for the concaveman package. Here's some code demonstrating my solution:

First, concaveman out of the box:

install.packages("concaveman")
library(concaveman)

mls <- st_multilinestring(list(rbind(c(2,2),c(1,3)),
                    rbind(c(0,0),c(1,1),c(2,1)),
                    rbind(c(3,4),c(0,4))))
mls.envelope <- mls %>%
    st_combine() %>%
    st_cast("POINT") %>%
    st_sf() %>%
    concaveman()
plot(mls.envelope)

But this isn't quite what I was looking for..:

So I iterate over the dataset, pairing adjacent lines and creating polygons from those pairs. This happens with every adjacent pair. Then it was necessary to drop the "pairs" that were formed from a single line at the beginning and end of the dataset. Finally st_union to merge all these polygons into a single feature. Cumbersome, but it seems to work.

Here's the code:

concavemanWrapper <- function(multilinestring){
  temp <- multilinestring %>% st_combine() %>%
    st_cast("LINESTRING") %>% st_sf() %>%
    mutate(set1 = ceiling(row_number()/2),
           set2 = floor(row_number()/2)) %>%
    st_cast("POINT")
  
  set1.polys <- map(unique(temp$set1),
                    ~ concaveman(temp[temp$set1 %in% .,])) %>% 
    map2(unique(temp$set1), ~ mutate(.x, set1 = .y)) %>% 
    reduce(rbind)
  
  set2.polys <- map(unique(temp$set2),
                    ~ concaveman(temp[temp$set2 %in% .,])) %>% 
    map2(unique(temp$set2), ~ mutate(.x, set2 = .y)) %>% 
    reduce(rbind)
   
  #Drop the single feature "polygon": 
  set1.polys <- filter(set1.polys, st_is_valid(set1.polys))
  set2.polys <- filter(set2.polys, st_is_valid(set2.polys))
    
  
  fullPolygon <- st_union(set1.polys,set2.polys) %>% st_union()
  return(fullPolygon)
}

mls.envelope2 <- concavemanWrapper(mls)

ggplot() + geom_sf(data=mls,size=2,col="black") +
  geom_sf(data=mls.envelope,fill=NA,linetype=2,col="red",size=1) +
  geom_sf(data=mls.envelope2,fill=NA,linetype=2,col="darkgreen",size=1.5) +
  theme_minimal()

The second MULTILINESTRING from the question (mls2) worked with concaveman out of the box and it also works with this wrapper function: