R - Voronoi polygons from lines

Question

I use R (and sf, mainly) and I was wondering if it's possible to create Voronoi polygons starting from LINESTRING objects instead of POINT. For example, given a sfc such as

# package
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

# example linestring
my_linestring_sfc <- st_sfc(
  st_linestring(matrix(c(-2, -2, -1, -1, 0, 0), ncol = 2, byrow = TRUE)), 
  st_linestring(matrix(c(2, -2, 1, -1, 0, 0), ncol = 2, byrow = TRUE)), 
  st_linestring(matrix(c(0, 0, 0, 1, 0, 2), ncol = 2, byrow = TRUE))
)

par(mar = rep(0, 4))
plot(my_linestring_sfc, reset = FALSE)
plot(st_cast(my_linestring_sfc, "POINT"), add = TRUE, pch = 16)

I can create the Voronoi polygons associated with its points using the following code:

my_voronoi <- st_voronoi(do.call("c", st_cast(my_linestring_sfc, "POINT"))) %>% 
  st_collection_extract()

plot(my_voronoi, reset = FALSE)
plot(my_linestring_sfc, add = TRUE)
plot(st_cast(my_linestring_sfc, "POINT"), add = TRUE, pch = 16)

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

and I can merge the polygons according to their intersection with the original linestrings. The problem is that this tessellation does not correspond to my "intended" result since the middle polygon is shared between the three lines. The ideal result should be as follows (which I created manually):

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

I read a few SO posts on this topic. This is a solution with Mathematica software while this is another solution with ArcGis, and they both use the Euclidean Allocation. In this other question they suggest "Segment Voronoi Diagrams" and CGAL. This looks like the most promising solution using PostGIS, but I have no idea how to translate that code into R. Is it possible to adapt those ideas to my problem?

Note: Check here since maybe it's possible to adapt etiennebr approach to the previous PostGIS answer.

Answer 1

Steps towards a solution, but the centre is not right:

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0

# example linestring
my_linestring_sfc <- st_sfc(
  st_linestring(matrix(c(-2, -2, -1, -1, 0, 0), ncol = 2, byrow = TRUE)), 
  st_linestring(matrix(c(2, -2, 1, -1, 0, 0), ncol = 2, byrow = TRUE)), 
  st_linestring(matrix(c(0, 0, 0, 1, 0, 2), ncol = 2, byrow = TRUE))
)
my_linestring_sf <- st_sf(data.frame(id = 1:3), geometry = my_linestring_sfc)
my_linestring_sf_buffer <- st_buffer(my_linestring_sf, 0.2) %>% 
  st_union()

par(mar = rep(0, 4))
plot(my_linestring_sfc, reset = FALSE)
p <- st_cast(my_linestring_sf, "POINT")
#> Warning in st_cast.sf(my_linestring_sf, "POINT"): repeating attributes for all
#> sub-geometries for which they may not be constant
plot(p, add = TRUE, pch = 16)
plot(my_linestring_sf_buffer, add = TRUE)

junctions = stplanr::rnet_get_nodes(my_linestring_sf)
sel_in_junctions <- paste(st_coordinates(p$geometry)[, 1], st_coordinates(p$geometry)[, 2]) %in% 
  paste(st_coordinates(junctions)[, 1], st_coordinates(junctions)[, 2])

p_not_junctions = p[!sel_in_junctions, ]

plot(p_all)
#> Error in plot(p_all): object 'p_all' not found
plot(junctions, add = TRUE, cex = 4)
plot(p_not_junctions, add = TRUE, cex = 2)
p_multi_not_junctions <- do.call("c", p_not_junctions$geometry)
plot(my_linestring_sf, add = TRUE)

v <- st_voronoi(do.call("c", p_not_junctions$geometry))%>% 
  st_collection_extract()
v_joined = st_join(st_sf(v), p_not_junctions)
plot(v_joined)

v_intersection = st_intersection(v_joined, my_linestring_sf_buffer)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
plot(v_intersection)

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

Answer 2

A rather late answer, but if you rasterize your lines on a grid you can then use EBImage to construct these generalized voronoi diagrams very quickly, and then polygonizing the raster.

EBImage's voronoization is explained here enter link description here under section 9.2. It uses a modified version of the Euclidean distance that includes a weighting term, you'll want to explode the weighting term (e.g. by setting lambda = 1e15).

Answer 3

Here's another solution. I don't think this is the best option because it relies on old sp code (that I need to rewrite in sf for stplanr) and is probably very slow. But is this getting towards the result you were hoping for?:

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0

    # example linestring
    my_linestring_sfc <- st_sfc(
        st_linestring(matrix(c(-2, -2, -1, -1, 0, 0), ncol = 2, byrow = TRUE)), 
        st_linestring(matrix(c(2, -2, 1, -1, 0, 0), ncol = 2, byrow = TRUE)), 
        st_linestring(matrix(c(0, 0, 0, 1, 0, 2), ncol = 2, byrow = TRUE))
    )
    my_linestring_sf <- st_sf(data.frame(id = 1:3), geometry = my_linestring_sfc)
    buffer_size <- 0.2
    my_linestring_sf_buffer <- st_buffer(my_linestring_sf, 0.2) %>% 
        st_union()

    junctions <- stplanr::rnet_get_nodes(my_linestring_sf)
    junction_buffers <- st_buffer(junctions, buffer_size / 10)
    junction_buffers <- st_buffer(junctions, buffer_size)
    my_linestring_sf_toptailed = stplanr::toptail_buff(
        as(my_linestring_sf, "Spatial"),
        buff = as(junction_buffers, "Spatial")
    ) %>% 
        st_as_sf()
    my_linestring_sf_toptailed$id = 1:3
    plot(my_linestring_sf_toptailed$geometry)
    plot(junction_buffers, add = TRUE)

    p_all = st_cast(my_linestring_sf_toptailed$geometry, "POINT") %>% 
        do.call(what = "c") %>% 
        st_union()

    v <- st_voronoi(p_all) %>% 
        st_collection_extract() %>% 
        st_sf()
    plot(v)
    plot(my_linestring_sf_toptailed, add = TRUE)

    v_joined <- st_join(v, my_linestring_sf_toptailed)
    v_aggregated <- aggregate(v_joined, by = list(v_joined$id), FUN = mean)
    v_intersection <- st_intersection(v_aggregated, my_linestring_sf_buffer)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries
    plot(v_intersection)

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