3

I have a large rivernetwork as a SpatialLines(DataFrame) and want to subset it so I get a SpatialLine with only the shortest/fastest path between two SpatialPoints(DataFrame) which propably need to be snapped onto the line first.

In the example Code, I want to get the closest line between the red dots, while moving on the lines. Also I want to move between the left red point and the right blue point, which needs to be snapped to the closest SpatialLine

library(sp)    
x <- c(1,5,6,8)
y1 <- c(1,3,4,7)
y2 <- c(5,5,5,2)

L <- SpatialLines(list(Lines(Line(cbind(x,y1)), ID="a"),Lines(Line(cbind(x,y2)), ID="b")))
P <- SpatialPoints(data.frame(x=c(1,8),y=c(1,2)))
P_snap <- SpatialPoints(data.frame(x=c(8),y=c(1)))

plot(L)
points(P,col="red")
points(P_snap,col="blue")
1
1

A way to approach this problem could be by switching from sp to sf objects, and then use sfnetworks and tidygraph which internally use igraph.

First we recreate the lines and points as sf objects.

x <- c(1,5,6,8)
y1 <- c(1,3,4,7)
y2 <- c(5,5,5,2)

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
L = st_sf(
  ID = c("a", "b"),
  geometry = c(
    st_sfc(st_linestring(cbind(x,y1))),
    st_sfc(st_linestring(cbind(x,y2)))
  )
)
P = st_as_sf(data.frame(x=c(1,8),y=c(1,2)), coords = c("x","y"))
P_snap = st_as_sf(data.frame(x=8,y=1), coords = c("x","y"))

Then we can create an sfnetwork object from the LINESTRINGs L.

library(sfnetworks)
(G = L %>%
  as_sfnetwork(directed = F))
#> # A sfnetwork with 4 nodes and 2 edges
#> #
#> # CRS:  NA 
#> #
#> # An unrooted forest with 2 trees with spatially explicit edges
#> #
#> # Node Data:     4 x 1 (active)
#> # Geometry type: POINT
#> # Dimension:     XY
#> # Bounding box:  xmin: 1 ymin: 1 xmax: 8 ymax: 7
#>   geometry
#>    <POINT>
#> 1    (1 1)
#> 2    (8 7)
#> 3    (1 5)
#> 4    (8 2)
#> #
#> # Edge Data:     2 x 4
#> # Geometry type: LINESTRING
#> # Dimension:     XY
#> # Bounding box:  xmin: 1 ymin: 1 xmax: 8 ymax: 7
#>    from    to ID                geometry
#>   <int> <int> <chr>         <LINESTRING>
#> 1     1     2 a     (1 1, 5 3, 6 4, 8 7)
#> 2     3     4 b     (1 5, 5 5, 6 5, 8 2)
G %>% class()
#> [1] "sfnetwork" "tbl_graph" "igraph"
plot(G)

This will return a network consisting of the two linestrings as edges, and nodes created at the linestring end points. However, routing will be difficult since there is no node that connects both edges. This will probably not be a problem in the larger river network, but for this particular case we can fix ti with the following code:

G = L %>%
  # Combine LINESTRINGS into a MULTILINESTRING geometry
  st_combine() %>%
  # Create a node where the MULTILINESTRINGs cross
  st_node() %>%
  # Cast back to LINESTRINGs
  st_cast('LINESTRING') %>%
  # Reset IDs
  st_sf(ID = c("a","a","b","b")) %>%
  # Create sfnetwork
  as_sfnetwork(directed = F)

Now, with help of tidygraph we can convert the network to its spatial shortest paths. Convert will result in a subset of the original network, consisting only of the spatial shortest path in this case. sfnetworks accepts sf POINTs as to and from arguments, while internally snapping this points to the nearest node. It also uses the edge geographical length as weight internally

