Non-overlapping buffers

[Robinlovelace]

This may not be the best place to solve this, but I can imagine non overlapping buffers, that represent the area in which each line/feature is the nearest one with no overlaps, would be useful. Reproducible example plus sketch of what I'm thinking below.

# question about buffer types
library(sf)
#> Linking to GEOS 3.7.1, GDAL 2.4.2, PROJ 5.2.0
library(stplanr)
l1 = stplanr::osm_net_example[1, ]
l = stplanr::osm_net_example[l1, ]
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
lb = geo_projected(shp = l, fun = st_buffer, dist = 10)
lb_flat = geo_projected(shp = l, fun = st_buffer, dist = 10, endCapStyle = "FLAT")
plot(st_geometry(l))
plot(st_geometry(lb), col = sf.colors(nrow(l), alpha = 0.5), add = TRUE)

plot(st_geometry(l))
plot(st_geometry(lb_flat), col = sf.colors(nrow(l), alpha = 0.5), add = TRUE)

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

Very rough sketch of solution:

[mpadge]

Have a look at the great code by @mdsumner in mapscanner that does the exact opposite of that - aggregates overlapping polygons, and returns contours for degrees of overlap. That must have exactly what you need buried in there somewhere.

[Robinlovelace]

Seems similar to this: https://github.com/r-spatial/sf/issues/824

[mdsumner]

ms_aggregate_polys() can't do that because those edges (the ones that split the overlap) don't exist in the layer

[agila5]

Hope it makes sense. I'm quite sure it's possible to recode it without dplyr/map2 but that's not the main problem now... I added some comments in the code to explain it.

# packages
library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
library(stplanr)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(purrr)

# data
osm_net_example <- st_transform(osm_net_example, 27700) 
l1 <- osm_net_example[1, ]
example <- osm_net_example[l1, ] %>% mutate(my_ID = as.character(seq_len(nrow(.)))) %>% 
  st_set_agr("constant")

# change example data to points
example_points <- st_cast(example, "POINT")
example_unique_points <- example_points %>% 
  filter(!duplicated(geometry))
# I don't need the points shared between two or more lines and, most important,
# I need a unique ID for each point later
example_unique_multipoints <- example_unique_points %>% group_by(my_ID) %>% summarise()

voronoi_polygons <- st_voronoi(
  # https://github.com/r-spatial/sf/issues/824
  x = do.call("c", st_geometry(example_unique_multipoints))
  ) %>% 
  st_collection_extract() %>% 
  st_set_crs(27700)

voronoi_polygons_for_lines <- voronoi_polygons %>% 
  st_set_crs(27700) %>% 
  # now I need to merge the polygons associated with the same line
  # https://github.com/r-spatial/sf/issues/1030
  st_cast("MULTIPOLYGON", ids = unlist(st_intersects(voronoi_polygons, example_unique_multipoints))) %>% 
  st_union(by_feature = TRUE)

par(mar = rep(0, 4))
plot(voronoi_polygons_for_lines, col = sf.colors(3), reset = FALSE)
plot(st_geometry(example), lwd = 4, add = TRUE, col = "black")

# Now calculate the buffers
example_buffer <- st_buffer(example, dist = 20)

# Now I'd like to take the intersection of each buffer with the corresponding
# voronoi polygon. The problem is that if I run
st_intersection(st_geometry(example_buffer), voronoi_polygons_for_lines)
#> Geometry set for 9 features 
#> geometry type:  POLYGON
#> dimension:      XY
#> bbox:           xmin: 430879.8 ymin: 434314.2 xmax: 430989.7 ymax: 434456.9
#> epsg (SRID):    27700
#> proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
#> First 5 geometries:
#> POLYGON ((430921 434400.7, 430938.2 434431.7, 4...
#> POLYGON ((430909.3 434382.8, 430915.3 434388.1,...
#> POLYGON ((430910.4 434381.6, 430910.7 434382.1,...
#> POLYGON ((430916.7 434357, 430916.2 434357.3, 4...
#> POLYGON ((430905.4 434379.3, 430909.3 434382.8,...
# I get the intersection of each buffer with each polygon

