Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

calc_network_catchment does not seem to show all paths #261

Closed
stevenpandrews opened this issue Jun 5, 2018 · 14 comments
Closed

calc_network_catchment does not seem to show all paths #261

stevenpandrews opened this issue Jun 5, 2018 · 14 comments

Comments

@stevenpandrews
Copy link

stevenpandrews commented Jun 5, 2018 •
edited

The calc_network_catchment function seems to only find shortest paths based on entire features. It won't break features into smaller parts. If a portion of a line would be within 800m, but the end segment is farther than 800m away, the entire feature won't be selected.

It also seems to miss some paths that are are within the maximpedance, but are not the fastest path.

I ran the code mostly as shown in the example:

library(stplanr)
library(sp)
library(rgdal)

data_dir <- system.file("extdata", package = "stplanr")
unzip(file.path(data_dir, 'smallsa1.zip'), exdir=tempdir())
unzip(file.path(data_dir, 'testcycleway.zip'), exdir=tempdir())
unzip(file.path(data_dir, 'sydroads.zip'), exdir=tempdir())
sa1income <- readOGR(tempdir(),"smallsa1")
testcycleway <- readOGR(tempdir(),"testcycleway")
sydroads <- readOGR(tempdir(),"roads")
sydnetwork <- SpatialLinesNetwork(sydroads)
network <- calc_network_catchment(
  sln = sydnetwork,
  polygonlayer = sa1income,
  targetlayer = testcycleway,
  calccols = c('Total'),
  maximpedance = 800,
  distance = 10,
  projection = 'austalbers',
  dissolve = TRUE
)

mapview::mapview(sydroads) + mapview::mapview(network, color = "green") + mapview::mapview(testcycleway, color = "red")

I changed the distance to 10 to make it more clear what's going on. It seems like the circled neighborhood should have more coverage. Because the nodes on either end of each street are already part of a different shortest path, the roads never get selected.
image

Is there a way to get all paths that are less than 800m away? I can probably figure out a way to break the line segments into smaller pieces.
@stevenpandrews stevenpandrews changed the title calc_network_catchment does not seem to be accurate calc_network_catchment does not seem to show all paths Jun 5, 2018
@Robinlovelace
Copy link
Member

Interesting and I confirm the issue, with distance = 50 and

mapview::mapview(sydroads) +
  mapview::mapview(network, col.regions = "green") +
  mapview::mapview(testcycleway, color = "red", lwd = 5)

See what you mean about the lack of coverage in that neighbourhood. Also would you expect more coverage around the middle section of the left line?

image

Note: I didn't write the function or define its purpose so may not be the best person to fix this issue. Thoughts @richardellison ?

Cheers for reporting this with a decent reproducible example in any case.

@stevenpandrews
Copy link
Author

stevenpandrews commented Jun 6, 2018 •
edited

Yes, I'd expect there to be more coverage in a few places on the left side and a few roads on the right side too. I might dig into the igraph package to see if there are other functions other than igraph::get.shortest.paths that would do the trick.

@richardellison
Copy link
Collaborator

richardellison commented Jun 7, 2018 via email

@mpadge
Copy link
Member

mpadge commented Jun 7, 2018

OSM mandates that all intersections exist as single points, so dodgr could help here by contracting the network graph to intersections only, then rebuilding into sf-line segments that terminate at every intersection. Let me know if that might help, and I'll paste code ...

@Robinlovelace
Copy link
Member

@mpadge yes please do - that could solve this problem and perhaps be part of a wider function like route_clean() or route_rebuild(). Many thanks + think now the issue can be seen as a question of input data preparation.

@stevenpandrews
Copy link
Author

stevenpandrews commented Jun 7, 2018 •
edited

I think there are two issues here: 1) Long lines that don't break on intersections, which dodgr might help with, and 2) the idea that an edge that isn't on a shortest path, but is within the specified impedance, won't be included in the output of the igraph::get.shortest.paths function. For example, given this graph:

library(magrittr)
library(igraph)

# Make the dataframe
df <- data.frame(
  A = c("A", "A", "C", "B", "C", "D"), 
  B = c("B", "C", "D", "D", "B", "E"))

# Make the graph and set the edge weights
df.g <- graph.data.frame(d = df, directed = FALSE) %>% 
  set_edge_attr("weight", value = c(0.10, 0.10, 0.40, 0.30, 0.20, 1.50))

plot(df.g, vertex.label = V(df.g)$name)

# Run a few examples
shortest.paths(df.g, "A")
get.shortest.paths(df.g, "A", output = "epath")
all_simple_paths(df.g, "A")

