Sf: st_node bug

[buxinator]

Hi,
i have a Shapefile (read with st_read) with missing Nodes at Intersections of lines. For that i want to use st_node("file") which should exactly solve this problem but the output does not change anything at all. Any idea why?

[Robinlovelace]

A reproducible example stating what you've done and what your would expect may help, as with all software issues. (Vested interest: I've not used this function but imagine it could be useful for my work so would like to see it action.)

[buxinator]

ok i try so be more specific:
for example i have two lines in one shapefile and they cross somewhere. The Problem is that there is no node at the Crossing. I need a node there because I want to use this lines for routing and the routing algorithm only jumps to another line if there is a node.
The description of st_node(x) is: "st_node add nodes to linear geometries at intersections without node"

Input:

Expected output:

What I've done is loading the shape file <- st_read("lines.shp") and file1 <- st_node(file)
but it is not adding any node to the lines

[Robinlovelace]

Interesting. The behaviour is reproduced below - I would expect points to be returned but am not sure if that's the intended behaviour or not having never used the function.

library(sf)
#> Linking to GEOS 3.5.1, GDAL 2.2.2, proj.4 4.9.2
l1 = st_linestring(matrix(c(0, 0, 1, 1), ncol = 2, byrow = TRUE))
l2 = st_linestring(matrix(c(0, 1, 1, 0), ncol = 2, byrow = TRUE))
plot(l1)
plot(l2, add = TRUE)

l = st_sfc(l1, l2)
st_node(l)
#> Geometry set for 2 features 
#> geometry type:  MULTILINESTRING
#> dimension:      XY
#> bbox:           xmin: 0 ymin: 0 xmax: 1 ymax: 1
#> epsg (SRID):    NA
#> proj4string:    NA
#> MULTILINESTRING ((0 0, 1 1))
#> MULTILINESTRING ((0 1, 1 0))

Would you expect just a point to be returned?

[buxinator]

No, a node and a point are different things. A point is a new geometry and a node belongs to the line. I made the nodes visible with QGIS, if you "toggle editing" you see the nodes.

thats yout l:

and thats st_node(l):

so nothing changes but you would expect a new node at intersection. In my case I need this node because if i want to route for example from the bottom left to the bottom right it is not possible without the node at intersection

[rsbivand]

Yes, if you try:

nl <- st_node(l)
plot(nl)
plot(st_cast(nl, "MULTIPOINT"), add=TRUE)

you see a different behaviour from:

library(rgeos)
gnl <- gNode(as(l, "Spatial"))
plot(st_as_sfc(gnl))
plot(st_cast(st_as_sfc(gnl), "MULTIPOINT"), add=TRUE)

which is noding as expected.

[buxinator]

gNode also Unions all lines and I loose all data behind the lines? Am I correct?

[Robinlovelace]

Bingo - thanks for showing the expected output Roger. For people too lazy to run the above code (often me!) see below the output for that final line which is as the original would expect. @buxinator You could for now use rgeos::gNode(). It seems like a bug that it doesn't work in sf.

@buxinator regarding your question: you don't loose any data if you create a new object and casting it converts it into points as required.

[Robinlovelace]

[rsbivand]

Well, there are two unnoded lines input, and four lines split at the node in the object returned (as parts of a MULTILINESTRING feature), so:

st_cast(st_as_sfc(gnl), "LINESTRING")

and

st_touches(st_cast(st_as_sfc(gnl), "LINESTRING"), l, sparse = FALSE)

may help disambiguate. The output line strings cannot cross each other, but from predicates on before/after, you could identify which output lines constitute your input line. My guess - unchecked - is that st_node() would work noding a MULTILINESTRING in a single feature.

Checked:

st_combine(l)
st_node(st_combine(l))
plot(l)
plot(st_cast(st_node(st_combine(l)), "MULTIPOINT"), add=TRUE)

and

nl <- st_node(st_combine(l))
st_touches(st_cast(nl, "LINESTRING"), l, sparse=FALSE)

[buxinator]

st_node is not working on MULTILINESTRING, i tried.

@Robinlovelace I don't get what you mean? If i have two specific IDs for those lines, how can keep them after using gNode?

[rsbivand]

st_node needs all the intersecting LINESTRING objects in a single MULTILINESTRING object, gNode doesn't, but puts all the output in a single object. Use predicates between the input object and the output to find out which is which. The assumption until now has been that st_node and gNode are for cleaning LINESTRING objects that cross each other before merging lines or building polygons. Your use is different, but you can use predicates to get what you need.

[Robinlovelace]

Below is an indication of where I think we're at with this - it can be done but the lines must by in a multilinestring with each LINESTRING within that representing different segments on the line network (I imagine you're using a route network - have your tried the stplanr function SpatialLInesNetwork() @buxinator?).

