Bug or user error in attempt to generate scenario file?

[Robinlovelace]

I'm hitting this issue when trying to generate a simple scenario. I guess it's a small mistake in my input. However, the error message could certainly be less cryptic:

# Notes for seminar


# Set-up ------------------------------------------------------------------

remotes::install_github("a-b-street/abstr")
#> Using github PAT from envvar GITHUB_PAT
#> Skipping install of 'abstr' from a github remote, the SHA1 (24baec9c) has not changed since last install.
#>   Use `force = TRUE` to force installation
library(tidyverse)
library(tmap)
library(pct)
tmap_mode("view")
#> tmap mode set to interactive viewing
sf::sf_use_s2(FALSE)
#> Spherical geometry (s2) switched off

u = "https://github.com/ITSLeeds/TDS/raw/master/messages/leeds-north.geojson"
study_region = sf::read_sf(u)
pct_region = pct::pct_regions[study_region, ]
#> old-style crs object detected; please recreate object with a recent sf::st_crs()
#> old-style crs object detected; please recreate object with a recent sf::st_crs()
#> old-style crs object detected; please recreate object with a recent sf::st_crs()
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
#> old-style crs object detected; please recreate object with a recent sf::st_crs()
#> old-style crs object detected; please recreate object with a recent sf::st_crs()
zones_all = pct::get_pct_zones(region = pct_region$region_name, geography = "msoa")
zones = zones_all[study_region, ]
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
plot(zones$geometry)
plot(study_region$geometry, add = TRUE, lwd = 5)

od_all = pct::get_od()
#> No region provided. Returning national OD data.
#> Rows: 2402201 Columns: 14
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr  (2): Area of residence, Area of workplace
#> dbl (12): All categories: Method of travel to work, Work mainly at or from h...
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> Rows: 7201 Columns: 6
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): MSOA11CD, MSOA11NM
#> dbl (4): BNGEAST, BNGNORTH, LONGITUDE, LATITUDE
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

od_study_region = od_all %>% 
  filter(geo_code1 %in% zones$geo_code) %>% 
  filter(geo_code2 %in% zones$geo_code) 

nrow(od_study_region) # 254k od pairs
#> [1] 254
sum(od_study_region$all) # 22k
#> [1] 22224

osm_region = osmextract::oe_get(zones, boundary = zones, layer = "multipolygons")
#> although coordinates are longitude/latitude, st_contains assumes that they are planar
#> The input place was matched with West Yorkshire. 
#> The chosen file was already detected in the download directory. Skip downloading.
#> Warning: The boundary is composed by more than one features. Selecting the
#> first.
#> Start with the vectortranslate operations on the input file!
#> 0...10...20...30...40...50...60...70...80...90...100 - done.
#> Finished the vectortranslate operations on the input file!
#> Reading layer `multipolygons' from data source 
#>   `/mnt/57982e2a-2874-4246-a6fe-115c199bc6bd/data/osm/geofabrik_west-yorkshire-latest.gpkg' 
#>   using driver `GPKG'
#> Simple feature collection with 5365 features and 25 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -2.173294 ymin: 53.51973 xmax: -1.198814 ymax: 53.96315
#> Geodetic CRS:  WGS 84
osm_region %>% 
  sample_n(100) %>% 
  qtm()
#> Linking to GEOS 3.10.1, GDAL 3.4.0, PROJ 8.2.0; sf_use_s2() is FALSE

desire_lines = od::od_to_sf(od_study_region, zones)
#> 0 origins with no match in zone ids
#> 0 destinations with no match in zone ids
#>  points not in od data removed.

names(desire_lines)
#>  [1] "geo_code1"     "geo_code2"     "all"           "from_home"    
#>  [5] "light_rail"    "train"         "bus"           "taxi"         
#>  [9] "motorbike"     "car_driver"    "car_passenger" "bicycle"      
#> [13] "foot"          "other"         "geo_name1"     "geo_name2"    
#> [17] "la_1"          "la_2"          "geometry"
names(desire_lines_leeds)
#>  [1] "area_of_residence" "area_of_workplace" "all"              
#>  [4] "home"              "metro"             "train"            
#>  [7] "bus"               "taxi"              "motorcycle"       
#> [10] "drive"             "passenger"         "bicycle"          
#> [13] "on_foot"           "other"             "geometry"

desire_lines = desire_lines %>%
  select(geo_code1, geo_code2, Walk = foot, Bike = bicycle, Drive = car_driver, all)
desire_lines_small = desire_lines %>% 
  filter(geo_code1 != geo_code2) %>% 
  top_n(n = 10, wt = all)
plot(desire_lines_small)

sum(desire_lines_small$all)
#> [1] 7574

# A/B Street scenario -----------------------------------------------------



ablines = abstr::ab_scenario(od = desire_lines_small, zones = zones, 
                             origin_buildings = osm_region,
                             destination_buildings = osm_region,
                             pop_var = "all")
#> Warning in st_centroid.sf(origin_buildings): st_centroid assumes attributes are
#> constant over geometries of x
#> Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
#> of_largest_polygon): st_centroid does not give correct centroids for longitude/
#> latitude data
#> Warning in st_centroid.sf(destination_buildings): st_centroid assumes attributes
#> are constant over geometries of x
#> Warning in st_centroid.sfc(st_geometry(x), of_largest_polygon =
#> of_largest_polygon): st_centroid does not give correct centroids for longitude/
#> latitude data
#> 0 origins with no match in zone ids
#> 0 destinations with no match in zone ids
#>  points not in od data removed.
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
#> although coordinates are longitude/latitude, st_intersects assumes that they are planar
#> Error in odc_new[od[[1]] == i, "ox"] <- sj_df$x[sel_sj_o]: replacement has length zero

^{Created on 2022-03-10 by the reprex package (v2.0.1)}