Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Non-overlapping buffers #362

Open
Robinlovelace opened this issue Nov 7, 2019 · 26 comments
Open

Non-overlapping buffers #362

Robinlovelace opened this issue Nov 7, 2019 · 26 comments

Comments

@Robinlovelace
Copy link
Member

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:

image

@mpadge
Copy link
Member

mpadge commented Nov 7, 2019

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
Copy link
Member Author

Seems similar to this: r-spatial/sf#824

@mdsumner
Copy link

mdsumner commented Nov 7, 2019

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

@agila5
Copy link
Collaborator

agila5 commented Nov 8, 2019

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
Copy link
Member Author

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
Copy link
Collaborator

agila5 commented Nov 8, 2019

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
Copy link

wengraf commented Apr 8, 2020

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
Copy link
Member Author

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
Copy link
Member Author

@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
Copy link

wengraf commented Apr 8, 2020

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

@agila5
Copy link
Collaborator

agila5 commented Apr 8, 2020

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
Copy link
Collaborator

agila5 commented Apr 9, 2020

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

@Robinlovelace
Copy link
Member Author

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
Copy link
Collaborator

agila5 commented Apr 11, 2020

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
Copy link
Member

mpadge commented Apr 11, 2020

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

@Robinlovelace
Copy link
Member Author

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
Copy link
Member Author

hypertidy/laridae#8

@mdsumner
Copy link

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
Copy link

oh, no RTriangle is unsuitable nvm

@Robinlovelace
Copy link
Member Author

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
Copy link

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
Copy link

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
Copy link

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
Copy link
Member Author

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
Copy link

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
Copy link
Member Author

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

6 participants