UrbanAnalyst / dodgr

Distances on Directed Graphs in R
https://urbananalyst.github.io/dodgr/
127 stars 16 forks source link

Calculate isodistance occupying the whole network without wt_profile #181

Closed demcortillas closed 2 years ago

demcortillas commented 2 years ago

Hello everyone, I hope you are all well

My problem is roughly to build a dodgr_isodistance from a point and convert it into a polygon.

library(dodgr)
library(sf)

osm      <- st_read(dsn="OSM/NET_VIII.shp") # comes from OSM, downloaded with Quick OSM and filtered by "highway".
graph1  <- weight_streetnet(osm, wt_profile = 1)
graph1  <- graph1[graph1$component==1,] 
graph1  <- dodgr_contract_graph(graph1)

vertices        <- dodgr_vertices(graph1)
vertices        <-st_as_sf(vertices, coords = c("x", "y"), crs = 4326, agr = "constant")
p_real <- st_sfc(st_point(c(-73.0511859291954,-36.82369513651161)),crs = 4326)

I             <- st_nearest_feature(p_real,vertices)
I_net      <- vertices[I,]$id
iso_dist  <- dodgr_isodists(graph1, from = I_net , dlim = 10000, heap = "BHeap")

fmat <- matrix (1, nrow = 1, ncol = dim(iso_dist)[1])
f    <- dodgr_flows_aggregate (graph1, 
                               from = iso_dist$from, 
                               to   = iso_dist$id,
                               flows = fmat)          %>% 
                               merge_directed_graph() %>%
                               dodgr_to_sf()
plot(f$geometry)

unexpected lines

When I plot the lines, the result only shows the "main routes" or what I would expect to be "directed" but not the streets that are in a normal CBD (many streets in a small area). Something like in the following image:

ideal lines

My question is the following: How can I occupy the entire network without segregating via wt_profile? A function that does apply to what I'm doing is:

library(spNetwork)
calc_isochrones()

But that function takes a long time, so I am discarding that idea. Thank you very much for your cooperation.

mpadge commented 2 years ago

Thanks @demcortillas. If i understand correctly, your question is just how to extract the isodistance profile as a polygon? If you read the docs for ?dodgr_isodists, you'll see that the return value is described as:

https://github.com/ATFutures/dodgr/blob/fa79c6bb559f7811347dce3d94f207d6ef5943f9/R/iso.R#L16-L20

I think that means that return value is actually what you want, and you don't need to do any of the sf stuff in your example above at all. Here's an example of how to directly plot an isodistance contour, starting with your initial code.

library (dodgr)
osm <- readRDS ("<local>/<osm>/<network>")
graph1  <- weight_streetnet(osm, wt_profile = 1)
#> Loading required namespace: geodist
#> Loading required namespace: dplyr
graph1  <- graph1[graph1$component==1,] 
graph1  <- dodgr_contract_graph(graph1)
vertices <- dodgr_vertices(graph1)

You then use the following code to select a single point and map it onto the network:

# vertices        <-st_as_sf(vertices, coords = c("x", "y"), crs = 4326, agr = "constant")
# p_real <- st_sfc(st_point(c(-73.0511859291954,-36.82369513651161)),crs = 4326)
# I             <- st_nearest_feature(p_real,vertices)
# I_net      <- vertices[I,]$id

But that is not needed as dodgr maps coordinates to nearest network point anyway, so instead you can just submit p_real as your single from value. For the above network I'll fake a value as:

set.seed (1)
p_real <- c (x = min (vertices$x) + runif (1) * diff (range (vertices$x)),
                   y = min (vertices$y) + runif (1) * diff (range (vertices$y)))

The call to isodists is then:

# iso_dist  <- dodgr_isodists(graph1, from = I_net , dlim = 10000, heap = "BHeap")
iso_dist <- dodgr_isodists (graph1, from = p, dlim = 10000)

That result can be used to plot a contour directly, for example using these lines:

plot (iso_dist$x, iso_dist$y, type = "l")
points (p [1], p [2], pch = 19, cex = 2)


But it seems from the image you included above that what you're actually trying to do with the following lines is work out every possible path within that isodistance contour, right?

fmat <- matrix (1, nrow = 1, ncol = dim(iso_dist)[1])
f    <- dodgr_flows_aggregate (graph1, 
                               from = iso_dist$from, 
                               to   = iso_dist$id,
                               flows = fmat)          %>% 
    merge_directed_graph() %>%
    dodgr_to_sf()