# I'm not sure how to solve that. For the moment I use map2
my_solution <- map2(
  .x = st_geometry(st_buffer(example, dist = 20)), 
  .y = voronoi_polygons_for_lines, 
  .f = st_intersection
  ) %>% 
  st_sfc()

# this is the result
plot(my_solution, col = sf.colors(length(my_solution), alpha = 0.5))

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

[Robinlovelace]

Nice solution! That works for sure, just thinking about how to deal with the artefact in the centre. 2 options:

1. create a circle with diameter equal to the buffer distance
2. or create artificial points that are very close the intersections associated with each line

The first is simpler and links to the idea of treating junctions differently than open roads.

[agila5]

One update on the first option:

# The code before here is the same as with the other comment. 

# I'm not sure how to solve that. For the moment I use map2
my_solution <- map2(
  .x = st_geometry(st_buffer(example, dist = 20)), 
  .y = voronoi_polygons_for_lines, 
  .f = st_intersection
) %>% 
  st_sfc(crs = 27700) # I forgot about this in the previous answer

# look for duplicated points.
duplicated_points <- st_geometry(example_points)[duplicated(st_geometry(example_points))][1]
duplicated_points_buffer <- st_buffer(duplicated_points, 20)

plot(my_solution, col = sf.colors(length(my_solution), alpha = 0.5), reset = FALSE)
plot(st_buffer(duplicated_points, 20), add = TRUE)

# take difference between each buffer and the buffers at the junctions
my_solution2 <- st_difference(my_solution, duplicated_points_buffer) %>% 
  c(duplicated_points_buffer)
plot(my_solution2, col = sf.colors(length(my_solution2), alpha = 0.5))

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

Actually I shouldn't look for duplicated points but for "junction points" (i.e. points that are repeated 3 or more times as @mpadge said in #357 ) but first I should check if 1) it's useful and it does make sense (for example, st_difference returns MULTIPOLYGONS instead of POLYGONS and I'm not sure if that's important) and 2) if it works in more difficult cases. Next week...

[wengraf]

I just want to document here that the above has been hugely helpful to me. I wanted to join a large point dataset (2.5m observations, all on-road but not geocoded to the road centre line itself) to OS OpenRoads, to make a point dataset with attributes from OpenRoads. My previous method using st_nn just took way too long. I have not done everything @agila5 has set out, because with a complete national road network and points that must have been on-road, I don't need the buffers - just a nearest neighbour match made by joining points to voronois. @agila5 : are you planning on adding this as a function?

[Robinlovelace]

I think we can move on that basis. If you're happy to test the resulting rnet_buffer() (or whatever it ends up being called) function let us know @agila5. Otherwise I'm happy to give it a first go.

[Robinlovelace]

@wengraf you've got an eye for names. rnet_buffer() sound good to you. See here for other rnet() function names: https://docs.ropensci.org/stplanr/reference/index.html

[wengraf]

That sounds good to me. Happy to test it on my point data.

[agila5]

Hi! As i said on twitter, thank you very much for testing those ideas. If I don't get extremely bad results with the other project I plan to submit a PR on saturday or sunday (my way for celebrating Easter 😅).

[agila5]

On a second thought I think that probably it's better to wait for sfnetworks integration. @Robinlovelace thoughts?

[Robinlovelace]

On a second thought I think that probably it's better to wait for sfnetworks integration. @Robinlovelace thoughts?

I think this is a purely geographic operation that doesn't rely, at present, on sfnetworks. Spatial network analysis could help identify the junction points, but we can identify them without sfnetworks already I think.

[agila5]

Actually I think that the problem is much more complicated than what I expected 😅 I created a new question on GIS SO that summarises the most significant problem. Questions and suggestions are welcome.

[mpadge]

Pinging @mdsumner for an answer to that So question...

[Robinlovelace]