library(tidygraph)
#> 
#> Attaching package: 'tidygraph'
#> The following object is masked from 'package:stats':
#> 
#>     filter
# Convert to spatial shortest path between points in the network
G_sp1 = G %>%
  convert(to_spatial_shortest_paths, from = P[1,], to = P[2,])
# Convert to spatial shortest path to snapped point
G_sp2 = G %>%
  convert(to_spatial_shortest_paths, from = P[1,], to = P_snap)

We can plot the result to see how it worked

par(mar = c(0,0,0,0), mfrow = c(1,2))

plot(G, col = "grey", reset = FALSE)
plot(P, pch = 8, col = "green", add = TRUE)
G_sp1 %>%
  plot(col = "blue", add = TRUE)

plot(G, col = "grey", reset = FALSE)
plot(P, pch = 8, col = "green", add = TRUE)
plot(P_snap, pch = 8, col = "red", add = TRUE)
G_sp2 %>%
  plot(col = "purple", add = TRUE)

If you still need to obtain the shortest path as a SpatialLine, it can be extracted as follows:

(L_sp1 = G_sp1 %>%
  st_geometry(active = "edges") %>%
  as_Spatial())
#> An object of class "SpatialLines"
#> Slot "lines":
#> [[1]]
#> An object of class "Lines"
#> Slot "Lines":
#> [[1]]
#> An object of class "Line"
#> Slot "coords":
#>          [,1] [,2]
#> [1,] 1.000000  1.0
#> [2,] 5.000000  3.0
#> [3,] 6.000000  4.0
#> [4,] 6.333333  4.5
#> 
#> 
#> 
#> Slot "ID":
#> [1] "ID1"
#> 
#> 
#> [[2]]
#> An object of class "Lines"
#> Slot "Lines":
#> [[1]]
#> An object of class "Line"
#> Slot "coords":
#>          [,1] [,2]
#> [1,] 6.333333  4.5
#> [2,] 8.000000  2.0
#> 
#> 
#> 
#> Slot "ID":
#> [1] "ID2"
#> 
#> 
#> 
#> Slot "bbox":
#>   min max
#> x   1 8.0
#> y   1 4.5
#> 
#> Slot "proj4string":
#> CRS arguments: NA
3

apologies for all these libraries, i'm not sure exactly which is needed... Good luck! this is kind of complicated. and not exactly sure how it works.

library(igraph)
library(RColorBrewer)
library(SpatialTools)
library(sp)
library(ANN)
library(rgeos)
library(maptools)
library(sp)
library(rgeos)
library(plyr)
library(FNN)
library(tidyr)
library(foreign)
library(dplyr)
library(plyr)
library(reshape)


buildTopo <- function(lines) {

  g = gIntersection(lines, lines)
  edges = do.call(rbind, lapply(g@lines[[1]]@Lines, function(ls) {
    as.vector(t(ls@coords))
  }))
  lengths = sqrt((edges[, 1] - edges[, 3])^2 + (edges[, 2] - edges[, 4])^2)

  froms = paste(edges[, 1], edges[, 2])
  tos = paste(edges[, 3], edges[, 4])

  graph = graph.edgelist(cbind(froms, tos), directed = FALSE)
  E(graph)$weight = lengths

  xy = do.call(rbind, strsplit(V(graph)$name, " "))

  V(graph)$x = as.numeric(xy[, 1])
  V(graph)$y = as.numeric(xy[, 2])
  return(graph)
}

routePoints <- function(graph, from, to) {
  require(FNN)
  xyg = cbind(V(graph)$x, V(graph)$y)

  ifrom = get.knnx(xyg, from, 1)$nn.index[1, 1]
  ito = get.knnx(xyg, to, 1)$nn.index[1, 1]

  p = get.shortest.paths(graph, ifrom, ito, output = "vpath")
  p[[1]]
}

buildTopo() use this function to create a topographically correct object. The input is a SpatialLinesDataFrame

Then use routePoints(ResultOfbuildTopo, FromCoord, ToCoord)

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.