Tweaks for routing #67
Replies: 10 comments 8 replies
-
Same example in download.file("https://github.com/ropensci/stplanr/releases/download/0.6.1/Example.zip", "Example.zip")
unzip("Example.zip")
library(geosphere)
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
library(stplanr)
roads = st_read("Example/Roads/Roads.shp")
#> Reading layer `Roads' from data source `/tmp/Rtmp4Xr6kn/reprex42744568cfd8/Example/Roads/Roads.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 47 features and 1 field
#> geometry type: LINESTRING
#> dimension: XY
#> bbox: xmin: 605246.6 ymin: 7255464 xmax: 615324.9 ymax: 7266851
#> projected CRS: SIRGAS 2000 / UTM zone 22S
points = st_read("Example/Points/Points.shp")
#> Reading layer `Points' from data source `/tmp/Rtmp4Xr6kn/reprex42744568cfd8/Example/Points/Points.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 32 features and 10 fields
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: -49.95058 ymin: -24.80171 xmax: -49.86651 ymax: -24.72459
#> geographic CRS: SIRGAS 2000
# Convert roads to coordinate system of points
roads_trf = st_transform(roads, st_crs(points))
# Convert to points to SpatialPointsDataframe
points_sp = as(points, "Spatial")
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
#> but +towgs84= values preserved
from = c(-49.95058, -24.77502)
to = c(-49.91084, -24.75200)
p = SpatialLinesNetwork(roads_trf, uselonglat = FALSE, tolerance = 0)
r = route_local(p, from, to)
plot(p)
plot(r$geometry, add = TRUE, col = "red", lwd = 5)
plot(points_sp[c(3, 4), ], add = TRUE)
r2 = route_local(sln = p, points[3, ], points[4, ])
plot(r2$geometry, add = TRUE, col = "blue", lwd = 3) Created on 2020-07-07 by the reprex package (v0.3.0) |
Beta Was this translation helpful? Give feedback.
-
The good news: it does seem to work with the large network while download.file("https://github.com/ropensci/stplanr/releases/download/0.6.1/Example.zip", "Example.zip")
unzip("Example.zip")
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
library(sfnetworks)
library(tidygraph)
#>
#> Attaching package: 'tidygraph'
#> The following object is masked from 'package:stats':
#>
#> filter
roads = st_read("Example/Roads/Roads.shp")
#> Reading layer `Roads' from data source `/tmp/Rtmp4Xr6kn/reprex4274da90a2/Example/Roads/Roads.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 47 features and 1 field
#> geometry type: LINESTRING
#> dimension: XY
#> bbox: xmin: 605246.6 ymin: 7255464 xmax: 615324.9 ymax: 7266851
#> projected CRS: SIRGAS 2000 / UTM zone 22S
points = st_read("Example/Points/Points.shp")
#> Reading layer `Points' from data source `/tmp/Rtmp4Xr6kn/reprex4274da90a2/Example/Points/Points.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 32 features and 10 fields
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: -49.95058 ymin: -24.80171 xmax: -49.86651 ymax: -24.72459
#> geographic CRS: SIRGAS 2000
summary(sf::st_geometry_type(roads))
#> GEOMETRY POINT LINESTRING POLYGON
#> 0 0 47 0
#> MULTIPOINT MULTILINESTRING MULTIPOLYGON GEOMETRYCOLLECTION
#> 0 0 0 0
#> CIRCULARSTRING COMPOUNDCURVE CURVEPOLYGON MULTICURVE
#> 0 0 0 0
#> MULTISURFACE CURVE SURFACE POLYHEDRALSURFACE
#> 0 0 0 0
#> TIN TRIANGLE
#> 0 0
# Convert roads to coordinate system of points
roads_trf = st_transform(roads, st_crs(points))
# Convert to points to SpatialPointsDataframe
points_sp = as(points, "Spatial")
#> Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on GRS80 ellipsoid in CRS definition,
#> but +towgs84= values preserved
from = c(-49.95058, -24.77502)
to = c(-49.91084, -24.75200)
net = as_sfnetwork(roads_trf, directed = FALSE) %>%
activate("edges") %>%
mutate(weight = edge_length())
p1 = sf::st_as_sf(data.frame(x = from[1], y = from[2]), coords = c("x", "y"), crs = sf::st_crs(net))
p2 = sf::st_as_sf(data.frame(x = to[1], y = to[2]), coords = c("x", "y"), crs = sf::st_crs(net))
r = net %>%
convert(to_spatial_shortest_paths, p1, p2)
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
r2 = net %>%
convert(to_spatial_shortest_paths, points[4, ], points[4, ])
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
plot(net)
plot(r, col = "red", lwd = 5, add = TRUE)
plot(points_sp[c(3, 4), ], add = TRUE)
#> Error in as.double(y): cannot coerce type 'S4' to vector of type 'double'
plot(r2, add = TRUE, col = "blue", lwd = 3) # with complete roads...
roads = st_read("Example/Complete Roads/complete_roads.shp")
#> Reading layer `complete_roads' from data source `/tmp/Rtmp4Xr6kn/reprex4274da90a2/Example/Complete Roads/complete_roads.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 305265 features and 12 fields
#> geometry type: MULTILINESTRING
#> dimension: XY
#> bbox: xmin: -54.61044 ymin: -26.65364 xmax: -48.09292 ymax: -22.52137
#> geographic CRS: WGS 84
points = st_read("Example/Points/Points.shp")
#> Reading layer `Points' from data source `/tmp/Rtmp4Xr6kn/reprex4274da90a2/Example/Points/Points.shp' using driver `ESRI Shapefile'
#> Simple feature collection with 32 features and 10 fields
#> geometry type: POINT
#> dimension: XY
#> bbox: xmin: -49.95058 ymin: -24.80171 xmax: -49.86651 ymax: -24.72459
#> geographic CRS: SIRGAS 2000
# Convert roads to coordinate system of points
net = as_sfnetwork(roads_trf, directed = FALSE) %>%
activate("edges") %>%
mutate(weight = edge_length())
summary(sf::st_geometry_type(roads))
#> GEOMETRY POINT LINESTRING POLYGON
#> 0 0 0 0
#> MULTIPOINT MULTILINESTRING MULTIPOLYGON GEOMETRYCOLLECTION
#> 0 305265 0 0
#> CIRCULARSTRING COMPOUNDCURVE CURVEPOLYGON MULTICURVE
#> 0 0 0 0
#> MULTISURFACE CURVE SURFACE POLYHEDRALSURFACE
#> 0 0 0 0
#> TIN TRIANGLE
#> 0 0
roads = sf::st_cast(roads, "LINESTRING")
#> Warning in st_cast.sf(roads, "LINESTRING"): repeating attributes for all sub-
#> geometries for which they may not be constant
roads_trf = st_transform(roads, st_crs(points))
system.time({
net = as_sfnetwork(roads_trf, directed = FALSE) %>%
activate("edges") %>%
mutate(weight = edge_length())
})
#> user system elapsed
#> 34.143 0.421 34.784
r = net %>%
convert(to_spatial_shortest_paths, p1, p2)
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> Warning in (function (graph, from, to = V(graph), mode = c("out", "all", : At
#> structural_properties.c:4597 :Couldn't reach some vertices
r2 = net %>%
convert(to_spatial_shortest_paths, points[4, ], points[4, ])
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
plot(net)
plot(r, col = "red", lwd = 5, add = TRUE)
plot(points_sp[c(3, 4), ], add = TRUE)
#> Error in as.double(y): cannot coerce type 'S4' to vector of type 'double'
plot(r2, add = TRUE, col = "blue", lwd = 3) Created on 2020-07-08 by the reprex package (v0.3.0) |
Beta Was this translation helpful? Give feedback.
-
Hi Robin! I think the problem in the first case is that you wrote r2 = net %>%
tidygraph::convert(to_spatial_shortest_paths, points[4, ], points[4, ]) instead of r2 = net %>%
tidygraph::convert(to_spatial_shortest_paths, points[3, ], points[4, ]) Indeed: # packages
library(sf)
#> Linking to GEOS 3.7.1, GDAL 2.2.2, PROJ 4.9.2
library(sfnetworks)
# data
download.file("https://github.com/ropensci/stplanr/releases/download/0.6.1/Example.zip", "Example.zip")
unzip("Example.zip")
# Simple example - sfnetworks approach ---------------------------------------
roads = st_read("Example/Roads/Roads.shp", quiet = TRUE)
points = st_read("Example/Points/Points.shp", quiet = TRUE)
roads_trf = st_transform(roads, st_crs(points))
points_sp = as(points, "Spatial")
# calculate the sfNetwork object
net = as_sfnetwork(roads_trf, directed = FALSE) %>%
activate("edges") %>%
dplyr::mutate(weight = edge_length())
# routing - 1
from = c(-49.95058, -24.77502)
to = c(-49.91084, -24.75200)
p1 = sf::st_as_sf(data.frame(x = from[1], y = from[2]), coords = c("x", "y"), crs = sf::st_crs(net))
p2 = sf::st_as_sf(data.frame(x = to[1], y = to[2]), coords = c("x", "y"), crs = sf::st_crs(net))
r = net %>%
tidygraph::convert(to_spatial_shortest_paths, p1, p2)
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
# routing - 2
r2 = net %>%
tidygraph::convert(to_spatial_shortest_paths, points[3, ], points[4, ])
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
# plot
par(mar = rep(0, 4))
plot(net)
plot(r, col = "red", lwd = 5, add = TRUE)
plot(r2, col = "blue", lwd = 5, add = TRUE) Created on 2020-07-08 by the reprex package (v0.3.0) The second example is more difficult and still it's not working properly with sfnetworks or stplanr. Indeed: # packages
library(sf)
#> Linking to GEOS 3.7.1, GDAL 2.2.2, PROJ 4.9.2
library(sfnetworks)
# data
download.file("https://github.com/ropensci/stplanr/releases/download/0.6.1/Example.zip", "Example.zip")
unzip("Example.zip")
# Simple example - sfnetworks approach ---------------------------------------
roads = st_read("Example/Complete Roads/complete_roads.shp", quiet = TRUE) %>%
st_cast("LINESTRING")
#> Warning in st_cast.sf(., "LINESTRING"): repeating attributes for all sub-
#> geometries for which they may not be constant
points = st_read("Example/Points/Points.shp", quiet = TRUE)
roads_trf = st_transform(roads, st_crs(points))
points_sp = as(points, "Spatial")
# calculate the sfNetwork object
system.time({
net <- as_sfnetwork(roads_trf, directed = FALSE) %>%
activate("edges") %>%
dplyr::mutate(weight = edge_length())
})
#> user system elapsed
#> 44.309 0.484 44.793
# routing - 1
from = c(-49.95058, -24.77502)
to = c(-49.91084, -24.75200)
p1 = sf::st_as_sf(data.frame(x = from[1], y = from[2]), coords = c("x", "y"), crs = sf::st_crs(net))
p2 = sf::st_as_sf(data.frame(x = to[1], y = to[2]), coords = c("x", "y"), crs = sf::st_crs(net))
# It seems there are some problems
st_shortest_paths(net, p1, p2)
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> although coordinates are longitude/latitude, st_nearest_feature assumes that they are planar
#> Warning in (function (graph, from, to = V(graph), mode = c("out", "all", : At
#> structural_properties.c:4597 :Couldn't reach some vertices
#> $vpath
#> $vpath[[1]]
#> + 1/480021 vertex, from 6ef6227:
#> [1] 259784
#>
#>
#> $epath
#> NULL
#>
#> $predecessors
#> NULL
#>
#> $inbound_edges
#> NULL Created on 2020-07-08 by the reprex package (v0.3.0) I have an idea why it doesn't work but I need some time to check it. btw: the red dot that you show in the previous comment is just the starting point instead of the route. |
Beta Was this translation helpful? Give feedback.
-
Thanks @agila5 I've updated the check-boxes also. Great to have a big routing example. |
Beta Was this translation helpful? Give feedback.
This comment was marked as off-topic.
This comment was marked as off-topic.
This comment was marked as off-topic.
This comment was marked as off-topic.
-
Hi @JovaniSouza, if you are still interested in this problem please check ropensci/stplanr#416 (comment). I will implement the new approach as soon as possible. Feedbacks welcome. As soon as we fix this problem, I could prepare a draft for a new vignette that describes sfnetworks + OSM data! |
Beta Was this translation helpful? Give feedback.
-
Hello @agila5. Sorry for the delay in replying, now that I saw it. Very good. If I have any feedback I will let you know. Thanks for letting me know. |
Beta Was this translation helpful? Give feedback.
-
Hi! I think that we can lock the discussion here since the only missing point that should be discussed is the same as #50. Reprex to show the suggested (at least from my point of view) approach: # packages
library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
library(tidygraph)
library(sfnetworks)
# data
download.file("https://github.com/ropensci/stplanr/releases/download/0.6.1/Example.zip", "Example.zip")
unzip("Example.zip")
# read-in and transform data
roads = st_read("Example/Complete Roads/complete_roads.shp", quiet = TRUE) %>%
st_cast("LINESTRING")
#> Warning in st_cast.sf(., "LINESTRING"): repeating attributes for all sub-
#> geometries for which they may not be constant
points = st_read("Example/Points/Points.shp", quiet = TRUE)
roads_trf = st_transform(roads, st_crs(points))
points_sp = as(points, "Spatial")
# create sfnetwork structure
roads_sfnetwork <- as_sfnetwork(roads, directed = FALSE, length_as_weight = TRUE)
# apply subdivision morpher
roads_sfnetwork <- convert(roads_sfnetwork, to_spatial_subdivision)
#> Warning: to_spatial_subdivision assumes attributes are constant over geometries
# estimate shortest path
roads_sfnetwork %>%
convert(
to_spatial_shortest_paths,
from = st_sfc(st_point(c(-49.95058, -24.77502)), crs = st_crs(roads_sfnetwork)),
to = st_sfc(st_point(c(-49.91084, -24.75200)), crs = st_crs(roads_sfnetwork))
)
#> 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
#> # A sfnetwork with 28 nodes and 27 edges
#> #
#> # CRS: WGS 84
#> #
#> # An unrooted tree with spatially explicit edges
#> #
#> # Node Data: 28 x 2 (active)
#> # Geometry type: POINT
#> # Dimension: XY
#> # Bounding box: xmin: -49.94907 ymin: -24.79465 xmax: -49.88282 ymax:
#> # -24.75606
#> geometry .tidygraph_node_index
#> <POINT [°]> <int>
#> 1 (-49.94291 -24.79282) 6049
#> 2 (-49.93858 -24.78435) 6112
#> 3 (-49.93878 -24.7764) 6113
#> 4 (-49.94312 -24.79364) 6940
#> 5 (-49.94282 -24.79284) 6941
#> 6 (-49.94317 -24.79409) 7426
#> # ... with 22 more rows
#> #
#> # Edge Data: 27 x 17
#> # Geometry type: LINESTRING
#> # Dimension: XY
#> # Bounding box: xmin: -49.94907 ymin: -24.79557 xmax: -49.88282 ymax:
#> # -24.75606
#> from to OBJECTID osm_id code fclass name ref oneway maxspeed layer
#> <int> <int> <dbl> <fct> <int> <fct> <fct> <fct> <fct> <int> <dbl>
#> 1 2 3 1381 18712~ 5115 terti~ Estr~ <NA> B 0 0
#> 2 4 5 1570 18676~ 5114 secon~ Aven~ <NA> F 40 0
#> 3 6 7 1681 18676~ 5114 secon~ Aven~ <NA> F 40 0
#> # ... with 24 more rows, and 6 more variables: bridge <fct>, tunnel <fct>,
#> # Shape_Leng <dbl>, geometry <LINESTRING [°]>, weight [m],
#> # .tidygraph_edge_index <int> Created on 2020-12-14 by the reprex package (v0.3.0) |
Beta Was this translation helpful? Give feedback.
-
@agila5 @Robinlovelace @loreabad6 Do you think that in I can imagine that for spatial shortest paths, as @agila5 pointed out, it is rarely intended to calculate paths without edge weights. Therefore it might make more sense to always use the edge length as default in that case. Hence, when the edges don't have a weight column and the weights argument is set to NULL, Good idea or too confusing? |
Beta Was this translation helpful? Give feedback.
-
Opening this issue to document some observations about routing with
sfnetworks
based on ropensci/stplanr#412It may be worth splitting-out sub issues but it's a great example of routing on a big network so seems worth documenting in one place.
sfnetworks
does not seem to generate a result for the second route shown below whilestplanr
does, not sure why (see next comment for reproducible solution instplanr
)structural_properties.c:4597 :Couldn't reach some vertices
, also affectingstplanr
.Reproducible example:
Created on 2020-07-07 by the reprex package (v0.3.0)
Beta Was this translation helpful? Give feedback.
All reactions