plot(f$geometry)

The problem is that isodistance contours in dodgr, in indeed as a general principle of what isodistances actually are, are the collection of points which are maximally far from the defined from point yet just within the specified isodistance. So each point is a single network point which is either (1) terminal, or (2) the final point at d < dlim that connects with a subsequent point which is > dlim. In finding that set of "iso-vertices", the algorithms collects candidate points, for example dA. If that point connects to a subsequent point, dB, with dist(from, dB) > dist(from, dA), then dA is dropped from the candidate list and replaced by dB. The result is that you end up with far fewer of those points than the total number of points actually traversed in searching for them. So you can not use an isodistance algorithm to generate the kind of full network coverage it seems like you're wanting.

But what (i think?) you seem to want nevetheless seems reasonable, so i'll have a go at generating it.

mpadge commented 2 years ago

@demcortillas This code aggregates flows along all paths to all points out to the specified isodistance. The easiest way to do this is to make an initial dodgr_distances() call, then select all points within desired limit and aggregate flows to all of those. I've been meaning to look over this code for a while now, and in doing so found a couple of minor issues. So make sure you install this dev version from GitHub and see at least the version number shown below.

library (dodgr)
packageVersion ("dodgr")
#> [1] '0.2.14.38'
osm <- readRDS ("<local>/<osm>/<network>")
graph1  <- weight_streetnet(osm, wt_profile = 1)
#> Loading required namespace: geodist
#> Loading required namespace: dplyr
graph1  <- graph1[graph1$component==1,] 
#graph1  <- dodgr_contract_graph(graph1)
vertices <- dodgr_vertices(graph1)
p <- c (x = mean (vertices$x),
        y = mean (vertices$y))

d <- dodgr_distances (graph1, from = p) [1, ]
dim (vertices); length (d) # same
#> [1] 196788      5
#> [1] 196788
dlim <- 1000
to <- vertices [which (as.vector (d) < dlim), ]
fmat <- matrix (1., ncol = nrow (to))

f <- dodgr_flows_aggregate (graph1, from = p, to = to$id, flows = fmat, contract = TRUE)
#> xy has no named columns; assuming order is x then y
# or, if you read data via st_read which strips away the OSM id values, you'll need this:
# f <- dodgr_flows_aggregate (graph1, from = p, to = to [, c ("x", "y")], flows = 1, contract = TRUE)

# plotting:
library (sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 8.2.0; sf_use_s2() is TRUE
f_sf <- dodgr_to_sf (f)
#> old-style crs object detected; please recreate object with a recent sf::st_crs()
# reduce to easy plotting
f_sf <- f_sf [f_sf$flow > 0, ]
# flows are generally distributed on log-scales:
f_sf$flow <- log10 (f_sf$flow)
plot (f_sf [, c ("flow", "geometry")])

Created on 2022-06-22 by the reprex package (v2.0.1)

Of course if you just want the network without the flow aggregation, then you can just use the above code to get the isochrone, convert that to an sf_polygon() object, use dodgr_to_sf() to convert your OSM network to sf form, and then use st_within() to trip the network to within the polygon.

demcortillas commented 2 years ago

Hello! thank you very much for your reply! i tried the instruction you wrote, but when i ran the following: f <- dodgr_flows_aggregate (graph1, from = p, to = to [, c ("x", "y")], flows = 1, contract = TRUE) The answer was:

xy has no named columns; assuming order is x then y Error in graph_full[[n]][indx_to_full] <- graph[[n]][indx_to_contr] : NAs are not allowed in subscripted assignments

I should mention that when reading the spatial object, I did so with one osm file and one not. Both gave the same response. Any idea? Thank you very much for your answers

mpadge commented 2 years ago

Yeah, sorry, that was a minor bug that you can either avoid by explicitly passing contract = FALSE to dodgr_flows_aggregate(), or just install latest version (>= 0.2.14.46) and it should work regardless.

demcortillas commented 2 years ago

Thank you very much, I installed the latest version from: remotes::install_github("ATFutures/dodgr") but the version you download is: '0.2.14.43'. Forgive my naivety, but how can I upgrade to version '0.2.14.46'?

mpadge commented 2 years ago

Sorry @demcortillas , i neglected to push my local changes. Should be okay now.