diff --git a/DESCRIPTION b/DESCRIPTION index f18ec2d..b2ead9d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -30,12 +30,12 @@ RoxygenNote: 7.1.1 Imports: sf, raster, - methods + methods, + pbapply Depends: R (>= 2.10) Suggests: geodist, - pbapply, terra, colorspace, knitr, @@ -43,6 +43,10 @@ Suggests: ceramic, bookdown, covr, - testthat + testthat, + osmextract, + stplanr, + dplyr, + rgdal VignetteBuilder: knitr Config/testthat/edition: 3 diff --git a/NAMESPACE b/NAMESPACE index 3109ab4..c71d8f7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -2,7 +2,6 @@ export(elevation_extract) export(elevations_get) -export(plot_dz) export(plot_slope) export(sequential_dist) export(slope_3d) diff --git a/R/plot_slope.R b/R/plot_slope.R index e3aa117..61801cc 100644 --- a/R/plot_slope.R +++ b/R/plot_slope.R @@ -1,17 +1,17 @@ #' Plot slope data for a 3d linestring with base R graphics #' -#' @param r A linestring with xyz dimensions +#' @inheritParams slope_raster #' @inheritParams plot_dz #' @inheritParams sequential_dist #' #' @export #' @examples #' plot_slope(lisbon_route_3d) -#' r = lisbon_road_segment_3d -#' plot_slope(r) -#' r = lisbon_road_segment_3d -plot_slope = function(r, - lonlat = sf::st_is_longlat(r), +#' routes = lisbon_road_segment_3d +#' plot_slope(routes) +#' routes = lisbon_road_segment_3d +plot_slope = function(routes, + lonlat = sf::st_is_longlat(routes), fill = TRUE, horiz = FALSE, p = ifelse( @@ -27,7 +27,7 @@ plot_slope = function(r, title = "Slope colors (percentage gradient)", s = 3:18, ncol = 4) { - dz = distance_z(r, lonlat = lonlat) + dz = distance_z(routes, lonlat = lonlat) plot_dz(dz$d, dz$z, fill = fill) } #' Plot a digital elevation profile based on xyz data @@ -48,14 +48,14 @@ plot_slope = function(r, #' @param ncol Number of columns in legend #' @param horiz Should the legend be horizontal (`FALSE` by default) #' @param ... Additional parameters to pass to legend -#' @export #' @examples -#' r = lisbon_road_segment_3d -#' m = sf::st_coordinates(r) +#' routes = lisbon_road_segment_3d +#' m = sf::st_coordinates(routes) #' d = cumsum(sequential_dist(m, lonlat = FALSE)) #' d = c(0, d) #' z = m[, 3] -#' plot_dz(d, z) +#' # not exported +#' # plot_dz(d, z) plot_dz = function(d, z, fill = TRUE, @@ -78,7 +78,7 @@ plot_dz = function(d, if (fill) { b = make_breaks(brks) pal = make_pal(p, b) - g = slope_vector(x = d, e = z) + g = slope_vector(x = d, elevations = z) colz = make_colz(g, b, pal) lapply(seq(d)[-(length(d))], function(i) { graphics::polygon( @@ -95,15 +95,9 @@ plot_dz = function(d, ncol = ncol, cex = cex) } } -# g = slope_matrix(m) -# abline(h = 0, lty = 2) -# points(x[-length(x)], gx, col = "red") -# points(x[-length(x)], gxy, col = "blue") -# title("Distance (in x coordinates) elevation profile", -# sub = "Points show calculated gradients of subsequent lines") -distance_z = function(r, lonlat) { - m = sf::st_coordinates(r) +distance_z = function(routes, lonlat) { + m = sf::st_coordinates(routes) d = cumsum(sequential_dist(m, lonlat = lonlat)) d = c(0, d) z = m[, 3] diff --git a/R/slope_get.R b/R/slope_get.R index 3dfdfbc..1f09a21 100644 --- a/R/slope_get.R +++ b/R/slope_get.R @@ -7,44 +7,35 @@ #' @examples #' # time-consuming, and see #' # https://github.com/ITSLeeds/slopes/runs/648128378#step:18:107 -#' r = cyclestreets_route -#' # e = elevations_get(r) +#' routes = cyclestreets_route +#' # e = elevations_get(routes) #' # class(e) #' # e #' # plot(e) -#' # plot(sf::st_geometry(r), add = TRUE) -elevations_get = function(r, file = NULL, ...) { +#' # plot(sf::st_geometry(routes), add = TRUE) +elevations_get = function(routes, file = NULL, ...) { if(requireNamespace("ceramic")) { - mid_ext = sf_mid_ext_lonlat(r) + mid_ext = sf_mid_ext_lonlat(routes) bw = max(c(mid_ext$width, mid_ext$height)) / 1 # buffer width e = ceramic::cc_elevation(loc = mid_ext$midpoint, buffer = bw, ...) } else { message("Install the package ceramic") } # issue: cannot convert CRS currently - # cr = sf::st_crs(r) + # cr = sf::st_crs(routes) # cr$wkt # raster::crs(e) = raster::crs("+init=epsg:3857") # raster::projectRaster(e, "+init=epsg:4326") - # raster::projectRaster(e, sf::st_crs(r)[[2]]) + # raster::projectRaster(e, sf::st_crs(routes)[[2]]) e } -# -# trace_projected = sf::st_transform(trace, 3857) -# plot(e) -# plot(trace_projected$geometry, add = TRUE) -# -# # aim: get max distance from centrepoint -# bb = sf::st_bbox(sf::st_transform(trace, 4326)) -# # max of those 2 and divide by 2 - -sf_mid_ext_lonlat = function(r) { +sf_mid_ext_lonlat = function(routes) { res = list() - if(!sf::st_is_longlat(r)) { - r = sf::st_transform(r, 4326) + if(!sf::st_is_longlat(routes)) { + routes = sf::st_transform(routes, 4326) } - bb = sf::st_bbox(r) + bb = sf::st_bbox(routes) res$midpoint = c(mean(c(bb[1], bb[3])), mean(c(bb[2], bb[4]))) res$width = geodist::geodist(c(x = bb[1], y = bb[2]), c(x = bb[3], y = bb[2])) res$height = geodist::geodist( diff --git a/R/slopes.R b/R/slopes.R index 107aebb..87cad1a 100644 --- a/R/slopes.R +++ b/R/slopes.R @@ -20,39 +20,39 @@ #' #' @param x Vector of locations #' @param d Vector of distances between points -#' @param e Elevations in same units as x (assumed to be metres) +#' @param elevations Elevations in same units as x (assumed to be metres) #' @export #' @examples #' x = c(0, 2, 3, 4, 5, 9) -#' e = c(1, 2, 2, 4, 3, 1) / 10 -#' slope_vector(x, e) +#' elevations = c(1, 2, 2, 4, 3, 1) / 10 +#' slope_vector(x, elevations) #' m = sf::st_coordinates(lisbon_road_segment) #' d = sequential_dist(m, lonlat = FALSE) -#' e = elevation_extract(m, dem_lisbon_raster) -#' slope_distance(d, e) -#' slope_distance_mean(d, e) -#' slope_distance_weighted(d, e) -slope_vector = function(x, e) { +#' elevations = elevation_extract(m, dem_lisbon_raster) +#' slope_distance(d, elevations) +#' slope_distance_mean(d, elevations) +#' slope_distance_weighted(d, elevations) +slope_vector = function(x, elevations) { d = diff(x) - e_change = diff(e) + e_change = diff(elevations) e_change / d } #' @rdname slope_vector #' @export -slope_distance = function(d, e) { - e_change = diff(e) +slope_distance = function(d, elevations) { + e_change = diff(elevations) e_change / d } #' @rdname slope_vector #' @export -slope_distance_mean = function(d, e) { - e_change = diff(e) +slope_distance_mean = function(d, elevations) { + e_change = diff(elevations) mean(abs(e_change) / d) } #' @rdname slope_vector #' @export -slope_distance_weighted = function(d, e) { - e_change = diff(e) +slope_distance_weighted = function(d, elevations) { + e_change = diff(elevations) stats::weighted.mean(abs(e_change) / d, d) } #' Calculate the gradient of line segments from a 3D matrix of coordinates @@ -74,23 +74,27 @@ slope_distance_weighted = function(d, e) { #' points(x[-length(x)], gxy, col = "blue") #' title("Distance (in x coordinates) elevation profile", #' sub = "Points show calculated gradients of subsequent lines") -slope_matrix = function(m, e = m[, 3], lonlat = TRUE) { +slope_matrix = function(m, elevations = m[, 3], lonlat = TRUE) { d = sequential_dist(m, lonlat = lonlat) - e_change = diff(e) + e_change = diff(elevations) g = e_change / d g } #' @rdname slope_matrix #' @export -slope_matrix_weighted = function(m, e = m[, 3], lonlat = TRUE) { - g1 = slope_matrix(m, e = e, lonlat = lonlat) +slope_matrix_weighted = function(m, elevations = m[, 3], lonlat = TRUE) { + g1 = slope_matrix(m, elevations = elevations, lonlat = lonlat) d = sequential_dist(m = m, lonlat = lonlat) stats::weighted.mean(abs(g1), d, na.rm = TRUE) } #' Calculate the sequential distances between sequential coordinate pairs #' -#' @param lonlat Are the coordinates in lon/lat order? `TRUE` by default +#' @param lonlat Are the coordinates in lon/lat (geographic) coordinates? +#' Set this to `FALSE` if you have projected data, e.g. with coordinates +#' representing distance in meters, not degrees. Lonlat coodinates are assumed +#' (`lonlat = TRUE`) is the default. In `slope_raster()` and `plot_slope()` +#' it is taken from the CRS of the routes (`sf::st_is_longlat(routes)`). #' @inheritParams slope_matrix #' @export #' @examples @@ -111,11 +115,7 @@ sequential_dist = function(m, lonlat = TRUE) { } slope_matrices = function(m_xyz_split, fun = slope_matrix_weighted, ...) { - if (requireNamespace("pbapply", quietly = TRUE)) { - slope_list = pbapply::pblapply(m_xyz_split, fun, ...) - } else { - slope_list = lapply(m_xyz_split, fun, ...) - } + slope_list = pbapply::pblapply(m_xyz_split, fun, ...) unlist(slope_list) } #' Calculate the gradient of line segments from a raster dataset @@ -128,16 +128,16 @@ slope_matrices = function(m_xyz_split, fun = slope_matrix_weighted, ...) { #' if the network is first split-up, e.g. using the function #' `stplanr::rnet_breakup_vertices()` from the #' [`stplanr`](https://docs.ropensci.org/stplanr/reference/) package. -#' **Note:** The `r` object must have a geometry type of `LINESTRING`. +#' **Note:** The `routes` object must have a geometry type of `LINESTRING`. #' The `sf::st_cast()` function can convert from `MULTILINESTRING` (and other) #' geometries to `LINESTRING`s as follows: -#' `r_linestring = sf::st_cast(r, "LINESTRING")`. +#' `r_linestring = sf::st_cast(routes, "LINESTRING")`. #' #' @inheritParams sequential_dist #' @inheritParams elevation_extract -#' @param r Routes, the gradients of which are to be calculated. +#' @param routes Routes, the gradients of which are to be calculated. #' The object must be of class `sf` with `LINESTRING` geometries. -#' @param e Raster overlapping with `r` and values representing elevations +#' @param dem Raster overlapping with `routes` and values representing elevations #' @param method The method of estimating elevation at points, #' passed to the `extract` function for extracting values from raster #' datasets. Default: `"bilinear"`. @@ -146,30 +146,30 @@ slope_matrices = function(m_xyz_split, fun = slope_matrix_weighted, ...) { #' @importFrom methods is #' @export #' @examples -#' r = lisbon_road_segments[1:3, ] -#' e = dem_lisbon_raster -#' (s = slope_raster(r, e)) -#' cor(r$Avg_Slope, s) +#' routes = lisbon_road_segments[1:3, ] +#' dem = dem_lisbon_raster +#' (s = slope_raster(routes, dem)) +#' cor(routes$Avg_Slope, s) slope_raster = function( - r, - e = NULL, - lonlat = sf::st_is_longlat(r), + routes, + dem = NULL, + lonlat = sf::st_is_longlat(routes), method = "bilinear", fun = slope_matrix_weighted, - terra = has_terra() && methods::is(e, "SpatRaster") + terra = has_terra() && methods::is(dem, "SpatRaster") ) { - stop_is_not_linestring(r) - r = sf::st_geometry(r) + stop_is_not_linestring(routes) + routes = sf::st_geometry(routes) # todo: split out this bit into slope_xyz function - m = sf::st_coordinates(r) + m = sf::st_coordinates(routes) # colnames(m) - z = elevation_extract(m, e, method = method, terra = terra) + z = elevation_extract(m, dem, method = method, terra = terra) m_xyz_df = data.frame(x = m[, "X"], y = m[, "Y"], z = z, L1 = m[, "L1"]) res = slope_xyz(m_xyz_df, fun = fun, lonlat = lonlat) res } -#' Extract slopes from xyz dataframe or sf objects +#' Extract slopes from xyz data frame or sf objects #' #' @param r_xyz An sf object with x, y, z dimensions #' @inheritParams slope_raster @@ -177,23 +177,19 @@ slope_raster = function( #' @export #' @examples #' r_xyz = lisbon_road_segment_3d -#' slope_xyz(r_xyz) -slope_xyz = function(r_xyz, fun = slope_matrix_weighted, lonlat = NULL) { - # browser() +#' slope_xyz(r_xyz, lonlat = FALSE) +slope_xyz = function(r_xyz, fun = slope_matrix_weighted, lonlat = TRUE) { if(inherits(r_xyz, "sf") | inherits(r_xyz, "sfc")) { - if(is.null(lonlat)) { - lonlat = sf::st_is_longlat(r_xyz) - } + lonlat = sf::st_is_longlat(r_xyz) r_xyz = as.data.frame(sf::st_coordinates(r_xyz)) } - if(is.null(lonlat)) { - lonlat = FALSE - } if("L1" %in% colnames(r_xyz)) { m_xyz_split = split(x = r_xyz, f = r_xyz[, "L1"]) res = slope_matrices(m_xyz_split, lonlat = lonlat, fun = fun) } else { - # todo: add content here + # todo: add content here if data frame was generated by sfheaders + # or another package that does not call id colums 'L1' by default + # (low priority) } res } @@ -202,15 +198,15 @@ slope_xyz = function(r_xyz, fun = slope_matrix_weighted, lonlat = NULL) { #' #' @param terra Should the `terra` package be used? #' `TRUE` by default if the package is installed *and* -#' if `e` is of class `SpatRast` +#' if `dem` is of class `SpatRast` #' @inheritParams slope_raster #' @inheritParams slope_matrix #' @export #' @examples #' m = sf::st_coordinates(lisbon_road_segments[1, ]) -#' e = dem_lisbon_raster -#' elevation_extract(m, e) -#' elevation_extract(m, e, method = "simple") +#' dem = dem_lisbon_raster +#' elevation_extract(m, dem) +#' elevation_extract(m, dem, method = "simple") #' # uncomment the following lines to test with terra (experimental) #' # u = paste0("https://github.com/ITSLeeds/slopes/", #' # "releases/download/0.0.0/dem_lisbon.tif" ) @@ -218,14 +214,14 @@ slope_xyz = function(r_xyz, fun = slope_matrix_weighted, lonlat = NULL) { #' # et = terra::rast("dem_lisbon.tif") #' # elevation_extract(m, et) elevation_extract = function(m, - e, + dem, method = "bilinear", - terra = has_terra() && methods::is(e, "SpatRaster") + terra = has_terra() && methods::is(dem, "SpatRaster") ) { if(terra) { - res = terra::extract(e, m[, 1:2], method = method)[[1]] + res = terra::extract(dem, m[, 1:2], method = method)[[1]] } else { - res = raster::extract(e, m[, 1:2], method = method) + res = raster::extract(dem, m[, 1:2], method = method) } # unlist(res) res @@ -236,65 +232,52 @@ elevation_extract = function(m, #' @inheritParams elevation_extract #' @export #' @examples -#' r = lisbon_road_segments[204, ] -#' e = dem_lisbon_raster -#' (r3d = slope_3d(r, e)) -#' sf::st_z_range(r) +#' routes = lisbon_road_segments[204, ] +#' dem = dem_lisbon_raster +#' (r3d = slope_3d(routes, dem)) +#' sf::st_z_range(routes) #' sf::st_z_range(r3d) #' plot(sf::st_coordinates(r3d)[, 3]) #' plot_slope(r3d) #' # r3d_get = slope_3d(cyclestreets_route) #' # plot_slope(r3d_get) slope_3d = function( - r, - e = NULL, + routes, + dem = NULL, method = "bilinear", - terra = has_terra() && methods::is(e, "SpatRaster") + terra = has_terra() && methods::is(dem, "SpatRaster") ) { - # if("geom" %in% names(r)) { - # rgeom = r$geom - # } else if("geometry" %in% names(r)) { - # rgeom = r$geometry - # } else { - # rgeom = sf::st_geometry(r) - # } - if(is.null(e)) { - e = elevations_get(r) - r_original = r # create copy to deal with projection issues - r = sf::st_transform(r, raster::crs(e)) - suppressWarnings({sf::st_crs(r) = sf::st_crs(r_original)}) - # plot(e) - # plot(r$geometry, add = TRUE) - m = sf::st_coordinates(r) + if(is.null(dem)) { + dem = elevations_get(routes) + r_original = routes # create copy to deal with projection issues + routes = sf::st_transform(routes, raster::crs(dem)) + suppressWarnings({sf::st_crs(routes) = sf::st_crs(r_original)}) + m = sf::st_coordinates(routes) mo = sf::st_coordinates(r_original) - z = as.numeric(elevation_extract(m, e, method = method, terra = terra)) + z = as.numeric(elevation_extract(m, dem, method = method, terra = terra)) m_xyz = cbind(mo[, 1:2], z) } else { - m = sf::st_coordinates(r) - z = as.numeric(elevation_extract(m, e, method = method, terra = terra)) + m = sf::st_coordinates(routes) + z = as.numeric(elevation_extract(m, dem, method = method, terra = terra)) m_xyz = cbind(m[, 1:2], z) } - n = nrow(r) + n = nrow(routes) if(n == 1) { - # currently only works for 1 line, to be generalised - rgeom3d_line = sf::st_linestring(m_xyz) - rgeom3d_sfc = sf::st_sfc(rgeom3d_line, crs = sf::st_crs(r)) - # message("Original geometry: ", ncol(rgeom[[1]])) - sf::st_geometry(r) = rgeom3d_sfc - # message("New geometry: ", ncol(r$geom[[1]])) + rgeom3d_sfc = sf::st_sfc(rgeom3d_line, crs = sf::st_crs(routes)) + sf::st_geometry(routes) = rgeom3d_sfc } else { linestrings = lapply(seq(n), function(i){ rgeom3d_line = sf::st_linestring(m_xyz[m[, 3] == i, ]) }) - rgeom3d_sfc = sf::st_sfc(linestrings, crs = sf::st_crs(r)) - sf::st_geometry(r) = rgeom3d_sfc + rgeom3d_sfc = sf::st_sfc(linestrings, crs = sf::st_crs(routes)) + sf::st_geometry(routes) = rgeom3d_sfc } - r + routes } -# terra = has_terra() && methods::is(e, "SpatRaster") -# terra + +# unexported function to check if the user has terra package installed has_terra = function() { res = requireNamespace("terra", quietly = TRUE) # print(res) @@ -310,9 +293,3 @@ stop_is_not_linestring = function(x) { "Only works with LINESTRINGs. Try converting with sf::st_cast() first." ) } -# # test: -# x = sf::st_cast(lisbon_road_segment, "MULTILINESTRING") -# is_linestring(lisbon_road_segment) -# is_linestring(x) -# stop_is_not_linestring(x) -# stop_is_not_linestring(lisbon_road_segment) diff --git a/README.Rmd b/README.Rmd index c16bcde..79aee5f 100644 --- a/README.Rmd +++ b/README.Rmd @@ -70,7 +70,7 @@ plot(sf::st_geometry(lisbon_road_segments), add = TRUE) Calculate the average gradient of each road segment as follows: ```{r} -lisbon_road_segments$slope = slope_raster(lisbon_road_segments, e = dem_lisbon_raster) +lisbon_road_segments$slope = slope_raster(lisbon_road_segments, dem = dem_lisbon_raster) summary(lisbon_road_segments$slope) ``` @@ -140,67 +140,9 @@ plot_slope(lisbon_route_3d) If you do not have a raster dataset representing elevations, you can automatically download them as follows. ```{r, message=FALSE, warning=FALSE} -lisbon_route_3d_auto = slope_3d(r = lisbon_route) +lisbon_route_3d_auto = slope_3d(lisbon_route) plot_slope(lisbon_route_3d_auto) ``` -# Performance -For this benchmark we will download the following small (< 100 kB) `.tif` file: - -```{r} -u = "https://github.com/ITSLeeds/slopes/releases/download/0.0.0/dem_lisbon.tif" -f = basename(u) -if(!file.exists(f)) { - download.file(u, f) -} -``` - -A benchmark can reveal how many route gradients can be calculated per second: - -```{r, results='hide'} -e = dem_lisbon_raster -r = lisbon_road_segments -et = terra::rast(f) -res = bench::mark(check = FALSE, - slope_raster = slope_raster(r, e), - slope_terra = slope_raster(r, et) -) -``` - -```{r} -res -``` - - -That is approximately - -```{r} -round(res$`itr/sec` * nrow(r)) -``` - -routes per second using the `raster` and `terra` (the default if installed, using `RasterLayer` and native `SpatRaster` objects) packages to extract elevation estimates from the raster datasets, respectively. - -The message: use the `terra` package to read-in DEM data for slope extraction if speed is important. - -To go faster, you can chose the `simple` method to gain some speed at the expense of accuracy: - -```{r, results='hide'} -e = dem_lisbon_raster -r = lisbon_road_segments -res = bench::mark(check = FALSE, - bilinear1 = slope_raster(r, e), - bilinear2 = slope_raster(r, et), - simple1 = slope_raster(r, e, method = "simple"), - simple2 = slope_raster(r, et, method = "simple") -) -``` - -```{r} -res -``` - -```{r, include=FALSE} -file.remove("dem_lisbon.tif") -``` diff --git a/README.md b/README.md index b40f4d0..a5a47c5 100644 --- a/README.md +++ b/README.md @@ -23,7 +23,6 @@ geometries and raster digital elevation model (DEM) datasets. - Install the development version from [GitHub](https://github.com/) with: ``` r @@ -74,7 +73,7 @@ plot(sf::st_geometry(lisbon_road_segments), add = TRUE) Calculate the average gradient of each road segment as follows: ``` r -lisbon_road_segments$slope = slope_raster(lisbon_road_segments, e = dem_lisbon_raster) +lisbon_road_segments$slope = slope_raster(lisbon_road_segments, dem = dem_lisbon_raster) summary(lisbon_road_segments$slope) #> Min. 1st Qu. Median Mean 3rd Qu. Max. #> 0.00000 0.01246 0.03534 0.05462 0.08251 0.27583 @@ -144,84 +143,10 @@ If you do not have a raster dataset representing elevations, you can automatically download them as follows. ``` r -lisbon_route_3d_auto = slope_3d(r = lisbon_route) +lisbon_route_3d_auto = slope_3d(lisbon_route) #> Preparing to download: 12 tiles at zoom = 15 from #> https://api.mapbox.com/v4/mapbox.terrain-rgb/ plot_slope(lisbon_route_3d_auto) ``` - -# Performance - -For this benchmark we will download the following small (< 100 kB) -`.tif` file: - -``` r -u = "https://github.com/ITSLeeds/slopes/releases/download/0.0.0/dem_lisbon.tif" -f = basename(u) -if(!file.exists(f)) { - download.file(u, f) -} -``` - -A benchmark can reveal how many route gradients can be calculated per -second: - -``` r -e = dem_lisbon_raster -r = lisbon_road_segments -et = terra::rast(f) -res = bench::mark(check = FALSE, - slope_raster = slope_raster(r, e), - slope_terra = slope_raster(r, et) -) -``` - -``` r -res -#> # A tibble: 2 x 6 -#> expression min median `itr/sec` mem_alloc `gc/sec` -#> -#> 1 slope_raster 49.4ms 50.4ms 19.9 5.72MB 4.97 -#> 2 slope_terra 51.7ms 52.6ms 18.9 2.17MB 5.41 -``` - -That is approximately - -``` r -round(res$`itr/sec` * nrow(r)) -#> [1] 5391 5127 -``` - -routes per second using the `raster` and `terra` (the default if -installed, using `RasterLayer` and native `SpatRaster` objects) packages -to extract elevation estimates from the raster datasets, respectively. - -The message: use the `terra` package to read-in DEM data for slope -extraction if speed is important. - -To go faster, you can chose the `simple` method to gain some speed at -the expense of accuracy: - -``` r -e = dem_lisbon_raster -r = lisbon_road_segments -res = bench::mark(check = FALSE, - bilinear1 = slope_raster(r, e), - bilinear2 = slope_raster(r, et), - simple1 = slope_raster(r, e, method = "simple"), - simple2 = slope_raster(r, et, method = "simple") -) -``` - -``` r -res -#> # A tibble: 4 x 6 -#> expression min median `itr/sec` mem_alloc `gc/sec` -#> -#> 1 bilinear1 49.7ms 50.7ms 19.8 5.72MB 8.47 -#> 2 bilinear2 52.1ms 52.6ms 18.9 2.12MB 4.73 -#> 3 simple1 42.9ms 43.6ms 22.8 2.05MB 4.56 -#> 4 simple2 50.2ms 51.2ms 19.4 2.12MB 4.84 -``` diff --git a/data-raw/slides.Rmd b/data-raw/slides.Rmd index 2d91e49..5b9dadd 100644 --- a/data-raw/slides.Rmd +++ b/data-raw/slides.Rmd @@ -61,7 +61,9 @@ background-size: 100%
# Structure: -## Uses +## Why slopes? + +-- -- @@ -69,6 +71,32 @@ background-size: 100% -- +--- + + +![](slope-edinburgh-bridge.png) + +??? + +- RL + +- My research into tools for prioritising cycling investment + +- UK not very hilly but models fail slightly in hilly areas + + +--- + +... + + +??? + +- Rosa motivation + +--- + + ## Some results diff --git a/data-raw/slope-validation.R b/data-raw/slope-validation.R index f6b9c2c..cdac482 100644 --- a/data-raw/slope-validation.R +++ b/data-raw/slope-validation.R @@ -33,7 +33,7 @@ usethis::use_data(magnolia_xy, overwrite = TRUE) sf::st_crs(magnolia_xy) # slope_raster(magnolia_xy) # fails -magnolia_xyz = slope_3d(r = magnolia_xy) +magnolia_xyz = slope_3d(magnolia_xy) magnolia_xyz$slopes = round(slope_xyz(magnolia_xyz)*100) summary(magnolia_xyz$SLOPE_PCT) summary(magnolia_xyz$slopes) diff --git a/man/elevation_extract.Rd b/man/elevation_extract.Rd index 41d6726..41f805b 100644 --- a/man/elevation_extract.Rd +++ b/man/elevation_extract.Rd @@ -6,15 +6,15 @@ \usage{ elevation_extract( m, - e, + dem, method = "bilinear", - terra = has_terra() && methods::is(e, "SpatRaster") + terra = has_terra() && methods::is(dem, "SpatRaster") ) } \arguments{ \item{m}{Matrix containing coordinates and elevations} -\item{e}{Raster overlapping with \code{r} and values representing elevations} +\item{dem}{Raster overlapping with \code{routes} and values representing elevations} \item{method}{The method of estimating elevation at points, passed to the \code{extract} function for extracting values from raster @@ -22,16 +22,16 @@ datasets. Default: \code{"bilinear"}.} \item{terra}{Should the \code{terra} package be used? \code{TRUE} by default if the package is installed \emph{and} -if \code{e} is of class \code{SpatRast}} +if \code{dem} is of class \code{SpatRast}} } \description{ Extract elevations from coordinates } \examples{ m = sf::st_coordinates(lisbon_road_segments[1, ]) -e = dem_lisbon_raster -elevation_extract(m, e) -elevation_extract(m, e, method = "simple") +dem = dem_lisbon_raster +elevation_extract(m, dem) +elevation_extract(m, dem, method = "simple") # uncomment the following lines to test with terra (experimental) # u = paste0("https://github.com/ITSLeeds/slopes/", # "releases/download/0.0.0/dem_lisbon.tif" ) diff --git a/man/elevations_get.Rd b/man/elevations_get.Rd index d79d853..b0e7d66 100644 --- a/man/elevations_get.Rd +++ b/man/elevations_get.Rd @@ -4,10 +4,10 @@ \alias{elevations_get} \title{Get elevation data from hosted maptile services} \usage{ -elevations_get(r, file = NULL, ...) +elevations_get(routes, file = NULL, ...) } \arguments{ -\item{r}{Routes, the gradients of which are to be calculated. +\item{routes}{Routes, the gradients of which are to be calculated. The object must be of class \code{sf} with \code{LINESTRING} geometries.} \item{file}{Where to save the resulting data if specified (not implemented)} @@ -20,10 +20,10 @@ Get elevation data from hosted maptile services \examples{ # time-consuming, and see # https://github.com/ITSLeeds/slopes/runs/648128378#step:18:107 -r = cyclestreets_route -# e = elevations_get(r) +routes = cyclestreets_route +# e = elevations_get(routes) # class(e) # e # plot(e) -# plot(sf::st_geometry(r), add = TRUE) +# plot(sf::st_geometry(routes), add = TRUE) } diff --git a/man/figures/README-dem-lisbon-1.png b/man/figures/README-dem-lisbon-1.png index 5e343bc..d0d9be1 100644 Binary files a/man/figures/README-dem-lisbon-1.png and b/man/figures/README-dem-lisbon-1.png differ diff --git a/man/figures/README-plot_slope-1.png b/man/figures/README-plot_slope-1.png index d4aca5d..2bd7381 100644 Binary files a/man/figures/README-plot_slope-1.png and b/man/figures/README-plot_slope-1.png differ diff --git a/man/figures/README-route-1.png b/man/figures/README-route-1.png index 316f210..0e06c9a 100644 Binary files a/man/figures/README-route-1.png and b/man/figures/README-route-1.png differ diff --git a/man/figures/README-slope-vis-1.png b/man/figures/README-slope-vis-1.png index ac4f788..e3c041c 100644 Binary files a/man/figures/README-slope-vis-1.png and b/man/figures/README-slope-vis-1.png differ diff --git a/man/figures/README-unnamed-chunk-10-1.png b/man/figures/README-unnamed-chunk-10-1.png index 6968c06..7ee65f7 100644 Binary files a/man/figures/README-unnamed-chunk-10-1.png and b/man/figures/README-unnamed-chunk-10-1.png differ diff --git a/man/plot_dz.Rd b/man/plot_dz.Rd index 88822fc..5662b43 100644 --- a/man/plot_dz.Rd +++ b/man/plot_dz.Rd @@ -57,10 +57,11 @@ plot_dz( Plot a digital elevation profile based on xyz data } \examples{ -r = lisbon_road_segment_3d -m = sf::st_coordinates(r) +routes = lisbon_road_segment_3d +m = sf::st_coordinates(routes) d = cumsum(sequential_dist(m, lonlat = FALSE)) d = c(0, d) z = m[, 3] -plot_dz(d, z) +# not exported +# plot_dz(d, z) } diff --git a/man/plot_slope.Rd b/man/plot_slope.Rd index e72da12..7536f2b 100644 --- a/man/plot_slope.Rd +++ b/man/plot_slope.Rd @@ -5,8 +5,8 @@ \title{Plot slope data for a 3d linestring with base R graphics} \usage{ plot_slope( - r, - lonlat = sf::st_is_longlat(r), + routes, + lonlat = sf::st_is_longlat(routes), fill = TRUE, horiz = FALSE, p = ifelse(test = requireNamespace("colorspace", quietly = TRUE), @@ -22,9 +22,14 @@ plot_slope( ) } \arguments{ -\item{r}{A linestring with xyz dimensions} +\item{routes}{Routes, the gradients of which are to be calculated. +The object must be of class \code{sf} with \code{LINESTRING} geometries.} -\item{lonlat}{Are the coordinates in lon/lat order? \code{TRUE} by default} +\item{lonlat}{Are the coordinates in lon/lat (geographic) coordinates? +Set this to \code{FALSE} if you have projected data, e.g. with coordinates +representing distance in meters, not degrees. Lonlat coodinates are assumed +(\code{lonlat = TRUE}) is the default. In \code{slope_raster()} and \code{plot_slope()} +it is taken from the CRS of the routes (\code{sf::st_is_longlat(routes)}).} \item{fill}{Should the profile be filled? \code{TRUE} by default} @@ -54,7 +59,7 @@ Plot slope data for a 3d linestring with base R graphics } \examples{ plot_slope(lisbon_route_3d) -r = lisbon_road_segment_3d -plot_slope(r) -r = lisbon_road_segment_3d +routes = lisbon_road_segment_3d +plot_slope(routes) +routes = lisbon_road_segment_3d } diff --git a/man/sequential_dist.Rd b/man/sequential_dist.Rd index 1fea943..0e49623 100644 --- a/man/sequential_dist.Rd +++ b/man/sequential_dist.Rd @@ -9,7 +9,11 @@ sequential_dist(m, lonlat = TRUE) \arguments{ \item{m}{Matrix containing coordinates and elevations} -\item{lonlat}{Are the coordinates in lon/lat order? \code{TRUE} by default} +\item{lonlat}{Are the coordinates in lon/lat (geographic) coordinates? +Set this to \code{FALSE} if you have projected data, e.g. with coordinates +representing distance in meters, not degrees. Lonlat coodinates are assumed +(\code{lonlat = TRUE}) is the default. In \code{slope_raster()} and \code{plot_slope()} +it is taken from the CRS of the routes (\code{sf::st_is_longlat(routes)}).} } \description{ Calculate the sequential distances between sequential coordinate pairs diff --git a/man/slope_3d.Rd b/man/slope_3d.Rd index 8e8ba10..c5431b1 100644 --- a/man/slope_3d.Rd +++ b/man/slope_3d.Rd @@ -5,17 +5,17 @@ \title{Take a linestring and add a third (z) dimension to its coordinates} \usage{ slope_3d( - r, - e = NULL, + routes, + dem = NULL, method = "bilinear", - terra = has_terra() && methods::is(e, "SpatRaster") + terra = has_terra() && methods::is(dem, "SpatRaster") ) } \arguments{ -\item{r}{Routes, the gradients of which are to be calculated. +\item{routes}{Routes, the gradients of which are to be calculated. The object must be of class \code{sf} with \code{LINESTRING} geometries.} -\item{e}{Raster overlapping with \code{r} and values representing elevations} +\item{dem}{Raster overlapping with \code{routes} and values representing elevations} \item{method}{The method of estimating elevation at points, passed to the \code{extract} function for extracting values from raster @@ -23,16 +23,16 @@ datasets. Default: \code{"bilinear"}.} \item{terra}{Should the \code{terra} package be used? \code{TRUE} by default if the package is installed \emph{and} -if \code{e} is of class \code{SpatRast}} +if \code{dem} is of class \code{SpatRast}} } \description{ Take a linestring and add a third (z) dimension to its coordinates } \examples{ -r = lisbon_road_segments[204, ] -e = dem_lisbon_raster -(r3d = slope_3d(r, e)) -sf::st_z_range(r) +routes = lisbon_road_segments[204, ] +dem = dem_lisbon_raster +(r3d = slope_3d(routes, dem)) +sf::st_z_range(routes) sf::st_z_range(r3d) plot(sf::st_coordinates(r3d)[, 3]) plot_slope(r3d) diff --git a/man/slope_matrix.Rd b/man/slope_matrix.Rd index 077b94d..f0dc221 100644 --- a/man/slope_matrix.Rd +++ b/man/slope_matrix.Rd @@ -5,16 +5,20 @@ \alias{slope_matrix_weighted} \title{Calculate the gradient of line segments from a 3D matrix of coordinates} \usage{ -slope_matrix(m, e = m[, 3], lonlat = TRUE) +slope_matrix(m, elevations = m[, 3], lonlat = TRUE) -slope_matrix_weighted(m, e = m[, 3], lonlat = TRUE) +slope_matrix_weighted(m, elevations = m[, 3], lonlat = TRUE) } \arguments{ \item{m}{Matrix containing coordinates and elevations} -\item{e}{Elevations in same units as x (assumed to be metres)} +\item{elevations}{Elevations in same units as x (assumed to be metres)} -\item{lonlat}{Are the coordinates in lon/lat order? \code{TRUE} by default} +\item{lonlat}{Are the coordinates in lon/lat (geographic) coordinates? +Set this to \code{FALSE} if you have projected data, e.g. with coordinates +representing distance in meters, not degrees. Lonlat coodinates are assumed +(\code{lonlat = TRUE}) is the default. In \code{slope_raster()} and \code{plot_slope()} +it is taken from the CRS of the routes (\code{sf::st_is_longlat(routes)}).} } \description{ Calculate the gradient of line segments from a 3D matrix of coordinates diff --git a/man/slope_raster.Rd b/man/slope_raster.Rd index 7f059b0..32a34f8 100644 --- a/man/slope_raster.Rd +++ b/man/slope_raster.Rd @@ -5,21 +5,25 @@ \title{Calculate the gradient of line segments from a raster dataset} \usage{ slope_raster( - r, - e = NULL, - lonlat = sf::st_is_longlat(r), + routes, + dem = NULL, + lonlat = sf::st_is_longlat(routes), method = "bilinear", fun = slope_matrix_weighted, - terra = has_terra() && methods::is(e, "SpatRaster") + terra = has_terra() && methods::is(dem, "SpatRaster") ) } \arguments{ -\item{r}{Routes, the gradients of which are to be calculated. +\item{routes}{Routes, the gradients of which are to be calculated. The object must be of class \code{sf} with \code{LINESTRING} geometries.} -\item{e}{Raster overlapping with \code{r} and values representing elevations} +\item{dem}{Raster overlapping with \code{routes} and values representing elevations} -\item{lonlat}{Are the coordinates in lon/lat order? \code{TRUE} by default} +\item{lonlat}{Are the coordinates in lon/lat (geographic) coordinates? +Set this to \code{FALSE} if you have projected data, e.g. with coordinates +representing distance in meters, not degrees. Lonlat coodinates are assumed +(\code{lonlat = TRUE}) is the default. In \code{slope_raster()} and \code{plot_slope()} +it is taken from the CRS of the routes (\code{sf::st_is_longlat(routes)}).} \item{method}{The method of estimating elevation at points, passed to the \code{extract} function for extracting values from raster @@ -30,7 +34,7 @@ datasets. Default: \code{"bilinear"}.} \item{terra}{Should the \code{terra} package be used? \code{TRUE} by default if the package is installed \emph{and} -if \code{e} is of class \code{SpatRast}} +if \code{dem} is of class \code{SpatRast}} } \description{ This function takes an \code{sf} representing routes over geographical space @@ -42,14 +46,14 @@ If calculating slopes associated with OSM data, the results may be better if the network is first split-up, e.g. using the function \code{stplanr::rnet_breakup_vertices()} from the \href{https://docs.ropensci.org/stplanr/reference/}{\code{stplanr}} package. -\strong{Note:} The \code{r} object must have a geometry type of \code{LINESTRING}. +\strong{Note:} The \code{routes} object must have a geometry type of \code{LINESTRING}. The \code{sf::st_cast()} function can convert from \code{MULTILINESTRING} (and other) geometries to \code{LINESTRING}s as follows: -\code{r_linestring = sf::st_cast(r, "LINESTRING")}. +\code{r_linestring = sf::st_cast(routes, "LINESTRING")}. } \examples{ -r = lisbon_road_segments[1:3, ] -e = dem_lisbon_raster -(s = slope_raster(r, e)) -cor(r$Avg_Slope, s) +routes = lisbon_road_segments[1:3, ] +dem = dem_lisbon_raster +(s = slope_raster(routes, dem)) +cor(routes$Avg_Slope, s) } diff --git a/man/slope_vector.Rd b/man/slope_vector.Rd index 01dd826..819caea 100644 --- a/man/slope_vector.Rd +++ b/man/slope_vector.Rd @@ -7,18 +7,18 @@ \alias{slope_distance_weighted} \title{Calculate the gradient of line segments from distance and elevation vectors} \usage{ -slope_vector(x, e) +slope_vector(x, elevations) -slope_distance(d, e) +slope_distance(d, elevations) -slope_distance_mean(d, e) +slope_distance_mean(d, elevations) -slope_distance_weighted(d, e) +slope_distance_weighted(d, elevations) } \arguments{ \item{x}{Vector of locations} -\item{e}{Elevations in same units as x (assumed to be metres)} +\item{elevations}{Elevations in same units as x (assumed to be metres)} \item{d}{Vector of distances between points} } @@ -41,12 +41,12 @@ on the resulting gradient estimate (see examples). } \examples{ x = c(0, 2, 3, 4, 5, 9) -e = c(1, 2, 2, 4, 3, 1) / 10 -slope_vector(x, e) +elevations = c(1, 2, 2, 4, 3, 1) / 10 +slope_vector(x, elevations) m = sf::st_coordinates(lisbon_road_segment) d = sequential_dist(m, lonlat = FALSE) -e = elevation_extract(m, dem_lisbon_raster) -slope_distance(d, e) -slope_distance_mean(d, e) -slope_distance_weighted(d, e) +elevations = elevation_extract(m, dem_lisbon_raster) +slope_distance(d, elevations) +slope_distance_mean(d, elevations) +slope_distance_weighted(d, elevations) } diff --git a/man/slope_xyz.Rd b/man/slope_xyz.Rd index ae7fce2..1260027 100644 --- a/man/slope_xyz.Rd +++ b/man/slope_xyz.Rd @@ -2,9 +2,9 @@ % Please edit documentation in R/slopes.R \name{slope_xyz} \alias{slope_xyz} -\title{Extract slopes from xyz dataframe or sf objects} +\title{Extract slopes from xyz data frame or sf objects} \usage{ -slope_xyz(r_xyz, fun = slope_matrix_weighted, lonlat = NULL) +slope_xyz(r_xyz, fun = slope_matrix_weighted, lonlat = TRUE) } \arguments{ \item{r_xyz}{An sf object with x, y, z dimensions} @@ -12,12 +12,16 @@ slope_xyz(r_xyz, fun = slope_matrix_weighted, lonlat = NULL) \item{fun}{The slope function to calculate per route, \code{slope_matrix_weighted} by default.} -\item{lonlat}{Are the coordinates in lon/lat order? \code{TRUE} by default} +\item{lonlat}{Are the coordinates in lon/lat (geographic) coordinates? +Set this to \code{FALSE} if you have projected data, e.g. with coordinates +representing distance in meters, not degrees. Lonlat coodinates are assumed +(\code{lonlat = TRUE}) is the default. In \code{slope_raster()} and \code{plot_slope()} +it is taken from the CRS of the routes (\code{sf::st_is_longlat(routes)}).} } \description{ -Extract slopes from xyz dataframe or sf objects +Extract slopes from xyz data frame or sf objects } \examples{ r_xyz = lisbon_road_segment_3d -slope_xyz(r_xyz) +slope_xyz(r_xyz, lonlat = FALSE) } diff --git a/vignettes/benchmark.Rmd b/vignettes/benchmark.Rmd new file mode 100644 index 0000000..7ee253f --- /dev/null +++ b/vignettes/benchmark.Rmd @@ -0,0 +1,84 @@ +--- +title: "Benchmarking slopes calculation" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Benchmarking slopes calculation} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + eval = FALSE +) +``` + +```{r setup} +library(slopes) +library(bench) +``` + +**Note: benchmarks do not currently evaluate on actions for faster build times.** + +# Performance + +For this benchmark we will download the following small (< 100 kB) `.tif` file: + +```{r} +u = "https://github.com/ITSLeeds/slopes/releases/download/0.0.0/dem_lisbon.tif" +f = basename(u) +if(!file.exists(f)) { + download.file(u, f, mode = "wb") +} +``` + +A benchmark can reveal how many route gradients can be calculated per second: + +```{r, results='hide'} +e = dem_lisbon_raster +r = lisbon_road_segments +et = terra::rast(f) +res = bench::mark(check = FALSE, + slope_raster = slope_raster(r, e), + slope_terra = slope_raster(r, et) +) +``` + +```{r} +res +``` + + +That is approximately + +```{r} +round(res$`itr/sec` * nrow(r)) +``` + +routes per second using the `raster` and `terra` (the default if installed, using `RasterLayer` and native `SpatRaster` objects) packages to extract elevation estimates from the raster datasets, respectively. + +The message: use the `terra` package to read-in DEM data for slope extraction if speed is important. + +To go faster, you can chose the `simple` method to gain some speed at the expense of accuracy: + +```{r, results='hide'} +e = dem_lisbon_raster +r = lisbon_road_segments +res = bench::mark(check = FALSE, + bilinear1 = slope_raster(r, e), + bilinear2 = slope_raster(r, et), + simple1 = slope_raster(r, e, method = "simple"), + simple2 = slope_raster(r, et, method = "simple") +) +``` + +```{r} +res +``` + + +```{r, include=FALSE} +file.remove("dem_lisbon.tif") +``` diff --git a/vignettes/roadnetworkcycling.Rmd b/vignettes/roadnetworkcycling.Rmd new file mode 100644 index 0000000..c3e7413 --- /dev/null +++ b/vignettes/roadnetworkcycling.Rmd @@ -0,0 +1,196 @@ +--- +title: "Example: gradients of a road network for a given city" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Example: gradients of a road network for a given city} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>", + fig.width=7, fig.height=5, fig.align = "center" +) +``` + +An example of the demand for data provided by the `slopes` package is a map showing gradients across São Paulo (Brazil, see image below), with a simplistic classification for cycling difficulty. + +![using slopes() to create a road network gradient for cycling for São Paulo (Brazil)](https://pbs.twimg.com/media/ErJ2dr8WMAIHwMn?format=jpg&name=small) + +This vignette will guide through the production of an interactive slope map for a road network, using `slopes`, `osmextract`, `sf`, `stplanr` and `tmap`. + +For the convenience of sample, we will use (Isle of Wight)[https://en.wikipedia.org/wiki/Isle_of_Wight] case, with 384 km2. See [Other examples] below. + +This will follow three steps: + +1. Download of road network from [OpenStreetMap](https://www.openstreetmap.org/) +2. Prepare the network +3. Compute slopes and export the map in html + + +## Extract the OSM network from geofabrik + +For this step you may use `osmextract` [package](https://itsleeds.github.io/osmextract/articles/osmextract.html) which downloads the most recent information available at OSM (https://download.geofabrik.de/index.html) and converts to _GeoPackage_ (.gpkg), the equivalent to _shapefile_. + +```{r setup1, warning=FALSE, message = FALSE} +library(dplyr) +library(sf) +# remotes::install_github("ITSLeeds/osmextract") +library(osmextract) +``` + +```{r get_iow, warning=FALSE, message = FALSE} +# get the network +iow_osm = oe_get("Isle of Wight", provider = "geofabrik", stringsAsFactors = FALSE, + quiet = FALSE, force_download = TRUE, force_vectortranslate = TRUE) # 7 MB + +# filter the major roads +iow_network = iow_osm %>% + dplyr::filter(highway %in% c('primary', "primary_link", 'secondary',"secondary_link", + 'tertiary', "tertiary_link", "trunk", "trunk_link", + "residential", "cycleway", "living_street", "unclassified", + "motorway", "motorway_link", "pedestrian", "steps", "track")) #remove: "service" +``` + +## Clean the road network + +These are optional steps that give better results, although they may slow down the process since they increase the number of segments present in the network. + +### Filter the unconnected segments + +The [`rnet_group()`](https://docs.ropensci.org/stplanr/reference/rnet_group.html) function from `stplanar` package assesses the connectivity of each segment assigns a group number (similar to a clustering process). Then we may filter the main group, the one with more connected segments. + +```{r setup2, warning=FALSE, message = FALSE} +# remotes::install_github("ropensci/stplanr") +library(stplanr) +``` + +```{r filter} +# filter unconnected roads +iow_network$group = rnet_group(iow_network) +iow_network_clean = iow_network %>% filter(group == 1) # the network with more connected segments +``` + +### Break the segments on vertices + +A very long segment will have an assigned average slope, but a very long segment can be broken into its nodes and have its own slope in each part of the segment. On one hand, we want the segments to break at their nodes. On the other hand, we don't want artificial *nodes* to be created where two lines cross, in particular where they have different **z** levels (_eg._ *brunels*: bridges and tunnels). + +The [`rnet_breakup_vertices`](https://docs.ropensci.org/stplanr/reference/rnet_breakup_vertices.html) from `stplanr` breaks the segments at their inner vertices, preserving the **brunels**. + +```{r breaking, warning=FALSE, message = FALSE} +iow_network_segments = rnet_breakup_vertices(iow_network_clean) +``` +In this case, there are more `r round(nrow(iow_network_segments)/nrow(iow_network_clean),2)`x segments than the original network. + +## Get slope values for each segment + +For this case we will use `slope_raster()` [function](https://itsleeds.github.io/slopes/reference/slope_raster.html), to retrieve the z values from a digital elevation model. This raster was obtained from STRM NASA mission. + +The **SRTM** (*Shuttle Radar Topography Mission*) NASA’s mission provides [freely available](https://gisgeography.com/srtm-shuttle-radar-topography-mission/) worldwide DEM, with a resolution of 25 to 30m and with a vertical accuracy of 16m - [more](https://www2.jpl.nasa.gov/srtm/). The resolution for USA might be better. + +Alternatively, **COPERNICUS** ESA's mission also provides [freely available](https://land.copernicus.eu/imagery-in-situ/eu-dem) DEM for all Europe, with a 25m resolution and a vertical accuracy of 7m - [more](https://land.copernicus.eu/user-corner/publications/eu-dem-flyer/view). + +Depending of how large is your road network, you can use `slope_3d()` [function](https://itsleeds.github.io/slopes/reference/slope_3d.html) - this will require a valid [Mapbox api key](https://docs.mapbox.com/api/overview/). + +```{r import_dem, message=FALSE} +# Import and plot DEM +u = "https://github.com/U-Shift/Declives-RedeViaria/releases/download/0.2/IsleOfWightNASA_clip.tif" +f = basename(u) +download.file(url = u, destfile = f, mode = "wb") +dem = raster::raster(f) +# res(dem) #27m of resolution +network = iow_network_segments + +library(raster) +plot(dem) +plot(sf::st_geometry(network), add = TRUE) #check if they overlay +``` + +All the required data is prepared to estimate the road segments' gradient. + +```{r slopes_values} +# Get the slope value for each segment (abs), using slopes package +library(slopes) +library(geodist) +network$slope = slope_raster(network, dem) +network$slope = network$slope*100 #percentage +summary(network$slope) #check the values +``` + +Half of the road segments in Isle of Wight have a gradient bellow `r round(median(network$slope),1)`%. + +We will adopt a simplistic qualitative classification for **cycling effort uphill**, and compare the number of segments in each class. + +```{r classify} +# Classify slopes +network$slope_class = network$slope %>% + cut( + breaks = c(0, 3, 5, 8, 10, 20, Inf), + labels = c("0-3: flat", "3-5: mild", "5-8: medium", "8-10: hard", + "10-20: extreme", ">20: impossible"), + right = F + ) +round(prop.table(table(network$slope_class))*100,1) +``` +It means that **`r round(prop.table(table(network$slope_class))[[1]]*100)`%** of the roads are flat or almost flat (0-3%) and about **`r round(prop.table(table(network$slope_class))[[1]]*100)+round(prop.table(table(network$slope_class))[[2]]*100)`%** of the roads are easily cyclable (0-5%). + +Now let us put this information on a map (see [here](https://rpubs.com/RobinLovelace/781081) for interactive version). + +```{r map, message = FALSE, eval=FALSE} +# more useful information +network$length = st_length(network) + +# make an interactive map +library(tmap) +palredgreen = c("#267300", "#70A800", "#FFAA00", "#E60000", "#A80000", "#730000") #color palette +# tmap_mode("view") +tmap_options(basemaps = leaflet::providers$CartoDB.Positron) #basemap + +slopemap = + tm_shape(network) + + tm_lines( + col = "slope_class", + palette = palredgreen, + lwd = 2, #line width + title.col = "Slope [%]", + popup.vars = c("Highway" = "highway", + "Length" = "length", + "Slope: " = "slope", + "Class: " = "slope_class"), + popup.format = list(digits = 1), + # id = "slope" + id = "name" #if it gets too memory consuming, delete this line + ) + +slopemap +``` + + +```{r export, echo=TRUE, eval=FALSE} +#export to html +tmap_save(slopemap, "html/SlopesIoW.html") + +# export information as geopackage +st_write(network, "shapefiles/SlopesIoW.gpkg", append=F) +``` + +#### Result: + +- [Isle of Wight (UK)](http://web.tecnico.ulisboa.pt/~rosamfelix/gis/declives/SlopesIoW.html) + +## Other examples + +- [São Paulo (Brazil)](http://web.tecnico.ulisboa.pt/~rosamfelix/gis/declives/DeclivesSaoPaulo.html) +- [Lisbon (Portugal)](http://web.tecnico.ulisboa.pt/~rosamfelix/gis/declives/DeclivesLisboa.html) +- [Oporto (Portugal)](http://web.tecnico.ulisboa.pt/~rosamfelix/gis/declives/DeclivesPorto_EU.html) +- Bristol (UK) +- Sheffield (UK) + + +```{r tidyup, include=FALSE} +rm(iow_osm,iow_network_clean,iow_network_segments, iow_network, slopemap) +file.remove(f) # remove the file, tidy up +``` + diff --git a/vignettes/slopes.Rmd b/vignettes/slopes.Rmd index eb5e934..c4cf946 100644 --- a/vignettes/slopes.Rmd +++ b/vignettes/slopes.Rmd @@ -1,8 +1,8 @@ --- -title: "An introduction to slopes" +title: "An introduction to the slopes package" output: bookdown::html_vignette2 vignette: > - %\VignetteIndexEntry{An introdution to slopes} + %\VignetteIndexEntry{An introduction to the slopes package} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} bibliography: slope-references.bib