1

I am trying to create a dataset of roads data to do shortest path analysis. To do this, my roads data needs to be connected: it should be possible to jump from one road to another. I am planning to use the R library shp2graph for this. This library contains a function nt.connect that finds the largest connected network. However, this function finds that my data hardly overlaps, despite appearing to overlap on visual inspection. The function returns the following image to visualise these separate networks: unconnected roads network

How can I connect these overlapping pieces of road to each other?

What I've tried:

  • Buffering the roads network. This returns a SpatialPolygon object which nt.connect does not accept: it wants a SpatialLines object. If I would convert the SpatialPolygon object to a SpatialLines object it would take the boundary of the Polygon, which would result in two roads right next to each other. I am worried this will mess with the subsequent analysis.

Some notes:

  • I am using OpenStreetMap data. The dataset (a subset of my full dataset) can be downloaded here: https://www.dropbox.com/sh/qu2fg3f9fy5tk7k/AABjv81Vv1R9c0RHuMJwIHH2a?dl=0
  • This is only a (very small) portion of the full roads network I want to use. Manually making sure they overlap is therefore not feasible
  • This library I am planning to use uses R. My work needs to be fully reproducible and this is difficult when working with a GUI.
  • Some pieces of the road are true "orphans": they do not overlap with the larger network. I would prefer if these can be incorporated into the main network, for example by drawing the shortest possible line to the full network. I would also be fine with just dropping these orphans.

Below is the code to check that the network is connected:

require(sf) 
require(sp)
require(shp2graph)

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
roads <- read_sf("../GIS Stack example", "roads")
roads_sp <- as(roads, "Spatial")
roads_connected <- nt.connect(roads_sp)
  • If you want to also ask about QGIS and Python then please do so in separate focused questions. – PolyGeo Sep 4 at 10:53
0

You can transform your road network into a raster whose values are 1 for pixels over the network and NA otherwise. Then, you can use the package gdistance to transform the raster into a transition matrix and then calculate the shortest path over it. Try

library(sf) 
library(sp)
library(raster)
library(gdistance)

# loading your road data:
roads <- read_sf("../GIS Stack example", "roads")
roads_sp <- as(roads, "Spatial")

# creating empty raster:
r <- raster()
res(r) <- 1 # decrease this to make the calculations more precise (and more computer demanding)

# re-projecting raods to raster's projection:
roads_sp <- spTransform(roads_sp, r@crs)

# rasterizing the rodas' dataset:
r <- setExtent(r, roads_sp@bbox)
r <- rasterize(roads_sp, r)
r[!is.na(r)] <- 1 # setting cost 0 to travel over roads_sp
plot(r) # plot of the road network as a raster

# transforming it into a transition matrix to calculate distances:
r <- transition(r, mean, directions = 8)
r <- geoCorrection(r, "c")

# to calculate the distance between two points:
pp <- SpatialPoints(coords = matrix(c(-11.0874,7.9956, -11.10635, 7.93558),2,2,byrow=T), proj4string = roads_sp@proj4string) # random two points
# plot of the two points for illustration:
plot(pp, add=T)

# calculating distances and plotting the shortest path:
distance <- shortestPath(r, pp[1,], pp[2,], output = "SpatialLines")
plot(distance, add=T)

The output is the image below. The step rasterize() can be computationally demanding and you might want to check the fasterize package as an alternative.enter image description here

  • Thanks! I'm running into the problem that sometimes points aren't exactly on the line, and then shortestPath complains that some vertices can't be reached: In get.shortest.paths(adjacencyGraph, indexOrigin, indexGoal) : At structural_properties.c:4597 :Couldn't reach some vertices – Paul Hofman Sep 17 at 12:07
  • Also, is there any way to also calculate distance along the route? I've been trying this by converting to a UTM CRS and then using rgeos::gLength, but I've been stopped by the problem in the previous comment – Paul Hofman Sep 17 at 12:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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