This issue is raising lots of interesting questions about what junctions are, how to prevent overlaps in buffers and much more... Looking forward to seeing what comes of this. I think there are simple ways of tackling this in the 'no new features' scenario using st_nearest_feature() but I think the 'joining to junctions and segments away from junctions' solution is more interesting and useful.

[Robinlovelace]

https://github.com/hypertidy/laridae/issues/8

[mdsumner]

hey that's pretty cool in SO ... if I get to this I'll be looking at RTriangle - it returns the V(oronoi) as well, something I tend not to explore but should.

(Something I glanced close to recently was a way to easily call RTriangle with any kind of input format - maybe I'l dig that out first)

[mdsumner]

oh, no RTriangle is unsuitable nvm

[Robinlovelace]

Do you mean RTriangle cannot do Segment Voronoi Diagrams: https://doc.cgal.org/Manual/3.1/doc_html/cgal_manual/Segment_Voronoi_diagram_2/Chapter_main.html ?

That would be ideal, but wonder if there is a way of hacking RTriangle to get it to produce the output we need...

[mdsumner]

I don't think so.

You could use raster::distance to have a go at whuber's image version, but you'll end up with nasty jaggy edges (now I'm thinking about isoband)

[mdsumner]

whuber's raster example inspired me so I had a go, it's not brilliant but might be of use

  library(sf)
#> Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3
  library(raster)
#> Warning: package 'raster' was built under R version 3.6.3
#> Loading required package: sp
#> Warning: package 'sp' was built under R version 3.6.3
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))
)

isoband_raster <- function(x, lo, hi, auto = FALSE) {
  if (auto) {
    breaks <- pretty(values(x))
    lo <- head(breaks, -1)
    hi <- tail(breaks, -1)
  }
  ## OMG: note the [[1]] to avoid the case of as.matrix(brick) which is not helpful ...
  b <- isoband::isobands(xFromCol(x), yFromRow(x), as.matrix(x[[1]]), levels_low = lo, levels_hi = hi)
  st_cast(sf::st_sf(lo = lo, hi = hi, geometry  = sf::st_sfc(isoband::iso_to_sfg(b), crs = projection(x))), 
          "POLYGON")
}

sf <- st_sf(g = my_linestring_sfc)
nxy <- c(256, 256)
r <- raster(sf, ncols = nxy[1], nrows = nxy[2])
d <- distance(rasterize(sf, r))
#> Warning in .couldBeLonLat(x, warnings = warnings): CRS is NA. Assuming it is
#> longitude/latitude
d[] <- scales::rescale(d[])
d[d == 0] <- -.1
plot(isoband_raster(d, 0, 1)$geometry, col = hcl.colors(3))
#> Warning in st_cast.sf(sf::st_sf(lo = lo, hi = hi, geometry =
#> sf::st_sfc(isoband::iso_to_sfg(b), : repeating attributes for all sub-geometries
#> for which they may not be constant

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

[fBedecarrats]

Thanks for these great reflections. I found an interesting solution here on QGIS, associated with a demo video, in the line of what @agila5 had proposed. I guess this demonstrates a procedure that could be implemented in R.

[Robinlovelace]

Hi @fBedecarrats thanks for sharing, this is indeed a nice solution. Variants with and without 'junction buffers' could be useful. From a road safety perspective many crashes happen at junctions which are not on one road segment or another so we'd want circular points at each intersection. Heads-up @agila5 I want to implement your solution above but first some refactoring of {stplanr} is in order, starting with removing dependencies on {rgeos} and {rgdal}...

[fBedecarrats]

Thanks for your feedback @Robinlovelace. My use case is somewhat different (affect gps sensed routes crowdsourced by the users of a mobile app to cycling facilities), but I share this need to treat specifically the intersections. I found a solution for non-overlapping buffers already implemented by @statnmap in {cartomisc} with a function named regional_seas. Its inteded purpose was to project buffers from polygons, but it works perfectly with lines. I tried editing the call to st_buffers() with an argument endCapStyle = "FLAT", but some results are weird.

[Robinlovelace]

Great thanks for the info. Heads-up: I may ping you here when I have things to test. For now my priority is #481 / #332