and suppose we were trying to find all of the edges within 1 unit from "A" by finding the shortest path to each node, we'd find out that AB, AC, ABD, ABDE are all the shortest paths (from get.shortest.paths(df.g, "A", output = "epath")). So even though ABC (and ACB) are within 1 unit, it's never covered by a shortest path. Same goes for CD. We would throw out ABDE (because it is more than 1 unit away when we get to E). This is understandable given that we're using the ends of lines as nodes.

It seems like we need a list of all the paths the algorithm searched on its way to the shortest path and choose all of them that are shorter than the impedance parameter. Using igraph::all_simple_paths() works on this simple example, but seems to run slowly on a very pared back road network. I was also having trouble getting the distance of all the paths out of all_simple_paths(df.g, "A").

@Robinlovelace
Copy link
Member

Good clarification @sandrews-CTPS - agree 1) is a data prep issue and 2) is more an issue associated with the function. I think we can make progress on both.

@richardellison
Copy link
Collaborator

richardellison commented Jun 7, 2018 via email

@mpadge
Copy link
Member

mpadge commented Jun 8, 2018

The dodgr code to strip a network back to line segments which terminate at every junction point is simply

> streetnet <- dodgr::dodgr_streetnet ("glebe sydney")
> nrow (streetnet)
[1] 690
> streetnet2 <- weight_streetnet (streetnet, wt_profile = 1) %>% # wt_profile irrelevant here
    dodgr::dodgr_to_sf (graph)

The streetnet2 object returned from dodgr_to_sf() is a list of $data and $geoms, where the latter is an sfc() object. Assembly is left to the user to avoid sf dependency, but is simply:

> streetnet2 <- sf::st_sf (data = g$dat, geometry = g$geoms)
> nrow (streetnet2)
[1] 2264

Cutting at juntions increases number of LINESTRING objects by just over threefold.

@Robinlovelace
Copy link
Member

Thanks for that - seems dodgr can add new nodes where they are needed, breaking the network down into its component segments. That's awesome - something I hadn't grasped about the package.

It tooks some wrestling to get what I think are comparative results based on OSM and dodr networks but eventually it seems to have paid-off as catchment based on the dodgr one seems way more detailed (see below).

@richardellison should we recommend users decompose their networks before doing SpatialLinesNetwork()?

@mpadge I got an error message with the below, should I add an issue to https://github.com/ATFutures/dodgr/issues ?:

weight_streetnet(sydroads_mini, wt_profile = 1) %>% dodgr_to_sf()
Error in rcpp_aggregate_to_sf(net, gc$graph, gc$edge_map) : 
  CHAR() can only be applied to a 'CHARSXP', not a 'NULL'

@sandrews-CTPS what does this solve your problem for now? Input into the documentation very welcome - give me a shout if you'd like pointers on that - but the input so far is appreciated and should (eventually) lead to improvements in stplanr in any case!

library(stplanr)
library(dodgr)
library(sf)
#> Linking to GEOS 3.5.1, GDAL 2.1.2, proj.4 4.9.3
library(rgdal)
#> Loading required package: sp
#> rgdal: version: 1.2-20, (SVN revision 725)
#>  Geospatial Data Abstraction Library extensions to R successfully loaded
#>  Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
#>  Path to GDAL shared files: /usr/share/gdal/2.1
#>  GDAL binary built with GEOS: TRUE 
#>  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
#>  Path to PROJ.4 shared files: (autodetected)
#>  Linking to sp version: 1.2-7
library(sp)

data_dir <- system.file("extdata", package = "stplanr")
unzip(file.path(data_dir, "smallsa1.zip"), exdir = tempdir())
unzip(file.path(data_dir, "testcycleway.zip"), exdir = tempdir())
unzip(file.path(data_dir, "sydroads.zip"), exdir = tempdir())
sa1income <- readOGR(tempdir(), "smallsa1")
#> OGR data source with driver: ESRI Shapefile 
#> Source: "/tmp/Rtmpjqg5nc", layer: "smallsa1"
#> with 638 features
#> It has 19 fields
#> Integer64 fields read as strings:  OBJECTID_1
testcycleway <- readOGR(tempdir(), "testcycleway")
#> OGR data source with driver: ESRI Shapefile 
#> Source: "/tmp/Rtmpjqg5nc", layer: "testcycleway"
#> with 2 features
#> It has 2 fields
#> Integer64 fields read as strings:  id
testcycleway_sf <- read_sf(tempdir(), "testcycleway")
sydroads <- readOGR(tempdir(), "roads")
#> OGR data source with driver: ESRI Shapefile 
#> Source: "/tmp/Rtmpjqg5nc", layer: "roads"
#> with 4451 features
#> It has 7 fields
#> Integer64 fields read as strings:  osm_id
sydroads_sf <- as(sydroads, "sf")
sydroads_mini <- sydroads_sf[testcycleway_sf[1, ], , op = st_is_within_distance, 
  dist = 100]