One perhaps naive question: how can I get st_cast() to output all the linestrings in the code below?

library(sf)
#> Linking to GEOS 3.5.1, GDAL 2.2.2, proj.4 4.9.2
# > Linking to GEOS 3.5.1, GDAL 2.2.2, proj.4 4.9.2
l1 = st_linestring(matrix(c(0, 0, 1, 1), ncol = 2, byrow = TRUE))
l2 = st_linestring(matrix(c(0, 1, 1, 0), ncol = 2, byrow = TRUE))
l = st_multilinestring(list(l1, l2))
plot(l, lwd = 5, col = "grey")
ln = st_node(l)
ln_linestring = st_cast(ln, "LINESTRING", group_or_split = TRUE)
#> Warning in st_cast.MULTILINESTRING(ln, "LINESTRING", group_or_split =
#> TRUE): keeping first linestring only
plot(ln_linestring, add = TRUE, col = 1:length(ln_linestring))  # was expecting all lines in there
p = st_cast(ln, "MULTIPOINT")
p
#> MULTIPOINT (0 0, 0.5 0.5, 0.5 0.5, 1 1, 0 1, 0.5 0.5, 0.5 0.5, 1 0)
plot(p, add = TRUE)

[Robinlovelace]

@buxinator regarding your question, fair point: there is no 1:1 relationship between nodes and edges on a road/graph network. I have been looking for a similar solution in the igraph package, so far to no avail so my first port of call would be st_join() to rejoin attributes to those points.

[Robinlovelace]

Update, I can see the answer to my question about st_cast() that the solution is in an answer to an issue raised by me in #114. The solution is in the below updated code chunk:

library(sf)
#> Linking to GEOS 3.5.1, GDAL 2.2.2, proj.4 4.9.2
l1 = st_linestring(matrix(c(0, 0, 1, 1), ncol = 2, byrow = TRUE))
l2 = st_linestring(matrix(c(0, 1, 1, 0), ncol = 2, byrow = TRUE))
l = st_multilinestring(list(l1, l2))
plot(l, lwd = 5, col = "grey")
ln = st_node(l)
ln_linestring = st_cast(st_sfc(ln), "LINESTRING", group_or_split = TRUE)
plot(ln_linestring, add = TRUE, col = 1:length(ln_linestring))
p = st_cast(ln, "MULTIPOINT")
p
#> MULTIPOINT (0 0, 0.5 0.5, 0.5 0.5, 1 1, 0 1, 0.5 0.5, 0.5 0.5, 1 0)
plot(p, add = TRUE)
p_many = st_cast(st_sfc(p), "POINT")
plot(p_many, col = 1:length(p_many), add = TRUE)

[buxinator]

@Robinlovelace I am using OSRM to set a local server the library(osrm) for R

[buxinator]

Thanks for all your help! If I find a satisfying solution for my problem, I will give an update

[Robinlovelace]

Great stuff - let us know how you get on.

[edzer]

Can someone here please summarise what the bug (see title) in current st_node is?

[Robinlovelace]

I'm not sure there is one after the latest code chunk which reproduces the behaviour of gNode(). Can you confirm that st_node() should only work on multilinestring geometries? Just checked the docs and currently this is not stated:

st_node add nodes to linear geometries at intersections without node

Perhaps updated documentation would solve the issue.

[edzer]

No: taken from the example section:

> st_node(st_linestring(rbind(c(0,0), c(1,1), c(0,1), c(1,0), c(0,0))))
# MULTILINESTRING ((0 0, 0.5 0.5), (0.5 0.5, 1 1, 0 1, 0.5 0.5), (0.5 0.5, 1 0, 0 0))

but also

> st_node(st_multilinestring(list(rbind(c(0,0), c(1,1), c(0,1), c(1,0), c(0,0)))))
# MULTILINESTRING ((0 0, 0.5 0.5), (0.5 0.5, 1 1, 0 1, 0.5 0.5), (0.5 0.5, 1 0, 0 0))

I think the point is that st_node works only on individual line geometries; maybe we should make that more clear?

[Robinlovelace]

I think the point is that st_node works only on individual line geometries; maybe we should make that more clear?

I agree!

[Robinlovelace]

A related question: documentation for many sf functions is aggregated into a 'metadocumentation' page, geos_unary being an example. I think it would help people find out about functions if there were more documentation pages that were about a specific function. Any objections to some documentations being split-out? Broader issues like unary operations could be mentioned with reference to vignettes/vignette sections.

[edzer]

Do you think it is more helpful to have 12 different help pages for the 12 functions in that page, or to have at least a few details about each of them in the Details section?

[Robinlovelace]

I'm confident 12 separate pages will be more useful for most people, especially for beginners.

[Robinlovelace]

Caveat: If they are well written and link to a central theme or a metadocumentation page such as geos_unary.