JeremyGelb / spNetwork

An R package to perform spatial analysis on networks.
GNU General Public License v2.0
34 stars 2 forks source link

Error: Error in coordsLines[[nearest_line_index[x]]] : attempt to select less than one element in get1index <real> #20

Closed angle-zhang closed 4 months ago

angle-zhang commented 5 months ago

I'm getting this error when trying to implement the code from this tutorial with my own data from OSMdata. Traceback is pasted below

Error in coordsLines[[nearest_line_index[x]]] : 
  attempt to select less than one element in get1index <real>
 traceback()
7: nrow(coordsLine)
6: vapply(2:nrow(coordsLine), function(x) {
       nearestPointOnSegment(coordsLine[(x - 1):x, ], coordsPoint)
   }, FUN.VALUE = c(0, 0, 0))
5: nearestPointOnLine(coordsLines[[nearest_line_index[x]]], coordsPoints[x, 
       ])
4: FUN(X[[i]], ...)
3: vapply(1:nrow(points), function(x) {
       return(nearestPointOnLine(coordsLines[[nearest_line_index[x]]], 
           coordsPoints[x, ]))
   }, FUN.VALUE = c(0, 0))
2: snapPointsToLines2(start_points, lines, "OID")
1: calc_isochrones(lines = pgh_network, start_points = df_center, 
       dists = c(500), weight = "length")

Code here:

pgh_bb <- getbb("Pittsburgh",  display_name_contains = "United States")
pgh <- opq(bbox =  pgh_bb) %>% 
  add_osm_feature(key = 'highway') %>% 
  osmdata_sf() %>% 
  osm_poly2line()

pgh_network <- pgh$osm_lines %>% 
  select(highway)

pgh_network$length <- as.numeric(st_length(pgh_network))
  centroid <- st_centroid(pgh_blocks)
  center_coords <- st_coordinates(centroid)[1:2,] #do for two right now
  df_center <- center_coords %>%
              data.frame() %>%
              mutate("OID" = row_number()) %>%
              st_as_sf(coords = c("X","Y"), crs = proj_crs)

  iso_results <- calc_isochrones(lines = pgh_network,
                                 start_points = df_center,
                                 dists = c(500),
                                 weight = "length"
  )

And here is a diagram of the study area image

JeremyGelb commented 5 months ago

Hello !

I can have a look at this bug, but I need the coordinates of the points you are using. The variable pgh_blocks is not defined in your code above.

angle-zhang commented 5 months ago

Hi Jeremy, thank you I appreciate the help. These are the two coordinates that I'm extracting from PGH_blocks: X Y [1,] 1374390 419121.5 [2,] 1357437 426328.9

JeremyGelb commented 5 months ago

The error is caused by your road network. It comes from OSM data and has the CRS 4326 (lon/lat). spNetwork expects projected data. I can see that your starting points are projected with the crs 2272. When I set everything in the same crs it works well :

library(spNetwork)
library(osmdata)
library(dplyr)
library(sf)

pgh_bb <- getbb("Pittsburgh",  display_name_contains = "United States")

pgh <- opq(bbox =  pgh_bb) %>%
  add_osm_feature(key = 'highway') %>%
  osmdata_sf() %>%
  osm_poly2line()

pgh_network <- pgh$osm_lines %>%
  select(highway) %>%
  st_transform(2272)

pgh_network$length <- as.numeric(st_length(pgh_network))

center_coords <- rbind(c(1374390, 419121.5),
                       c(1357437, 426328.9)
                       )

df_center <- center_coords %>%
  data.frame() %>%
  mutate("OID" = row_number()) %>%
  st_as_sf(coords = c("X1","X2"), crs = 2272)

iso_results <- calc_isochrones(lines = pgh_network,
                               start_points = df_center,
                               dists = c(500),
                               weight = "length"
)

iso_results |>
  mutate(distance = as.factor(distance)) |>
  tm_shape() +
  tm_lines(col = "distance", palette = "viridis") +
  tm_shape(df_center) +
  tm_bubbles(col = "red")

image