# get network from dodgr: what's the difference?
streetnet <- dodgr::dodgr_streetnet("glebe sydney")
streetnet_mini = streetnet[sydroads_mini, ]
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
# setdiff(names(streetnet_mini), names(sydroads_mini)) # quite a lot!

# characterise street network:
l = st_length(streetnet_mini)
w = l/mean(l)
plot(streetnet_mini["osm_id"], lwd = w)

# clean the street network:
nrow(streetnet_mini)
#> [1] 171
streetnet_mini$highway = "tertiary"  # attempt to make the next bit work...
graph <- weight_streetnet(streetnet_mini)
streetnet_dodgr <- dodgr::dodgr_to_sf(graph) %>% st_sf()
nrow(streetnet_dodgr)/nrow(streetnet_mini)  # 4 fold increase
#> [1] 4.02924

l = st_length(streetnet_dodgr)
w = l/mean(l)
plot(st_geometry(streetnet_dodgr), lwd = w, col = 1:nrow(streetnet_dodgr))

# sln1 <- SpatialLinesNetwork(sydroads_mini) # fails
sln1 <- SpatialLinesNetwork(sydroads)
sln2 <- SpatialLinesNetwork(as(streetnet_dodgr, "Spatial"))

network1 <- calc_network_catchment(sln = sln1, polygonlayer = sa1income, targetlayer = testcycleway, 
  calccols = c("Total"), maximpedance = 800, distance = 10, projection = "austalbers", 
  dissolve = TRUE)
#> Warning in igraph::get.shortest.paths(sln@g, x, output = "epath"): At
#> structural_properties.c:4597 :Couldn't reach some vertices

#> Warning in igraph::get.shortest.paths(sln@g, x, output = "epath"): At
#> structural_properties.c:4597 :Couldn't reach some vertices

#> Warning in igraph::get.shortest.paths(sln@g, x, output = "epath"): At
#> structural_properties.c:4597 :Couldn't reach some vertices

#> Warning in igraph::get.shortest.paths(sln@g, x, output = "epath"): At
#> structural_properties.c:4597 :Couldn't reach some vertices

#> Warning in igraph::get.shortest.paths(sln@g, x, output = "epath"): At
#> structural_properties.c:4597 :Couldn't reach some vertices

#> Warning in igraph::get.shortest.paths(sln@g, x, output = "epath"): At
#> structural_properties.c:4597 :Couldn't reach some vertices

#> Warning in igraph::get.shortest.paths(sln@g, x, output = "epath"): At
#> structural_properties.c:4597 :Couldn't reach some vertices

#> Warning in igraph::get.shortest.paths(sln@g, x, output = "epath"): At
#> structural_properties.c:4597 :Couldn't reach some vertices

#> Warning in igraph::get.shortest.paths(sln@g, x, output = "epath"): At
#> structural_properties.c:4597 :Couldn't reach some vertices

#> Warning in igraph::get.shortest.paths(sln@g, x, output = "epath"): At
#> structural_properties.c:4597 :Couldn't reach some vertices

#> Warning in igraph::get.shortest.paths(sln@g, x, output = "epath"): At
#> structural_properties.c:4597 :Couldn't reach some vertices

network2 <- calc_network_catchment(sln = sln2, polygonlayer = sa1income, targetlayer = testcycleway[1, 
  ], calccols = c("Total"), maximpedance = 800, distance = 10, projection = "austalbers", 
  dissolve = TRUE)

# network1 = network1[network2, ]

plot(network2)
plot(network1, col = "red", add = T)

