Demonstrate routing on large network represented as an sf object in memory #412

JovaniSouza · 2020-07-07T20:59:11Z

I would like to know why it gives error when I use a shapefile with many roads. I am using the code below to make the path between the points. When I use the Roads.shp file, it works correctly. However when I use the complete_roads.shp file, which are several roads, it doesn't work, it gives an error. Could you help me solve it? Shapefiles can be downloaded from the following website: https://drive.google.com/file/d/1krpRiU2RBPWQiVz0KOaHn0JY5zCAY2eH/view?usp=sharing

Thank you so much!

library(geosphere)
library(sf)
library(stplanr)

roads <- st_read("Example/Roads/Roads.shp")
points <- st_read("Example/Points/Points.shp")

# 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")

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)

Image generated by code above

The text was updated successfully, but these errors were encountered:

Robinlovelace · 2020-07-07T21:31:16Z

I can reproduce the error and am not sure what's causing it:

setwd("/tmp")
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/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/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)

# with complete roads...
roads <- st_read("Example/Complete Roads/complete_roads.shp")
#> Reading layer `complete_roads' from data source `/tmp/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/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)
#> Error in coord_matches_sf(as.matrix(nodecoords %>% dplyr::select(.data$X, : Mat::rows(): indices out of bounds or incorrectly used

^{Created on 2020-07-07 by the reprex package (v0.3.0)}

Robinlovelace · 2020-07-07T21:32:41Z

The coord_matches_sf() function was written by @richardellison. Any ideas Richard?

Robinlovelace · 2020-07-07T21:56:35Z

Update: the geometry in the complete roads dataset is not LINESTRING.

Robinlovelace · 2020-07-07T22:00:51Z

@JovaniSouza can you try running this line before creating the spatial network?:

roads = sf::st_cast(roads, "LINESTRING")

Robinlovelace · 2020-07-07T23:12:47Z

Update: I can get this to work with the sfnetworks package, but it's not fast:

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)}

JovaniSouza · 2020-07-08T01:57:44Z

Thansk for reply Robin. What do you suggest, use the sfnetworks package? I tried to install by remotes :: install_github ("luukvdmeer / sfnetworks"), but it gaves error.

agila5 · 2020-07-08T07:43:01Z

Hi @JovaniSouza, please open a new issue in the sfnetworks repository and add as many details as possible (e.g. OS, R version, sf version, igraph version, tidygraph version, dplyr version and so on)

Robinlovelace · 2020-07-08T08:29:15Z

It's a great use case. Are you happy for me to close the issue here as a solution has been found in another package?

JovaniSouza · 2020-07-08T12:58:52Z

Thanks @agila5, I I will do this and @Robinlovelace thank you very much for your help, I would like to use your package, I found it very interesting, unfortunately it didn't work, so I will use the package you suggested to me. If you solve this problem please let me know. Thank you!

Robinlovelace · 2020-07-08T15:03:37Z

I will keep this issue open to track progress on interfaces to sfnetworks and other packages to make this happen. One other suggestion @JovaniSouza, have you tried using route_dodgr()?

JovaniSouza · 2020-07-08T15:06:56Z

No @Robinlovelace,
I haven't tried, is it from your package?

Robinlovelace · 2020-07-08T15:10:28Z

I haven't tried, is it from your package?

Yes, give it a try and see the documentation in the link in my previous comment.

agila5 mentioned this issue Jul 8, 2020

rnet_breakup_vertices is way too slow for real network data #416

Closed

Robinlovelace changed the title ~~Errorrr route_local~~ Demonstrate routing on large network represented as an sf object in memory Jul 9, 2020

Robinlovelace added the enhancement label Jul 9, 2020

Robinlovelace added this to the 1.0 milestone Jul 9, 2020

agila5 mentioned this issue Oct 5, 2020

Morpher that breaks edge if it shares point with other edge luukvdmeer/sfnetworks#73

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Demonstrate routing on large network represented as an sf object in memory #412

Demonstrate routing on large network represented as an sf object in memory #412

JovaniSouza commented Jul 7, 2020 •

edited

Loading

Robinlovelace commented Jul 7, 2020 •

edited

Loading

Robinlovelace commented Jul 7, 2020 •

edited

Loading

Robinlovelace commented Jul 7, 2020

Robinlovelace commented Jul 7, 2020

Robinlovelace commented Jul 7, 2020

JovaniSouza commented Jul 8, 2020

agila5 commented Jul 8, 2020

Robinlovelace commented Jul 8, 2020

JovaniSouza commented Jul 8, 2020

Robinlovelace commented Jul 8, 2020

JovaniSouza commented Jul 8, 2020

Robinlovelace commented Jul 8, 2020

Demonstrate routing on large network represented as an sf object in memory #412

Demonstrate routing on large network represented as an sf object in memory #412

Comments

JovaniSouza commented Jul 7, 2020 • edited Loading

Robinlovelace commented Jul 7, 2020 • edited Loading

Robinlovelace commented Jul 7, 2020 • edited Loading

Robinlovelace commented Jul 7, 2020

Robinlovelace commented Jul 7, 2020

Robinlovelace commented Jul 7, 2020

JovaniSouza commented Jul 8, 2020

agila5 commented Jul 8, 2020

Robinlovelace commented Jul 8, 2020

JovaniSouza commented Jul 8, 2020

Robinlovelace commented Jul 8, 2020

JovaniSouza commented Jul 8, 2020

Robinlovelace commented Jul 8, 2020

JovaniSouza commented Jul 7, 2020 •

edited

Loading

Robinlovelace commented Jul 7, 2020 •

edited

Loading

Robinlovelace commented Jul 7, 2020 •

edited

Loading