Unable to purrr::map multiple sfc_point to osrmIsochrone function in R

Question

Relatively new to purrr question: I'm walking through the great tutorial available here and have been successful up until it requires generating isochrones for each major trauma centre per the following provided code:

mtc_iso <- major_trauma_centres %>%
  mutate(iso = map(geometry,
                   osrmIsochrone,
                   breaks = seq(0, 150, by = 30),
                   returnclass = "sf")) %>%
  st_drop_geometry() %>%
  unnest(iso) %>%
  st_set_geometry("geometry") %>%
  st_transform(crs = 27700)

The error generated is:

Error in `stopifnot()`:
ℹ In argument: `iso = map(geometry, osrmIsochrone, breaks = seq(0, 150, by = 30), returnclass = "sf")`.
Caused by error in `map()`:
ℹ In index: 1.
Caused by error:
! "loc" should be a vector of coordinates, a data.frame or a matrix of coordinates, an sfc POINT object or an sf POINT object.

I've examined the major_trauma_centres object and its geometry element and they appear to conform to loc's requirements, yet the error remains:

> class(major_trauma_centres)
[1] "sf"         "tbl_df"     "tbl"        "data.frame"
> class(major_trauma_centres$geometry)
[1] "sfc_POINT" "sfc"     

I've also stepped through providing individual rows to the osrmIsochrone function with success, but I'm having no luck in trying to use the map function across the whole list of geometries.

Any assistance would be greatly appreciated!

Answer 1

I'd guess osrmIsochrone has changed since, another indicator is the use of deprecated returnclass argument. When using map() like that, it cycles through every element of geometry and osrmIsochrone() receives a sfg POINT:

> str(major_trauma_centres$geometry[[1]])
 'XY' num [1:2] 0.139 52.174
> class(major_trauma_centres$geometry[[1]])
[1] "XY"    "POINT" "sfg" 

And it will fail during inherits() checks.

While map() can be used, perhaps consider a rowwise aproach with nested tibbles instead. We'll first nest the geometry column to create 1 sf object per row, those can be passed directly to osrmIsochrone() when we have changed grouping to row-wise:

library(dplyr)
library(tidyr)
library(sf)
library(osrm)
library(mapview)

# 2 first features from linked sample:
geojson <- 
r"(
{
  "type": "FeatureCollection",
  "name": "major_trauma_centres",
  "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },
  "features": [
  { "type": "Feature", "properties": { "name": "ADDENBROOKE'S HOSPITAL", "org_id": "RGT01" }, 
    "geometry": { "type": "Point", "coordinates": [ 0.139114480160586, 52.173740681989244 ] } 
  },
  { "type": "Feature", "properties": { "name": "DERRIFORD HOSPITAL", "org_id": "RK950" }, 
    "geometry": { "type": "Point", "coordinates": [ -4.113684475426827, 50.416719948823165 ] } 
  }]
}
)"
major_trauma_centres <- read_sf(geojson)
major_trauma_centres
#> Simple feature collection with 2 features and 2 fields
#> Geometry type: POINT
#> Dimension:     XY
#> Bounding box:  xmin: -4.113684 ymin: 50.41672 xmax: 0.1391145 ymax: 52.17374
#> Geodetic CRS:  WGS 84
#> # A tibble: 2 × 3
#>   name                   org_id             geometry
#>   <chr>                  <chr>           <POINT [°]>
#> 1 ADDENBROOKE'S HOSPITAL RGT01  (0.1391145 52.17374)
#> 2 DERRIFORD HOSPITAL     RK950  (-4.113684 50.41672)

mtc_iso <- major_trauma_centres %>%
  nest(sf = geometry) %>% 
  rowwise() %>% 
  # # A tibble: 2 × 3
  # # Rowwise: 
  #   name                   org_id sf          
  #   <chr>                  <chr>  <list>      
  # 1 ADDENBROOKE'S HOSPITAL RGT01  <sf [1 × 1]>
  # 2 DERRIFORD HOSPITAL     RK950  <sf [1 × 1]> 
  
  mutate(sf = list(osrmIsochrone(sf, breaks = seq(0, 150, by = 30)))) %>% 
  # # A tibble: 2 × 3
  # # Rowwise: 
  #   name                   org_id sf          
  #   <chr>                  <chr>  <list>      
  # 1 ADDENBROOKE'S HOSPITAL RGT01  <sf [5 × 4]>
  # 2 DERRIFORD HOSPITAL     RK950  <sf [5 × 4]>
  
  unnest(sf) %>% 
  # # A tibble: 10 × 6
  #    name                   org_id    id isomin isomax                                                       geometry
  #    <chr>                  <chr>  <int>  <dbl>  <dbl>                                             <MULTIPOLYGON [°]>
  #  1 ADDENBROOKE'S HOSPITAL RGT01      1      0     30 (((-0.1396775 52.29685, -0.1948542 52.23069, -0.1619805 52.11…
  #  2 ADDENBROOKE'S HOSPITAL RGT01      2     30     60 (((0.4178975 52.66022, 0.2320392 52.61955, 0.06204679 52.5708…
  # .. ...
  st_sf() %>% 
  st_transform(crs = 27700)

mtc_iso
#> Simple feature collection with 10 features and 5 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: 122723 ymin: 2911.149 xmax: 665720.2 ymax: 413292.3
#> Projected CRS: OSGB36 / British National Grid
#> # A tibble: 10 × 6
#>    name                   org_id    id isomin isomax                    geometry
#>  * <chr>                  <chr>  <int>  <dbl>  <dbl>          <MULTIPOLYGON [m]>
#>  1 ADDENBROOKE'S HOSPITAL RGT01      1      0     30 (((526960.8 268154.6, 5233…
#>  2 ADDENBROOKE'S HOSPITAL RGT01      2     30     60 (((563626.1 309688.1, 5511…
#>  3 ADDENBROOKE'S HOSPITAL RGT01      3     60     90 (((562689.8 337536.3, 5618…
#>  4 ADDENBROOKE'S HOSPITAL RGT01      4     90    120 (((487023.5 373144.4, 4745…
#>  5 ADDENBROOKE'S HOSPITAL RGT01      5    120    150 (((461692.4 413292.3, 4576…
#>  6 DERRIFORD HOSPITAL     RK950      1      0     30 (((256930.5 74023.6, 24695…
#>  7 DERRIFORD HOSPITAL     RK950      2     30     60 (((257811 106269.7, 254748…
#>  8 DERRIFORD HOSPITAL     RK950      3     60     90 (((323479.7 120954, 310482…
#>  9 DERRIFORD HOSPITAL     RK950      4     90    120 (((322451.3 155866.8, 3108…
#> 10 DERRIFORD HOSPITAL     RK950      5    120    150 (((363210.5 201731, 350329…
mapview(mtc_iso, z = "id")

^{Created on 2023-06-30 with reprex v2.0.2}

For map(), something like this would work:

major_trauma_centres %>%
  mutate(iso = map(geometry, \(pnt) unclass(pnt) %>% osrmIsochrone(breaks = seq(0, 150, by = 30))))

This will pass each POINT feature though unclass() first and osrmIsochrone() will receive a vector of coordinates, one of the acceptable formats for loc argument.