16

I have in my possession four sets of shapefiles:
1. Locations (points), around 1000
2. Roads (lines), each road has a speed attribute
3. Terrain (polygons), each zone also has a speed attribute
4. Rivers (lines), each river has a crossing delay time attribute

How can I construct a Voronoi diagram for my location using travel time as a distance metric? I could have used pgRouting if I had roads only, but I am not sure how to approach different off-road terrain and rivers. I'm willing to code this, but I've spent enough time trying to come up with an algorithm, but to no avail.

My current plan is to dump it all into a sufficiently large raster and create 1000 distance maps out of it, then use them to calculate the Voronoi cells, but I have a suspicion that a sufficiently large raster will take too long to process. Is there a better vector-based algorithm?

Improve this question
6
  • 2
    I describe an extremely efficient solution in a recently awarded patent viewable at google.com/patents/….
    – whuber
    Jun 24, 2013 at 16:43
  • @whuber Oh, interesting. It seems you're also rasterizing the region. I am not sure if I can use a patented algorithm, though.
    – Alexey
    Jun 25, 2013 at 9:55
  • 1
    @dassouki That's not how patents work :-(. You have to demonstrate priority and exclusiveness, which means you can't publish beforehand. I do have presentation slides, though: look for the "raster-vector hybrids" section at the end (part 4, pp 37-45 inclusive). To fully appreciate slide 38 ("Why not put the two together?"), though, you should look at the first three parts :-).
    – whuber
    Jun 25, 2013 at 13:35
  • 2
    It appears that Whuber answered this question as a comment. @whuber please consider adding your comment as an answer.
    – kttii
    Oct 17, 2016 at 16:16
  • 1
    @kttii Thank you for reminding me of this thread. I only gave a reference: that doesn't qualify as an answer. Someday if I get the time I'll see about providing an actual answer.
    – whuber
    Oct 17, 2016 at 17:45

2 Answers 2

Reset to default
2

This question gave me some ideas to write a post about it, trying to get the R tools to perform such a thing. Here you can check all details. It is in spanish but you can use the translation of the browser.

This is the resume of what I tryed:

  1. Convert all layers into points by sampling the lines layer (streets & rivers) and by sampling the polygons (areas with speed restrictions) by a regular grid. At this points you can either change the speed of streen by the polygon one or keep both points with different speeds.
  2. Merge the speedpoints layers.
  3. Create an empty raster grid and rasterize the points using it.
  4. Calculate the accumulated cost with the gdistance library.
  5. (optional but as you asked for Voronoi...) Convert back the raster into points and calculate the polygons with st_voronoi.
  6. Consolidate this process into a functiton that can by used in parallel over all you points.

Here, speeds is the raster of speed values:

# Read raw data
locs <- st_read("./data/20230126/data-raw/LEIP/LEIP.shp")
roads <- st_read("./data/20230126/data-raw/Carreteras_IMD/Carreteras_IMD.shp")
terrain <- st_read("./data/20230126/data-raw/AEIP/AEIP.shp")
rios <- st_read("./data/20230126/data-raw/Rios/Rios.shp")

# Create random speeds for features
r <- r %>% mutate(SPEED = sample(50, size = nrow(r))) %>% st_cast("LINESTRING")
w <- w %>% mutate(SPEED = sample(500, size = nrow(w))/10) %>% st_cast("LINESTRING")
t <- t %>% mutate(SPEED = sample(50, size = nrow(t)))


wpoints <- w %>% st_cast("POINT")
rpoints <- r %>% st_cast("POINT")
tpoints <- st_sample(t, size = 500, type = "regular") %>% st_as_sf() %>% st_join(t, join = st_nearest_feature)

# merge all points with speeds
speedpoints <- bind_rows(wpoints, rpoints, tpoints) %>% dplyr::select(SPEED)

#función para empty raster (raster base)
library(raster)
library(fasterize)
rasterbase <- function(layer, res){
    r <- raster(extent(layer), res)
    res(r) <- c(res, res)
    crs(r) <- crs(layer)
    values(r) <- NA
    return(r)
}
raster <- rasterbase(l, 50)

g = st_combine(st_geometry(speedpoints))
x <- st_voronoi(g, bOnlyEdges = FALSE)
v = st_as_sf(st_collection_extract(x)) %>%
    st_join(speedpoints, join = st_nearest_feature) %>% 
    st_make_valid() %>% 
    st_cast("POLYGON") %>% 
    st_crop(l)

speeds <- fasterize(v, raster, field = "SPEED")
        
tr <- transition(speeds, transitionFunction=mean, directions=8)
    trc <- geoCorrection(tr, type="c", multpl=FALSE, scl=TRUE)

Here is the function for travel time for one starting point:

traveltime <- function(p, trc){
    # p <- p <- l[1,]
    ac <- accCost(trc, st_coordinates(p))
    ttras <- round(ac)
    ttpol <- polygonizator(ttras) %>% rename(TIME = layer)
    ttpol$ID <- p$ID
    return(ttpol)
}

And here, lsel is the layer with Locations (points), around 1000:

lsel <- l
library(parallel)
cl <- parallel::makeCluster(detectCores(), type="FORK")
doParallel::registerDoParallel(cl)
TTIMES <- parLapply(cl, 1:nrow(lsel), function(x){traveltime(lsel[x,], trc)})
stopCluster(cl)

enter image description here

enter image description here

enter image description here

May not be a perfect solution, but as an aproximation can be useful. May be others can develop this further.

Hope it helps.

Improve this answer
0

These days, file up ArcGIS Online or ArcGIS Pro, and use the travel-time buffer tool.

Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.