Session info 
devtools::session_info()
#> Session info -------------------------------------------------------------
#>  setting  value                       
#>  version  R version 3.5.0 (2018-04-23)
#>  system   x86_64, linux-gnu           
#>  ui       X11                         
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  tz       Etc/UTC                     
#>  date     2018-06-08
#> Packages -----------------------------------------------------------------
#>  package      * version    date      
#>  assertthat     0.2.0      2017-04-11
#>  backports      1.1.2      2017-12-13
#>  base         * 3.5.0      2018-05-20
#>  bindr          0.1.1      2018-03-13
#>  bindrcpp       0.2.2      2018-03-29
#>  class          7.3-14     2015-08-30
#>  classInt       0.2-3      2018-04-16
#>  compiler       3.5.0      2018-05-20
#>  curl           3.2        2018-03-28
#>  datasets     * 3.5.0      2018-05-20
#>  DBI            1.0.0      2018-05-02
#>  devtools       1.13.5     2018-02-18
#>  digest         0.6.15     2018-01-28
#>  dodgr        * 0.1.0.099  2018-06-08
#>  dplyr          0.7.5      2018-05-19
#>  e1071          1.6-8      2017-02-02
#>  evaluate       0.10.1     2017-06-24
#>  foreign        0.8-70     2018-04-23
#>  formatR        1.5        2017-04-25
#>  geosphere      1.5-7      2017-11-05
#>  glue           1.2.0      2017-10-29
#>  graphics     * 3.5.0      2018-05-20
#>  grDevices    * 3.5.0      2018-05-20
#>  grid           3.5.0      2018-05-20
#>  htmltools      0.3.6      2017-04-28
#>  httr           1.3.1      2017-08-20
#>  igraph         1.2.1      2018-03-10
#>  jsonlite       1.5        2017-06-01
#>  knitr          1.20       2018-02-20
#>  lattice        0.20-35    2017-03-25
#>  lubridate      1.7.4      2018-04-11
#>  lwgeom         0.1-4      2018-01-28
#>  magrittr       1.5        2014-11-22
#>  maptools       0.9-2      2017-03-25
#>  memoise        1.1.0      2017-04-21
#>  methods      * 3.5.0      2018-05-20
#>  mime           0.5        2016-07-07
#>  openxlsx       4.0.17     2017-03-23
#>  osmdata        0.0.7      2018-05-17
#>  pillar         1.2.2      2018-04-26
#>  pkgconfig      2.0.1      2017-03-21
#>  purrr          0.2.4      2017-10-18
#>  R.methodsS3    1.7.1      2016-02-16
#>  R.oo           1.22.0     2018-04-22
#>  R.utils        2.6.0      2017-11-05
#>  R6             2.2.2      2017-06-17
#>  raster         2.6-7      2017-11-13
#>  rbenchmark     1.0.0      2012-08-30
#>  Rcpp           0.12.17    2018-05-18
#>  RcppParallel   4.4.0      2018-03-02
#>  rgdal        * 1.2-20     2018-05-07
#>  rgeos          0.3-26     2017-10-31
#>  rlang          0.2.0      2018-02-20
#>  rmarkdown      1.9        2018-03-01
#>  rprojroot      1.3-2      2018-01-03
#>  rvest          0.3.2      2016-06-17
#>  sf           * 0.6-3      2018-05-17
#>  sp           * 1.3-1      2018-06-05
#>  spData         0.2.8.9    2018-05-25
#>  spDataLarge    0.2.6.5    2018-05-25
#>  stats        * 3.5.0      2018-05-20
#>  stplanr      * 0.2.4.9000 2018-05-25
#>  stringi        1.2.2      2018-05-02
#>  stringr        1.3.1      2018-05-10
#>  tibble         1.4.2      2018-01-22
#>  tidyselect     0.2.4      2018-02-26
#>  tools          3.5.0      2018-05-20
#>  udunits2       0.13       2016-11-17
#>  units          0.5-1      2018-01-08
#>  utils        * 3.5.0      2018-05-20
#>  withr          2.1.2      2018-03-15
#>  xml2           1.2.0      2018-01-24
#>  yaml           2.1.19     2018-05-01
#>  source                                
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  local                                 
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  local                                 
#>  CRAN (R 3.5.0)                        
#>  local                                 
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  Github (ATFutures/dodgr@51d30fd)      
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  local                                 
#>  local                                 
#>  local                                 
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  local                                 
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  cran (@1.0.0)                         
#>  CRAN (R 3.5.0)                        
#>  cran (@4.4.0)                         
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  cran (@1.3-1)                         
#>  Github (nowosad/spData@7b53933)       
#>  Github (nowosad/spDataLarge@bc058ad)  
#>  local                                 
#>  Github (robinlovelace/stplanr@276c68e)
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  local                                 
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  local                                 
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)                        
#>  CRAN (R 3.5.0)

@mpadge
Copy link
Member

mpadge commented Jun 8, 2018

@Robinlovelace yes, please open an issue in dodgr with reprex

@stevenpandrews
Copy link
Author

@Robinlovelace I'm happy that "long line" problem is on its way towards being solved. I think I'll keep trying to solve the issue with the shortest paths "missing" some edges. I think I can improve what's there, but it's going to take a little time. Perhaps we can get a list of all of the nodes that were on the shortest paths, and select all edges connecting those nodes too. I think that works out logically.

I've never actually contributed to a github project so I'm not super confident about how the workflow works.

@Robinlovelace
Copy link
Member

Robinlovelace commented Jun 8, 2018 •
edited

@sandrews-CTPS You've already contributed by raising a decent reproducible issue. Congrats 👍 in terms of more contributing you can clone a version of it, play around, e.g. by changing R/SpatialLinesNetwork.R and then rebuilding, with reference to http://r-pkgs.had.co.nz/

This was referenced Jun 8, 2018
Closed
Closed
@Robinlovelace
Copy link
Member

Closing for now. If you have any further comments/questions on this @stevenpandrews please let us know and we'll try to fix them.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
4 participants
@Robinlovelace @mpadge @richardellison @stevenpandrews