Isochrone and isodistance resolution

[pasipasi123]

dplyr_1.1.4 sf_1.0-15 mapview_2.11.2 osmdata_0.2.5 dodgr_0.3.0.015

Thanks for the great package. I'm wondering if it's possible to enhance or increase the resolution of dodgr_isochrone() (and also isodistance) results? I'm analyzing walking and cycling in an urban environment, where short distances are relevant. By resolution I mean a similar measure that's available in osrm::osrmIsochrone() (link) that basically increases the the number of polygon vertices to create a more well defined polygon. IIRC osrm creates a regular reachable point matrix within the catchment area with given resolution/density to create the isochrone polygon. I don't have osrm examples right now because I don't have the router server running (which is why I'm thankful for the dodgr package).

Example code:

library(dodgr)
library(osmdata)
library(mapview)
library(sf)
library(dplyr)

dat <- opq("Espoo Leppävaara") %>%
  add_osm_feature(key = "highway") %>%
  osmdata_sc()
graph <- weight_streetnet(dat, wt_profile = "foot")

x <- dodgr_isochrones(graph, from = "10594643693", tlim = c(150, 300, 600, 900, 1200))

x %>% 
  st_as_sf(coords = c("x", "y"), crs = 4326) %>% 
  group_by(tlim) %>% 
  summarise(do_union = FALSE) %>% 
  st_cast("POLYGON") %>% 
  mapview()

This produces polygons that are fairly crude given the density of the street network (dodgr_to_sf(graph) %>% mapview()).

Creating a one minute walk isochrone does not seem to be possible, because only a few points are created, and the two minute isochrone shows that the some routes needed for reaching the endpoints are not within the polygon area. I'm guessing only added points along the route and also between the routes could help here.

x2 <- dodgr_isochrones(graph, from = "10594643693", tlim = c(120)) 

x2 %>% 
  st_as_sf(coords = c("x", "y"), crs = 4326) %>% 
  group_by(tlim) %>% 
  summarise(do_union = FALSE) %>% 
  st_cast("POLYGON") %>% 
  mapview()

[mpadge]

Good question. The current vertices are those in the network taken just before one step further would exceed the limits, so in that sense are the minimal set of absolute limit vertices in the network. Extra points could be included by expanding these exact limit points to include additional approximate limit points. Let me think about how that might be done and get back to you.

[mpadge]

@pasipasi123 Sorry about delayed response here. I've checked out several ways of potentially filling in the polygons in cases like yours, but none of them end up any better than simply calculating a convex hull around the points. sf provides the sf_convex_hull() function to do that "properly", by taking the earth's ellipsoid into account, but the following code shows that simply calculating a mathematical hull directly from the points is sufficiently accurate:

library (dodgr)
packageVersion ("dodgr")
#> [1] '0.3.0.17'

f <- "junk.Rds"
if (!file.exists (f)) {
    dat <- osmdata::opq("Espoo Leppävaara") |>
        osmdata::add_osm_feature(key = "highway") |>
        osmdata::osmdata_sc(quiet = FALSE)
    saveRDS (dat, f)
} else {
    dat <- readRDS (f)
}
graph <- weight_streetnet(dat, wt_profile = "foot")
#> Loading required namespace: geodist
#> Loading required namespace: dplyr

from <- "10594643693"
v <- dodgr_vertices (graph)
v_from <- v [match (from, v$id), ]
x <- dodgr_isochrones(graph, from = from, tlim = 120)
p <- sfheaders::sf_point (x [, c ("x", "y")]) |>
    sf::st_set_crs (4326) |>
    sf::st_combine ()
p_sf <- sf::st_sf (geometry = sf::st_cast (p, "POLYGON"))
hull_sf <- sf::st_sf (geometry = sf::st_convex_hull (p))
h <- grDevices::chull (x [, c ("x", "y")])
hull <- sfheaders::sf_polygon (x [h, c ("x", "y")]) [, "geometry"]
sf::st_crs (hull) <- 4326
dat <- rbind (p_sf, hull, hull_sf)
plot (dat, col = "transparent", lwd = 1, border = c ("red", "blue", "green"))

sf::st_area (hull_sf) - sf::st_area (hull)
#> -1.251829e-08 [m^2]
sf::st_perimeter (hull_sf) - sf::st_perimeter (hull)
#> -6.434675e-11 [m]

^{Created on 2024-05-17 with reprex v2.1.0}

Any other way to extract "well behaved" polygons using the street network is (1) far from trivial, and (2) unavoidably dependent on algorithm design with the unavoidable consequence of results being harder to interpret, and therefore likely harder to use.

---

My proposal here is to add a convex parameter to the function with a default of FALSE for current behaviour, but if set to TRUE, then the hull is reduced to the convex hull only, like the green line shown above.

[mpadge]

The above commits now enable this:

library (dodgr)
packageVersion ("dodgr")
#> [1] '0.3.0.19'

f <- "junk.Rds"
if (!file.exists (f)) {
    dat <- osmdata::opq("Espoo Leppävaara") |>
        osmdata::add_osm_feature(key = "highway") |>
        osmdata::osmdata_sc(quiet = FALSE)
    saveRDS (dat, f)
} else {
    dat <- readRDS (f)
}
graph <- weight_streetnet(dat, wt_profile = "foot")
#> Loading required namespace: geodist
#> Loading required namespace: dplyr

from <- "10594643693"

v <- dodgr_vertices (graph)
tlim <- 120 * 1:2
x <- dodgr_isochrones(graph, from = from, tlim = tlim)
table (x$tlim)
#> 
#> 120 240 
#>  32  18
xc <- dodgr_isochrones(graph, from = from, tlim = tlim, convex = TRUE)
table (xc$tlim)
#> 
#> 120 240 
#>   8  11

one_sf_poly <- function (x, tlim) {
    res <- sfheaders::sf_point (x [x$tlim == tlim, c ("x", "y")]) |>
        sf::st_set_crs (4326) |>
        sf::st_combine ()
    sf::st_sf (geometry = sf::st_cast (res, "POLYGON"))
}
dat <- rbind (
    one_sf_poly (x, 120), one_sf_poly (x, 240),
    one_sf_poly (xc, 120), one_sf_poly (xc, 240)
)
plot (dat, col = "transparent", lwd = 1, lty = c (2, 2, 1, 1), border = c ("red", "red", "blue", "blue"))

^{Created on 2024-05-17 with reprex v2.1.0}

[mpadge]

Re-opening, because the convex parameter alone still doesn't resolve this issue. The current PR #226 produces the following plot of convex hull isodistance contours for 3 dlim values:

Each successive value should entirely contain all prior ones, yet that clearly doesn't happen here. The whole iso routines might be better re-written to simply use all points (like current dodgr_isoverts() function), and to calculate hulls around those?

[mpadge]

Those commits finish implementation of new concavity parameter which enables arbitrarily concave polygonal hulls to be fitted to iso results:

library(dodgr)
library(osmdata)
library(mapview)
library(sf)
library(dplyr)

f <- "el.Rds"
if (!file.exists (f)) {
    dat <- opq("Espoo Leppävaara") %>%
        add_osm_feature(key = "highway") %>%
        osmdata_sc()
    saveRDS (dat, f)
} else {
    dat <- readRDS (f)
}
graph <- weight_streetnet(dat, wt_profile = "foot")

# Note the new concavity parameter:
x <- dodgr_isochrones(graph, from = "10594643693", tlim = c(150, 300, 600, 900, 1200), concavity = 0.5)

m <- x %>% 
    st_as_sf(coords = c("x", "y"), crs = 4326) %>% 
    group_by(tlim) %>% 
    summarise(do_union = FALSE) %>% 
    st_cast("POLYGON") %>%
    mapview ()

# Lines to plot mapview result as png:
map2img <- function (m) {
    fout = tempfile(fileext = ".png")
    withr::with_envvar (
        c ("CHROMOTE_CHROME" = "/usr/bin/brave"),
        mapshot2(m, file = fout)
    )
    png::readPNG(fout)
}
grid::grid.raster (map2img (m))

x2 <- dodgr_isochrones(graph, from = "10594643693", tlim = c(120), concavity = 0.5) 
m2 <- x2 %>% 
    st_as_sf(coords = c("x", "y"), crs = 4326) %>% 
    group_by(tlim) %>% 
    summarise(do_union = FALSE) %>% 
    st_cast("POLYGON") %>%
    mapview ()
grid::grid.raster (map2img (m2))

^{Created on 2024-05-24 with reprex v2.1.0}