st_network_cost - How to index "to" argument when "duplicated values in argument 'to' were removed. "

[zackarno]

Hello again!

When using the st_network_cost function sometimes I want to visualize some particular shortest paths returned from from the cost matrix. This is straightforward when nrow(cost_matrix)= nrow(from_argument) & ncol(cost_matrix)=nrow(to_argument). However, often, the ncol(cost_matrix) < nrow(to_argument). This seems to make it difficult to index the closest value in the to_argument for a given shortest path for subsequent visualization

I believe the issue arises due to this warning: "duplicated values in argument 'to' were removed" and the code below
if (any(duplicated(to))) {
warning("Duplicated values in argument 'to' were removed.",
call. = FALSE)
to = unique(to)
}

It makes sense to me to remove the duplicated "to" nodes prior to the calculation, but I am wondering what is the best way to use the returned cost_matrix to to index the closest "to" argument for a particular shortest path.

I hope this makes sense, thanks again for the great work on this!

[luukvdmeer]

Hi and thanks for your question! I looked again at the documentation of st_network_cost() and noticed this behaviour is indeed not documented good enough. I will add some lines about this.

The main problem is the following: st_network_cost() is a wrapper around igraph::distances(). The difference is that the latter only accepts indices/names of existing nodes in the from and to argument, while in st_network_cost() you can provide spatial point geometries at any location. Internally, it will find then the nearest node to each of these points, and forward the indices of these nearest nodes to igraph::distances(). However, igraph::distances() simply does not allow duplicated node indices in its to argument. That means we are forced to remove them before calling igraph::distances(). Hence, removing duplicated "to" points is not really a design choice we made when writing sfnetworks, but rather a requirement from igraph (interestingly they don't put this limitation on the "from" argument). This can be a bit annoying when you are an sfnetworks user, since you don't know beforehand if some of your provided point geometries will have the same node as nearest node.

The only way around this currently, is that you do the step of finding the indices of the nearest nodes yourself. In that way you have the overview of which of your provided "to" points are actually missing in the cost matrix.

library(sfnetworks)
#> Registered S3 method overwritten by 'spatstat.geom':
#>   method     from
#>   print.boxx cli
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

# Create a network.
net = as_sfnetwork(roxel, directed = FALSE, length_as_weight = TRUE)

# Create a set of points to calculate the cost matrix with.
p1 = st_sfc(st_point(c(7.5321, 51.9533)), crs = 4326)
p2 = st_sfc(st_point(c(7.5322, 51.9532)), crs = 4326)
p3 = st_sfc(st_point(c(7.5337, 51.9555)), crs = 4326)
pts = c(p1, p2, p3)

# Find the nearest feature in the nodes table of the network to each point.
# Note that the first two points have the same nearest node.
nearest = st_nearest_feature(pts, st_geometry(net, "nodes"))
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
nearest
#> [1] 4 4 1

# Retrieve the cost matrix.
# Note that we retrieve a 3x2 matrix instead of a 3x3 one.
# That is: igraph forced us to remove duplicated "to" indices, but *not* duplicated "from" indices.
costs = st_network_cost(net, from = nearest, to = nearest)
#> Warning: Duplicated values in argument 'to' were removed.
costs
#>          [,1]     [,2]
#> [1,]   0.0000 354.6676
#> [2,]   0.0000 354.6676
#> [3,] 354.6676   0.0000

# Map each of your original point geometries to a column index in the cost matrix.
# The second point will have the cost values towards it stored in the first column, since it has the same nearest node as the first point.
col_idxs = match(nearest, unique(nearest))
col_idxs
#> [1] 1 1 2

# So if we want to obtain the travel cost from the third point to the second point, we find is as follows:
costs[3, col_idxs[2]]
#> [1] 354.6676

^{Created on 2021-05-13 by the reprex package (v2.0.0)}

Not sure if this exactly answers your question, but it may help you a step further.

[zackarno]

Hey,

Thanks for the detailed answer and helpful work around! I now understand how the limitation on the to argument in igraph is creating this issue. After reading through your workaround I developed a slightly modified work-flow so that I would be able to look at the cost distance matrix, check a few random points using st_network_paths (to make sure both functions returned the same shortest distance), and then visualize the shortest route. I copied the example below. While the results seem good with both functions returning the same shortest distance, one thing I found a bit strange was that it really does not look like the shortest route when plotted (in code below). Perhaps its a node snapping issue that would be fixed by st_network_blend, but I'm not positive.

library(sfnetworks)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(tidygraph)
#> 
#> Attaching package: 'tidygraph'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(tidyverse)


set.seed(799)
facilities<-st_sample(x = st_bbox(roxel) %>%
                        st_as_sfc(),
                      size = 100) %>% 
  st_as_sf()
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar

sites<-st_sample(x = st_bbox(roxel) %>%
                   st_as_sfc(),size = 10) %>% 
  st_as_sf()
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar


plot(sites, col="red")
plot(facilities,col="blue",add=T)
plot(roxel, col="grey",add=T)
#> Warning in plot.sf(roxel, col = "grey", add = T): ignoring all but the first
#> attribute

net<-as_sfnetwork(roxel,directed=FALSE,length_as_weight = TRUE)

net_main <- net %>% 
  filter(group_components()==1)

#add col with closest node
sites$closest_net_node<- sfnetworks:::set_path_endpoints(x = net_main,p = sites)
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar
facilities$closest_net_node<-sfnetworks:::set_path_endpoints(x = net_main,p = facilities)
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar

# cost distance matrix
cost_many<-st_network_cost(x = net_main,from = sites, to = facilities)
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar
#> Warning: Duplicated values in argument 'to' were removed.
cost_many %>% dim()
#> [1] 10 80

# create unique cost index
cost_idx<- unique(facilities$closest_net_node)
length(cost_idx)
#> [1] 80

# for this example lets look at site # 5 & its closest facility
which.min(cost_many[5,])
#> [1] 24
cost_many[5,][24]
#> [1] 274.6754

# node of interest
noi<-cost_idx[which.min(cost_many[5,])]

# now i want to compare it to st_network_paths to make sure i get the same distance
# and so that I am able to visualize the shortest route
path  <- sfnetworks::st_network_paths(x = net_main,
                                      from =sites %>% slice(5),
                                      to=facilities %>% 
                                        filter(closest_net_node==noi)
)
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_points assumes that they are planar


nps <- path %>% 
  pull(node_paths)# %>% unlist()

nps<- nps[[1]] %>% 
  unlist()

# calculate the total edge lengths of shortest calculated route - its the same!
net_main %>% 
  activate("nodes") %>% 
  slice(nps) %>% 
  activate("edges") %>% 
  as_tibble() %>% 
  st_drop_geometry() %>% 
  summarise(total_length=sum(weight))
#> # A tibble: 1 x 1
#>   total_length
#>            [m]
#> 1     274.6754

# visualize the route
path_edge<-net_main %>% 
  activate("nodes") %>% 
  slice(nps) %>% 
  activate("edges") %>% 
  st_as_sf()

plot(path_edge,add=T, col="green",lwd=3)
#> Warning in plot.sf(path_edge, add = T, col = "green", lwd = 3): ignoring all but
#> the first attribute

^{Created on 2021-05-20 by the reprex package (v2.0.0)}

[zackarno]

One tag along question:

I have a much larger network dataset that I am working with and I have been trying to do the same thing as above basically. Again, I've created a large cost distance matrix, used the work around above to sort out the issue with unique to arguments, and I want to use st_network_paths to first a.) confirm the results of st_network_cost, and b.) visualize a few of the routes. However, every time I do this the distance returned after summing edge lengths returned from st_network_paths is slightly different than the distance in the cost matrix. This being said the difference is only on the order of several hundred meters and the network is made from the entire OSM road network for a country (so its a relatively minor difference). I would very much like to make reprex, but im not sure where to start as the network is quite complex and large. Therefore any suggestions you might have related making this reprex or ideas of how this discrepancy could arise would be much appreciated!

[zackarno]

Here is my best shot at something like a reprex. While I greatly simplified the network to just a path of interest to demonstrate the issue, I could not figure out how to how use datapasta or dput to copy in the sf linestrings as code. Therefore, I have attached the simplified path/network file here- path.zip

In the example below I am showing that st_network_cost is returning a different shortest distance than the sum of all edge lengths returned by st_network_paths. However, I think that they should be exactly the same. Am I missing something or is there anything you can think of that would cause the discrepancy?

library(tidyverse)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(sfnetworks)

path<-read_rds("../Documents/path.rds")
# This is a tiny subset of my overall network. The linestring geometries copied here represent the shortest path on my network between two random points
# I tried my best to simplify the issue by subsetting the network like this, but its still quite a large and complex geometry at least when its written out like this.

# make start and end poitns
start_pt<-st_sfc(st_point(c(427377.9, 3920015)), crs = 32642)
end_pt<-st_sfc(st_point(c(476692.4, 3974736)), crs = 32642)

# plot path with start and end points
plot(st_geometry(path))
plot(st_geometry(start_pt),col="red" ,add=T)
plot(st_geometry(end_pt) , col="green" ,add=T)

net<-sfnetworks::as_sfnetwork(path,directed=FALSE, length_as_weight=TRUE)

# check out the total distance calculated here
result1<-st_network_cost(x = net,from = start_pt,to = end_pt)


# the original network used was already subset to just the shortest path- but lets just pretend we have a bigger
# network as we usually do and want to extract the specific shortest path
shortest_path<-st_network_paths(x = net,from = start_pt,to = end_pt)

# extract node_paths
nps<-shortest_path %>% 
  pull(node_paths) %>% unlist()

# the cost object above returns a different distance than when calculated from the result of st_network_paths
result2<-net %>% 
  activate("nodes") %>% 
  slice(nps) %>% 
  activate("edges") %>% 
  as_tibble() %>% 
  st_drop_geometry() %>% 
  summarise(total_length=sum(weight))

# not sure why these are not the same
result1
#>          [,1]
#> [1,] 106154.4
result2
#> # A tibble: 1 x 1
#>   total_length
#>            [m]
#> 1     116917.3

^{Created on 2021-05-21 by the reprex package (v2.0.0)}

[agila5]

Hi @zackarno! I'm sorry for the late reply; I totally missed this question. Anyway, I think that it's quite difficult to answer your problem without adding a lot of details, but I will try my best. First, I noticed that if you want to compare the two approaches (i.e. st_network_cost and st_network_paths), the easiest solution is to subset the net object using the id(s) of the edges instead of the nodes. For example:

library(sf)
library(tidygraph)
library(sfnetworks)

path <- readRDS("C:/Users/Utente/Desktop/path.RDS")

start_pt <- st_sfc(st_point(c(427377.9, 3920015)), crs = 32642)
end_pt <- st_sfc(st_point(c(476692.4, 3974736)), crs = 32642)

net <- sfnetworks::as_sfnetwork(path, directed = FALSE, length_as_weight = TRUE)

# check out the total distance calculated here
st_network_cost(x = net,from = start_pt, to = end_pt)
#>          [,1]
#> [1,] 106154.4

# estimate shortest path
shortest_path <- st_network_paths(net, start_pt, end_pt)

# extract edges
idx <- shortest_path %>% pull(edge_paths) %>% unlist()

# subset
net %E>% slice(idx) %>% st_as_sf("edges") %>% st_length() %>% sum()
#> 106154.4 [m]

^{Created on 2021-06-13 by the reprex package (v2.0.0)}

I think we can now conclude that this is not a bug in st_network_paths, but it's something related to subsetting sfnetwork objects using the id(s) of the nodes. Then, I noticed that the net object has several duplicated edges and loops:

# packages
library(sf)
library(tidygraph)
library(sfnetworks)

# data
path <- readRDS("C:/Users/Utente/Desktop/path.RDS")

# duplicated edges: 
idx_duplicated <- which(unlist(lengths(st_equals(path)) > 1L))
path[idx_duplicated[c(1, 4, 2, 3)], ]
#> Simple feature collection with 4 features and 0 fields
#> Geometry type: LINESTRING
#> Dimension:     XY
#> Bounding box:  xmin: 435098.7 ymin: 3924449 xmax: 451618 ymax: 3937324
#> Projected CRS: WGS 84 / UTM zone 42N
#> Registered S3 method overwritten by 'cli':
#>   method     from         
#>   print.boxx spatstat.geom
#> # A tibble: 4 x 1
#>                               geometry
#>                       <LINESTRING [m]>
#> 1 (435100.9 3924449, 435098.7 3924460)
#> 2 (435100.9 3924449, 435098.7 3924460)
#> 3   (451618 3937319, 451607.4 3937324)
#> 4   (451618 3937319, 451607.4 3937324)

# loops
net <- sfnetworks::as_sfnetwork(path, directed = FALSE, length_as_weight = TRUE)
par(mar = rep(0, 4))
net %N>% slice(c(12, 11, 10, 68, 67, 43, 42, 9, 6, 5)) %>% plot()

^{Created on 2021-06-13 by the reprex package (v2.0.0)}

The main point here is that subsetting an sfnetwork object using the id(s) of the nodes instead of the id(s) of the edges preserve all edges that go through those nodes (also the duplicated edges and the loops). Hence, I think that there is no "bug" here and you can ignore the small discrepancies. You just need to preprocess the input data (see also the second vignette). It's really easy to fix the problem with duplicated segments (just check the vignette), but I'm not sure how to remove the edges that create a loop (@luukvdmeer can you check here? In any case, I think this is not a big problem).

[luukvdmeer]

Just as a (late) update: since version 0.6 the st_network_cost() does not remove duplicated nodes in the 'to' argument anymore. See https://luukvdmeer.github.io/sfnetworks/news/index.html#sfnetworks-v060-coerde. Hence, the resulting cost matrix will have dimensions nxm where n is the number of provided 'from' nodes and m the number of provided 'to' nodes (i.e. ncol(cost_matrix) == nrow(to_argument)). As you would intuitively expect. Of course, if e.g. 'to' nodes m1 and m2 are the same node (or a different node at the same location) they will of course have the same cost value for each 'from' node ni.