From af3a7b8b0abfb04a13431d54e3e3d501faf04564 Mon Sep 17 00:00:00 2001 From: RFlx Date: Mon, 22 Jun 2026 13:43:55 +0100 Subject: [PATCH 1/4] extend elevation_add to support point geometries with configurable Z-coordinate embedding and elevation columns #21 --- R/slopes.R | 103 ++++++++++++++++++++++++++++++++------- man/elevation_add.Rd | 63 +++++++++++++++++++----- man/route_to_segments.Rd | 6 +-- 3 files changed, 139 insertions(+), 33 deletions(-) diff --git a/R/slopes.R b/R/slopes.R index a6b3934..2cb5b59 100644 --- a/R/slopes.R +++ b/R/slopes.R @@ -230,42 +230,110 @@ elevation_extract <- function(m, dem, method = "bilinear", terra = NULL) { terra::extract(dem, m[, 1:2], method = method)[[1]] } -#' Add elevation data to route linestrings +#' Add elevation data to routes linestrings or points #' -#' Adds elevation (Z) coordinates to linestring geometries using DEM data. +#' Adds elevation data to sf objects using a digital elevation model (DEM). #' -#' @param routes An sf object containing linestring geometries -#' @param dem A SpatRaster object containing elevation data (default: NULL for automatic download) -#' @param method Method for raster extraction (default: "bilinear") +#' For **linestring** geometries, the function attaches elevation as the Z +#' coordinate of each vertex, returning an XYZ linestring sf object. +#' +#' For **point** geometries, the function embeds the Z coordinate directly in +#' the geometry by default (returning POINT Z features), consistent with the +#' linestring behaviour. An `elevation` column can additionally be added by +#' setting `add_column = TRUE`. +#' +#' @param routes An sf object containing linestring or point geometries. +#' @param dem A SpatRaster object containing elevation data +#' (default: `NULL` for automatic download via `elevation_get()`). +#' @param method Method for raster extraction (default: `"bilinear"`). +#' @param add_z For **point** geometries only: if `TRUE` (the default), the Z +#' coordinate is embedded directly in the point geometry, returning POINT Z +#' features. Set to `FALSE` to keep the original XY geometry. Ignored for +#' linestrings (Z is always added to linestring vertices). +#' @param add_column For **point** geometries only: if `TRUE`, an `elevation` +#' column is added to the returned sf object in addition to (or instead of, +#' when `add_z = FALSE`) the Z geometry. Default: `FALSE`. Ignored for +#' linestrings. #' @param terra Deprecated. Ignored; terra is always used. -#' @return An sf object with XYZ linestring geometries +#' @return An sf object. For linestrings: XYZ linestring geometries. For points: +#' POINT Z geometry (when `add_z = TRUE`) and/or an `elevation` column (when +#' `add_column = TRUE`). #' @export #' @examples #' library(sf) -#' routes = lisbon_road_network[204, ] -#' dem = dem_lisbon() -#' (r3d = elevation_add(routes, dem)) -#' st_z_range(routes) +#' # Linestring usage: +#' routes <- lisbon_road_network[204, ] +#' dem <- dem_lisbon() +#' (r3d <- elevation_add(routes, dem)) #' st_z_range(r3d) #' plot(st_coordinates(r3d)[, 3]) #' plot_slope(r3d) +#' # Point usage — Z embedded in geometry by default: +#' pts <- sf::st_cast(sf::st_geometry(lisbon_road_network[204, ]), "POINT") +#' pts <- sf::st_sf(id = seq_along(pts), geometry = pts) +#' (pts_z <- elevation_add(pts, dem)) +#' sf::st_z_range(pts_z) +#' # Also add an elevation column: +#' (pts_z_col <- elevation_add(pts, dem, add_column = TRUE)) +#' pts_z_col$elevation +#' # Only an elevation column, keep XY geometry: +#' (pts_col <- elevation_add(pts, dem, add_z = FALSE, add_column = TRUE)) #' \dontrun{ #' # Get elevation data (requires internet connection, ceramic pkg, and API key): #' if (requireNamespace("ceramic", quietly = TRUE)) { -#' r3d_get = elevation_add(cyclestreets_route) +#' r3d_get <- elevation_add(cyclestreets_route) #' plot_slope(r3d_get) #' } #' } -elevation_add <- function(routes, dem = NULL, method = "bilinear", terra = NULL) { +elevation_add <- function(routes, dem = NULL, method = "bilinear", + add_z = TRUE, add_column = FALSE, + terra = NULL) { if (!is.null(terra)) { .Deprecated(msg = "The 'terra' argument is deprecated and ignored. terra is always used.") } stopifnotsf(routes) + + geom_types <- unique(sf::st_geometry_type(routes)) + + # ── POINT branch ──────────────────────────────────────────────────────────── + if (all(geom_types %in% c("POINT", "MULTIPOINT"))) { + if (is.null(dem)) { + dem <- elevation_get(routes) + r_original <- routes + routes <- sf::st_transform(routes, terra::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, dem, method = method)) + routes <- r_original # restore original CRS / coords + } else { + m <- sf::st_coordinates(routes) + z <- as.numeric(elevation_extract(m, dem, method = method)) + mo <- m + } + if (add_column) { + routes$elevation <- z + } + if (add_z) { + m_xyz <- cbind(mo[, 1:2], z) + points_z <- lapply(seq_len(nrow(m_xyz)), function(i) { + sf::st_point(m_xyz[i, ], dim = "XYZ") + }) + sf::st_geometry(routes) <- sf::st_sfc(points_z, crs = sf::st_crs(routes)) + } + return(routes) + } + + # ── LINESTRING branch ──────────────────────────────────────────────────────── if (is.null(dem)) { - dem <- elevation_get(routes) # returns SpatRaster + dem <- elevation_get(routes) # returns SpatRaster r_original <- routes routes <- sf::st_transform(routes, terra::crs(dem)) - suppressWarnings({sf::st_crs(routes) <- sf::st_crs(r_original)}) + 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, dem, method = method)) @@ -299,9 +367,9 @@ stopifnotsf <- function(x, arg_name = "routes") if (!methods::is(x, "sf")) stop( #' @return An sf object with one LINESTRING feature per vertex-to-vertex segment. #' @export #' @examples -#' route_xyz = elevation_add(lisbon_route, dem = dem_lisbon()) -#' segs = route_to_segments(route_xyz) -#' segs$slope = slope_xyz(segs) +#' route_xyz <- elevation_add(lisbon_route, dem = dem_lisbon()) +#' segs <- route_to_segments(route_xyz) +#' segs$slope <- slope_xyz(segs) #' summary(segs$slope) route_to_segments <- function(route_xyz) { coords <- sf::st_coordinates(route_xyz) @@ -309,4 +377,3 @@ route_to_segments <- function(route_xyz) { segs <- lapply(seq_len(n - 1), function(i) sf::st_linestring(coords[i:(i + 1), 1:3])) sf::st_sf(geometry = sf::st_sfc(segs, crs = sf::st_crs(route_xyz))) } - diff --git a/man/elevation_add.Rd b/man/elevation_add.Rd index 2e3b863..b2efeea 100644 --- a/man/elevation_add.Rd +++ b/man/elevation_add.Rd @@ -2,38 +2,77 @@ % Please edit documentation in R/slopes.R \name{elevation_add} \alias{elevation_add} -\title{Add elevation data to route linestrings} +\title{Add elevation data to routes linestrings or points} \usage{ -elevation_add(routes, dem = NULL, method = "bilinear", terra = NULL) +elevation_add( + routes, + dem = NULL, + method = "bilinear", + add_z = TRUE, + add_column = FALSE, + terra = NULL +) } \arguments{ -\item{routes}{An sf object containing linestring geometries} +\item{routes}{An sf object containing linestring or point geometries.} -\item{dem}{A SpatRaster object containing elevation data (default: NULL for automatic download)} +\item{dem}{A SpatRaster object containing elevation data +(default: \code{NULL} for automatic download via \code{elevation_get()}).} -\item{method}{Method for raster extraction (default: "bilinear")} +\item{method}{Method for raster extraction (default: \code{"bilinear"}).} + +\item{add_z}{For \strong{point} geometries only: if \code{TRUE} (the default), the Z +coordinate is embedded directly in the point geometry, returning POINT Z +features. Set to \code{FALSE} to keep the original XY geometry. Ignored for +linestrings (Z is always added to linestring vertices).} + +\item{add_column}{For \strong{point} geometries only: if \code{TRUE}, an \code{elevation} +column is added to the returned sf object in addition to (or instead of, +when \code{add_z = FALSE}) the Z geometry. Default: \code{FALSE}. Ignored for +linestrings.} \item{terra}{Deprecated. Ignored; terra is always used.} } \value{ -An sf object with XYZ linestring geometries +An sf object. For linestrings: XYZ linestring geometries. For points: +POINT Z geometry (when \code{add_z = TRUE}) and/or an \code{elevation} column (when +\code{add_column = TRUE}). } \description{ -Adds elevation (Z) coordinates to linestring geometries using DEM data. +Adds elevation data to sf objects using a digital elevation model (DEM). +} +\details{ +For \strong{linestring} geometries, the function attaches elevation as the Z +coordinate of each vertex, returning an XYZ linestring sf object. + +For \strong{point} geometries, the function embeds the Z coordinate directly in +the geometry by default (returning POINT Z features), consistent with the +linestring behaviour. An \code{elevation} column can additionally be added by +setting \code{add_column = TRUE}. } \examples{ library(sf) -routes = lisbon_road_network[204, ] -dem = dem_lisbon() -(r3d = elevation_add(routes, dem)) -st_z_range(routes) +# Linestring usage: +routes <- lisbon_road_network[204, ] +dem <- dem_lisbon() +(r3d <- elevation_add(routes, dem)) st_z_range(r3d) plot(st_coordinates(r3d)[, 3]) plot_slope(r3d) +# Point usage — Z embedded in geometry by default: +pts <- sf::st_cast(sf::st_geometry(lisbon_road_network[204, ]), "POINT") +pts <- sf::st_sf(id = seq_along(pts), geometry = pts) +(pts_z <- elevation_add(pts, dem)) +sf::st_z_range(pts_z) +# Also add an elevation column: +(pts_z_col <- elevation_add(pts, dem, add_column = TRUE)) +pts_z_col$elevation +# Only an elevation column, keep XY geometry: +(pts_col <- elevation_add(pts, dem, add_z = FALSE, add_column = TRUE)) \dontrun{ # Get elevation data (requires internet connection, ceramic pkg, and API key): if (requireNamespace("ceramic", quietly = TRUE)) { - r3d_get = elevation_add(cyclestreets_route) + r3d_get <- elevation_add(cyclestreets_route) plot_slope(r3d_get) } } diff --git a/man/route_to_segments.Rd b/man/route_to_segments.Rd index 3d9ab22..541aa7d 100644 --- a/man/route_to_segments.Rd +++ b/man/route_to_segments.Rd @@ -19,8 +19,8 @@ one per consecutive vertex pair. Useful for computing per-segment slopes with \code{\link[=slope_xyz]{slope_xyz()}}. } \examples{ -route_xyz = elevation_add(lisbon_route, dem = dem_lisbon()) -segs = route_to_segments(route_xyz) -segs$slope = slope_xyz(segs) +route_xyz <- elevation_add(lisbon_route, dem = dem_lisbon()) +segs <- route_to_segments(route_xyz) +segs$slope <- slope_xyz(segs) summary(segs$slope) } From 5ef3176e12fefffe2fcbcb278ed071d18d514ea1 Mon Sep 17 00:00:00 2001 From: RFlx Date: Mon, 22 Jun 2026 13:44:13 +0100 Subject: [PATCH 2/4] #21 add tests with elevation to points --- tests/testthat/Rplots.pdf | Bin 3611 -> 3611 bytes tests/testthat/test-slopes.R | 37 +++++++++++++++++++++++++++++++++++ 2 files changed, 37 insertions(+) diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf index 2835fe94d255fe17c60cb395fe98a73d1778ff8d..45f1c389d5edf9d04b31ed04cb16eebea237e73b 100644 GIT binary patch delta 47 wcmbO&Gh1eYxr&jIp|P=rv8g7PzHfetOJYf?f`*Hgk%5t!ff-D0WAp)D04L21UjP6A delta 47 wcmbO&Gh1eYxr(8gp}Dbvk&z~szHfetOJYf?f`*Hgk%5t!ff-D0WAp)D04FaDRR910 diff --git a/tests/testthat/test-slopes.R b/tests/testthat/test-slopes.R index a763c4b..01a1f08 100644 --- a/tests/testthat/test-slopes.R +++ b/tests/testthat/test-slopes.R @@ -86,3 +86,40 @@ test_that("slope_* functions work", { expect_error(slopes:::stop_is_not_linestring(1)) }) + +test_that("elevation_add() works with POINT geometries", { + dem <- dem_lisbon() + pts <- sf::st_cast(sf::st_geometry(lisbon_road_network[204, ]), "POINT") + pts <- sf::st_sf(id = seq_along(pts), geometry = pts) + + # Default: add_z = TRUE, no elevation column + pts_z <- elevation_add(pts, dem) + expect_s3_class(pts_z, "sf") + expect_equal(as.character(sf::st_geometry_type(pts_z)[1]), "POINT") + expect_true(!is.na(sf::st_z_range(pts_z)[1])) # has Z coordinate + expect_false("elevation" %in% names(pts_z)) + expect_equal( + sf::st_z_range(pts_z), + c(86, 92), + ignore_attr = TRUE, + tolerance = 10 + ) + + # add_column = TRUE: Z in geometry AND elevation column + pts_z_col <- elevation_add(pts, dem, add_column = TRUE) + expect_true("elevation" %in% names(pts_z_col)) + expect_true(!is.na(sf::st_z_range(pts_z_col)[1])) # still has Z + expect_equal( + round(pts_z_col$elevation[1:3], 2), + c(92.31, 91.93, 91.60) + ) + + # add_z = FALSE, add_column = TRUE: XY geometry with elevation column + pts_col <- elevation_add(pts, dem, add_z = FALSE, add_column = TRUE) + expect_true("elevation" %in% names(pts_col)) + expect_true(is.null(sf::st_z_range(pts_col))) # XY: no Z coordinate + expect_equal( + round(pts_col$elevation[1:3], 2), + c(92.31, 91.93, 91.60) + ) +}) From e1a50ad74b71d9e83d4d48e2c0cf6ef2225f11cf Mon Sep 17 00:00:00 2001 From: RFlx Date: Mon, 22 Jun 2026 13:57:54 +0100 Subject: [PATCH 3/4] reorganize elevation functions in elevation.R module and update related tests --- R/elevation.R | 245 ++++++++++++++++++++++++++++++++ R/slope_get.R | 44 +----- R/slopes.R | 149 +------------------ R/z-functions.R | 32 ----- tests/testthat/Rplots.pdf | Bin 3611 -> 3611 bytes tests/testthat/test-elevation.R | 88 ++++++++++++ tests/testthat/test-get.R | 8 -- tests/testthat/test-slopes.R | 57 +------- tests/testthat/test-z.R | 16 --- 9 files changed, 336 insertions(+), 303 deletions(-) create mode 100644 R/elevation.R delete mode 100644 R/z-functions.R create mode 100644 tests/testthat/test-elevation.R delete mode 100644 tests/testthat/test-get.R delete mode 100644 tests/testthat/test-z.R diff --git a/R/elevation.R b/R/elevation.R new file mode 100644 index 0000000..52bbec3 --- /dev/null +++ b/R/elevation.R @@ -0,0 +1,245 @@ +# Elevation get/extract/add and Z helper functions + +# ── Z value helpers ───────────────────────────────────────────────────────── + +z_value <- function(x) { + coords <- sf::st_coordinates(x) + if ("Z" %in% colnames(coords)) { + return(coords[, "Z"]) + } else { + stop("No Z coordinates found in the input data") + } +} + +z_mean <- function(x) { + mean(z_value(x), na.rm = TRUE) +} + +z_min <- function(x) { + min(z_value(x), na.rm = TRUE) +} + +z_max <- function(x) { + max(z_value(x), na.rm = TRUE) +} + +z_start <- function(x) { + z_vals <- z_value(x) + z_vals[1] +} + +z_end <- function(x) { + z_vals <- z_value(x) + z_vals[length(z_vals)] +} + +z_elevation_change_start_end <- function(x) { + z_end(x) - z_start(x) +} + +z_direction <- function(x) { + sign(z_elevation_change_start_end(x)) +} + +z_cumulative_difference <- function(x) { + z <- z_value(x) + sum(abs(diff(z)), na.rm = TRUE) +} + +# ── elevation_get ──────────────────────────────────────────────────────────── + +#' Get elevation data for routes +#' +#' Downloads elevation data using the ceramic package for given routes. +#' Returns a `SpatRaster` object (terra package). +#' +#' @param routes An sf object containing linestring geometries +#' @param ... Additional arguments passed to ceramic::cc_elevation +#' @return A SpatRaster covering the routes +#' @export +elevation_get <- function(routes, ...) { + if (!requireNamespace("ceramic", quietly = TRUE)) { + stop("Install the package ceramic to use elevation_get().") + } + mid_ext <- sf_mid_ext_lonlat(routes) + bw <- max(c(mid_ext$width, mid_ext$height)) / 1 # buffer width + suppressWarnings({ + e <- ceramic::cc_elevation(loc = mid_ext$midpoint, buffer = bw, ...) + }) + crs_routes <- sf::st_crs(routes) + terra::project(e, y = crs_routes$wkt) +} + +#' Extract midpoint and extent from routes in lonlat +#' +#' Internal helper function to get midpoint and extent of routes in lon/lat +#' coordinates. +#' +#' @param routes An sf object containing linestring geometries +#' @return A list with midpoint coordinates and width/height dimensions +sf_mid_ext_lonlat <- function(routes) { + res <- list() + if (!sf::st_is_longlat(routes)) { + routes <- sf::st_transform(routes, 4326) + } + 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( + c(x = bb[1], y = bb[2]), + c(x = bb[1], y = bb[4]) + ) + res +} + +# ── elevation_extract ──────────────────────────────────────────────────────── + +#' Extract elevation values from coordinates +#' +#' Extracts elevation values from a DEM raster at specified coordinate locations. +#' Accepts both `SpatRaster` (terra) and legacy `Raster*` (raster) objects; +#' legacy objects are automatically converted to `SpatRaster`. +#' +#' @param m Matrix or sf object with coordinates +#' @param dem A SpatRaster (or legacy RasterLayer) containing elevation data +#' @param method Method for raster extraction (default: "bilinear") +#' @param terra Deprecated. Ignored; terra is always used. +#' @return Numeric vector of elevation values +#' @export +elevation_extract <- function(m, dem, method = "bilinear", terra = NULL) { + if (!is.null(terra)) { + .Deprecated(msg = "The 'terra' argument is deprecated and ignored. terra is always used.") + } + if (any(grepl(pattern = "sf", class(m)))) m <- sf::st_coordinates(m) + if (!methods::is(dem, "SpatRaster")) { + if (requireNamespace("terra", quietly = TRUE)) { + message("Converting legacy Raster* object to SpatRaster. Consider using terra::rast() directly.") + dem <- terra::rast(dem) + } else { + stop("terra package is required. Install it with: install.packages('terra')") + } + } + terra::extract(dem, m[, 1:2], method = method)[[1]] +} + +# ── elevation_add ──────────────────────────────────────────────────────────── + +#' Add elevation data to routes linestrings or points +#' +#' Adds elevation data to sf objects using a digital elevation model (DEM). +#' +#' For **linestring** geometries, the function attaches elevation as the Z +#' coordinate of each vertex, returning an XYZ linestring sf object. +#' +#' For **point** geometries, the function embeds the Z coordinate directly in +#' the geometry by default (returning POINT Z features), consistent with the +#' linestring behaviour. An `elevation` column can additionally be added by +#' setting `add_column = TRUE`. +#' +#' @param routes An sf object containing linestring or point geometries. +#' @param dem A SpatRaster object containing elevation data +#' (default: `NULL` for automatic download via `elevation_get()`). +#' @param method Method for raster extraction (default: `"bilinear"`). +#' @param add_z For **point** geometries only: if `TRUE` (the default), the Z +#' coordinate is embedded directly in the point geometry, returning POINT Z +#' features. Set to `FALSE` to keep the original XY geometry. Ignored for +#' linestrings (Z is always added to linestring vertices). +#' @param add_column For **point** geometries only: if `TRUE`, an `elevation` +#' column is added to the returned sf object in addition to (or instead of, +#' when `add_z = FALSE`) the Z geometry. Default: `FALSE`. Ignored for +#' linestrings. +#' @param terra Deprecated. Ignored; terra is always used. +#' @return An sf object. For linestrings: XYZ linestring geometries. For points: +#' POINT Z geometry (when `add_z = TRUE`) and/or an `elevation` column (when +#' `add_column = TRUE`). +#' @export +#' @examples +#' library(sf) +#' # Linestring usage: +#' routes <- lisbon_road_network[204, ] +#' dem <- dem_lisbon() +#' (r3d <- elevation_add(routes, dem)) +#' st_z_range(r3d) +#' plot(st_coordinates(r3d)[, 3]) +#' plot_slope(r3d) +#' # Point usage — Z embedded in geometry by default: +#' pts <- sf::st_cast(sf::st_geometry(lisbon_road_network[204, ]), "POINT") +#' pts <- sf::st_sf(id = seq_along(pts), geometry = pts) +#' (pts_z <- elevation_add(pts, dem)) +#' sf::st_z_range(pts_z) +#' # Also add an elevation column: +#' (pts_z_col <- elevation_add(pts, dem, add_column = TRUE)) +#' pts_z_col$elevation +#' # Only an elevation column, keep XY geometry: +#' (pts_col <- elevation_add(pts, dem, add_z = FALSE, add_column = TRUE)) +#' \dontrun{ +#' # Get elevation data (requires internet connection, ceramic pkg, and API key): +#' if (requireNamespace("ceramic", quietly = TRUE)) { +#' r3d_get <- elevation_add(cyclestreets_route) +#' plot_slope(r3d_get) +#' } +#' } +elevation_add <- function(routes, dem = NULL, method = "bilinear", + add_z = TRUE, add_column = FALSE, + terra = NULL) { + if (!is.null(terra)) { + .Deprecated(msg = "The 'terra' argument is deprecated and ignored. terra is always used.") + } + stopifnotsf(routes) + + geom_types <- unique(sf::st_geometry_type(routes)) + + # ── POINT branch ──────────────────────────────────────────────────────────── + if (all(geom_types %in% c("POINT", "MULTIPOINT"))) { + if (is.null(dem)) { + dem <- elevation_get(routes) + r_original <- routes + routes <- sf::st_transform(routes, terra::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, dem, method = method)) + routes <- r_original # restore original CRS / coords + } else { + m <- sf::st_coordinates(routes) + z <- as.numeric(elevation_extract(m, dem, method = method)) + mo <- m + } + if (add_column) { + routes$elevation <- z + } + if (add_z) { + m_xyz <- cbind(mo[, 1:2], z) + points_z <- lapply(seq_len(nrow(m_xyz)), function(i) { + sf::st_point(m_xyz[i, ], dim = "XYZ") + }) + sf::st_geometry(routes) <- sf::st_sfc(points_z, crs = sf::st_crs(routes)) + } + return(routes) + } + + # ── LINESTRING branch ──────────────────────────────────────────────────────── + if (is.null(dem)) { + dem <- elevation_get(routes) # returns SpatRaster + r_original <- routes + routes <- sf::st_transform(routes, terra::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, dem, method = method)) + m_xyz <- cbind(mo[, 1:2], z) + } else { + m <- sf::st_coordinates(routes) + z <- as.numeric(elevation_extract(m, dem, method = method)) + m_xyz <- cbind(m[, 1:2], z) + } + n <- nrow(routes) + linestrings <- lapply(seq(n), function(i) sf::st_linestring(m_xyz[m[, "L1"] == i, ])) + rgeom3d_sfc <- sf::st_sfc(linestrings, crs = sf::st_crs(routes)) + sf::st_geometry(routes) <- rgeom3d_sfc + routes +} diff --git a/R/slope_get.R b/R/slope_get.R index 7b7d57e..c7e5ed5 100644 --- a/R/slope_get.R +++ b/R/slope_get.R @@ -1,46 +1,4 @@ -#' Get elevation data for routes -#' -#' Downloads elevation data using the ceramic package for given routes. -#' Returns a `SpatRaster` object (terra package). -#' -#' @param routes An sf object containing linestring geometries -#' @param ... Additional arguments passed to ceramic::cc_elevation -#' @return A SpatRaster covering the routes -#' @export -elevation_get = function(routes, ...) { - if (!requireNamespace("ceramic", quietly = TRUE)) { - stop("Install the package ceramic to use elevation_get().") - } - mid_ext = sf_mid_ext_lonlat(routes) - bw = max(c(mid_ext$width, mid_ext$height)) / 1 # buffer width - suppressWarnings({ - e = ceramic::cc_elevation(loc = mid_ext$midpoint, buffer = bw, ...) - }) - crs_routes = sf::st_crs(routes) - terra::project(e, y = crs_routes$wkt) -} - -#' Extract midpoint and extent from routes in lonlat -#' -#' Internal helper function to get midpoint and extent of routes in lon/lat coordinates. -#' -#' @param routes An sf object containing linestring geometries -#' @return A list with midpoint coordinates and width/height dimensions -sf_mid_ext_lonlat = function(routes) { - res = list() - if(!sf::st_is_longlat(routes)) { - routes = sf::st_transform(routes, 4326) - } - 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( - c(x = bb[1], y = bb[2]), - c(x = bb[1], y = bb[4]) - ) - res -} - +# Raster / matrix utility functions for slope data #' Convert slope matrix to SpatRaster #' #' Converts a slope matrix or a legacy RasterLayer to a SpatRaster (terra). diff --git a/R/slopes.R b/R/slopes.R index 2cb5b59..b1c1e4c 100644 --- a/R/slopes.R +++ b/R/slopes.R @@ -202,154 +202,6 @@ slope_xyz <- function(route_xyz, fun = slope_matrix_weighted, lonlat = TRUE, dir } } -#' Extract elevation values from coordinates -#' -#' Extracts elevation values from a DEM raster at specified coordinate locations. -#' Accepts both `SpatRaster` (terra) and legacy `Raster*` (raster) objects; -#' legacy objects are automatically converted to `SpatRaster`. -#' -#' @param m Matrix or sf object with coordinates -#' @param dem A SpatRaster (or legacy RasterLayer) containing elevation data -#' @param method Method for raster extraction (default: "bilinear") -#' @param terra Deprecated. Ignored; terra is always used. -#' @return Numeric vector of elevation values -#' @export -elevation_extract <- function(m, dem, method = "bilinear", terra = NULL) { - if (!is.null(terra)) { - .Deprecated(msg = "The 'terra' argument is deprecated and ignored. terra is always used.") - } - if (any(grepl(pattern = "sf", class(m)))) m <- sf::st_coordinates(m) - if (!methods::is(dem, "SpatRaster")) { - if (requireNamespace("terra", quietly = TRUE)) { - message("Converting legacy Raster* object to SpatRaster. Consider using terra::rast() directly.") - dem <- terra::rast(dem) - } else { - stop("terra package is required. Install it with: install.packages('terra')") - } - } - terra::extract(dem, m[, 1:2], method = method)[[1]] -} - -#' Add elevation data to routes linestrings or points -#' -#' Adds elevation data to sf objects using a digital elevation model (DEM). -#' -#' For **linestring** geometries, the function attaches elevation as the Z -#' coordinate of each vertex, returning an XYZ linestring sf object. -#' -#' For **point** geometries, the function embeds the Z coordinate directly in -#' the geometry by default (returning POINT Z features), consistent with the -#' linestring behaviour. An `elevation` column can additionally be added by -#' setting `add_column = TRUE`. -#' -#' @param routes An sf object containing linestring or point geometries. -#' @param dem A SpatRaster object containing elevation data -#' (default: `NULL` for automatic download via `elevation_get()`). -#' @param method Method for raster extraction (default: `"bilinear"`). -#' @param add_z For **point** geometries only: if `TRUE` (the default), the Z -#' coordinate is embedded directly in the point geometry, returning POINT Z -#' features. Set to `FALSE` to keep the original XY geometry. Ignored for -#' linestrings (Z is always added to linestring vertices). -#' @param add_column For **point** geometries only: if `TRUE`, an `elevation` -#' column is added to the returned sf object in addition to (or instead of, -#' when `add_z = FALSE`) the Z geometry. Default: `FALSE`. Ignored for -#' linestrings. -#' @param terra Deprecated. Ignored; terra is always used. -#' @return An sf object. For linestrings: XYZ linestring geometries. For points: -#' POINT Z geometry (when `add_z = TRUE`) and/or an `elevation` column (when -#' `add_column = TRUE`). -#' @export -#' @examples -#' library(sf) -#' # Linestring usage: -#' routes <- lisbon_road_network[204, ] -#' dem <- dem_lisbon() -#' (r3d <- elevation_add(routes, dem)) -#' st_z_range(r3d) -#' plot(st_coordinates(r3d)[, 3]) -#' plot_slope(r3d) -#' # Point usage — Z embedded in geometry by default: -#' pts <- sf::st_cast(sf::st_geometry(lisbon_road_network[204, ]), "POINT") -#' pts <- sf::st_sf(id = seq_along(pts), geometry = pts) -#' (pts_z <- elevation_add(pts, dem)) -#' sf::st_z_range(pts_z) -#' # Also add an elevation column: -#' (pts_z_col <- elevation_add(pts, dem, add_column = TRUE)) -#' pts_z_col$elevation -#' # Only an elevation column, keep XY geometry: -#' (pts_col <- elevation_add(pts, dem, add_z = FALSE, add_column = TRUE)) -#' \dontrun{ -#' # Get elevation data (requires internet connection, ceramic pkg, and API key): -#' if (requireNamespace("ceramic", quietly = TRUE)) { -#' r3d_get <- elevation_add(cyclestreets_route) -#' plot_slope(r3d_get) -#' } -#' } -elevation_add <- function(routes, dem = NULL, method = "bilinear", - add_z = TRUE, add_column = FALSE, - terra = NULL) { - if (!is.null(terra)) { - .Deprecated(msg = "The 'terra' argument is deprecated and ignored. terra is always used.") - } - stopifnotsf(routes) - - geom_types <- unique(sf::st_geometry_type(routes)) - - # ── POINT branch ──────────────────────────────────────────────────────────── - if (all(geom_types %in% c("POINT", "MULTIPOINT"))) { - if (is.null(dem)) { - dem <- elevation_get(routes) - r_original <- routes - routes <- sf::st_transform(routes, terra::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, dem, method = method)) - routes <- r_original # restore original CRS / coords - } else { - m <- sf::st_coordinates(routes) - z <- as.numeric(elevation_extract(m, dem, method = method)) - mo <- m - } - if (add_column) { - routes$elevation <- z - } - if (add_z) { - m_xyz <- cbind(mo[, 1:2], z) - points_z <- lapply(seq_len(nrow(m_xyz)), function(i) { - sf::st_point(m_xyz[i, ], dim = "XYZ") - }) - sf::st_geometry(routes) <- sf::st_sfc(points_z, crs = sf::st_crs(routes)) - } - return(routes) - } - - # ── LINESTRING branch ──────────────────────────────────────────────────────── - if (is.null(dem)) { - dem <- elevation_get(routes) # returns SpatRaster - r_original <- routes - routes <- sf::st_transform(routes, terra::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, dem, method = method)) - m_xyz <- cbind(mo[, 1:2], z) - } else { - m <- sf::st_coordinates(routes) - z <- as.numeric(elevation_extract(m, dem, method = method)) - m_xyz <- cbind(m[, 1:2], z) - } - n <- nrow(routes) - linestrings <- lapply(seq(n), function(i) sf::st_linestring(m_xyz[m[, "L1"] == i, ])) - rgeom3d_sfc <- sf::st_sfc(linestrings, crs = sf::st_crs(routes)) - sf::st_geometry(routes) <- rgeom3d_sfc - routes -} - # Utility functions has_terra <- function() requireNamespace("terra", quietly = TRUE) is_linestring <- function(x) unique(sf::st_geometry_type(x)) == "LINESTRING" @@ -377,3 +229,4 @@ route_to_segments <- function(route_xyz) { segs <- lapply(seq_len(n - 1), function(i) sf::st_linestring(coords[i:(i + 1), 1:3])) sf::st_sf(geometry = sf::st_sfc(segs, crs = sf::st_crs(route_xyz))) } + diff --git a/R/z-functions.R b/R/z-functions.R deleted file mode 100644 index 89d9116..0000000 --- a/R/z-functions.R +++ /dev/null @@ -1,32 +0,0 @@ - -z_value <- function(x) { - coords <- sf::st_coordinates(x) - if("Z" %in% colnames(coords)) { - return(coords[, "Z"]) - } else { - stop("No Z coordinates found in the input data") - } -} - -z_mean <- function(x) { - mean(z_value(x), na.rm = TRUE) -} - -z_min <- function(x) { - min(z_value(x), na.rm = TRUE) -} - -z_max <- function(x) { - max(z_value(x), na.rm = TRUE) -} - -z_start <- function(x) { - z_vals <- z_value(x) - z_vals[1] -} - -z_end <- function(x) { - z_vals <- z_value(x) - z_vals[length(z_vals)] -} - diff --git a/tests/testthat/Rplots.pdf b/tests/testthat/Rplots.pdf index 45f1c389d5edf9d04b31ed04cb16eebea237e73b..6353f1376c8a6d0dcd2600c1f975febfe002f11b 100644 GIT binary patch delta 23 bcmbO&Gh1eYHH)dKk;O!N2{64e`T#EgPr?SP delta 23 bcmbO&Gh1eYHH)!@vFSv62{64e`T#EgPqGH7 diff --git a/tests/testthat/test-elevation.R b/tests/testthat/test-elevation.R new file mode 100644 index 0000000..4223301 --- /dev/null +++ b/tests/testthat/test-elevation.R @@ -0,0 +1,88 @@ +test_that("functions to get elevations from sf objects work", { + if(nchar(Sys.getenv("MAPBOX_API_KEY")) < 8) + skip(message = "Skipping test, MAPBOX token in .Renviron needed") + r = cyclestreets_route + e = elevation_get(r) + expect_true(methods::is(e, "SpatRaster")) +}) + +test_that("elevation_extract() works", { + m = sf::st_coordinates(lisbon_road_segment) + e = elevation_extract(m, dem_lisbon()) + expect_identical(round(e[1:3], 2), c(92.31, 91.93, 91.60)) + e = elevation_extract(lisbon_road_segment, dem_lisbon()) + expect_identical(round(e[1:3], 2), c(92.31, 91.93, 91.60)) +}) + +test_that("elevation_add() works with LINESTRING geometries", { + e = dem_lisbon() + r = lisbon_road_network[204, ] + r3d = elevation_add(r, e) + expect_equal( + sf::st_z_range(r3d$geom), + c(86, 92), + ignore_attr = TRUE, + tolerance = 10 + ) + if(nchar(Sys.getenv("MAPBOX_API_KEY")) < 8) + skip(message = "Skipping test, MAPBOX token in .Renviron needed") + r3d2 = elevation_add(r) + expect_equal( + sf::st_z_range(r3d2$geom), + c(86, 92), + ignore_attr = TRUE, + tolerance = 10 + ) +}) + +test_that("elevation_add() works with POINT geometries", { + dem <- dem_lisbon() + pts <- sf::st_cast(sf::st_geometry(lisbon_road_network[204, ]), "POINT") + pts <- sf::st_sf(id = seq_along(pts), geometry = pts) + + # Default: add_z = TRUE, no elevation column + pts_z <- elevation_add(pts, dem) + expect_s3_class(pts_z, "sf") + expect_equal(as.character(sf::st_geometry_type(pts_z)[1]), "POINT") + expect_true(!is.na(sf::st_z_range(pts_z)[1])) # has Z coordinate + expect_false("elevation" %in% names(pts_z)) + expect_equal( + sf::st_z_range(pts_z), + c(86, 92), + ignore_attr = TRUE, + tolerance = 10 + ) + + # add_column = TRUE: Z in geometry AND elevation column + pts_z_col <- elevation_add(pts, dem, add_column = TRUE) + expect_true("elevation" %in% names(pts_z_col)) + expect_true(!is.na(sf::st_z_range(pts_z_col)[1])) # still has Z + expect_equal( + round(pts_z_col$elevation[1:3], 2), + c(92.31, 91.93, 91.60) + ) + + # add_z = FALSE, add_column = TRUE: XY geometry with elevation column + pts_col <- elevation_add(pts, dem, add_z = FALSE, add_column = TRUE) + expect_true("elevation" %in% names(pts_col)) + expect_true(is.null(sf::st_z_range(pts_col))) # XY: no Z coordinate + expect_equal( + round(pts_col$elevation[1:3], 2), + c(92.31, 91.93, 91.60) + ) +}) + +test_that("Functions on Z values work", { + x = slopes::lisbon_route_3d + expect_equal(class(z_value(x)), "numeric") + x = slopes::lisbon_route + expect_error(z_value(x)) + expect_error(z_start(x)) + expect_error(z_end(x)) + expect_error(z_mean(x)) + expect_error(z_min(x)) + expect_error(z_max(x)) + expect_error(z_elevation_change_start_end(x)) + expect_error(z_direction(x)) + expect_error(z_cumulative_difference(x)) +}) diff --git a/tests/testthat/test-get.R b/tests/testthat/test-get.R deleted file mode 100644 index 3bec16d..0000000 --- a/tests/testthat/test-get.R +++ /dev/null @@ -1,8 +0,0 @@ -test_that("functions to get elevations from sf objects work", { - if(nchar(Sys.getenv("MAPBOX_API_KEY")) < 8) - skip(message = "Skipping test, MAPBOX token in .Renviron needed") - r = cyclestreets_route - e = elevation_get(r) - expect_true(methods::is(e, "SpatRaster")) -}) - diff --git a/tests/testthat/test-slopes.R b/tests/testthat/test-slopes.R index 01a1f08..a5fc1a8 100644 --- a/tests/testthat/test-slopes.R +++ b/tests/testthat/test-slopes.R @@ -13,9 +13,6 @@ test_that("slope_* functions work", { expect_true(sequential_right) e = elevation_extract(m, dem_lisbon()) - expect_identical(round(e[1:3], 2), c(92.31, 91.93, 91.60)) - e = elevation_extract(lisbon_road_segment, dem_lisbon()) - expect_identical(round(e[1:3], 2), c(92.31, 91.93, 91.60)) s = slope_distance(d, e) expect_identical(round(s[1:3], 3), c(-0.047, -0.041, -0.025)) @@ -64,62 +61,10 @@ test_that("slope_* functions work", { r_xyz = lisbon_road_segment_3d expect_equal(slope_xyz(r_xyz), 0.0950132312274622, ignore_attr = TRUE) - r = lisbon_road_network[204, ] - r3d = elevation_add(r, e) - expect_equal( - sf::st_z_range(r3d$geom), - c(86, 92), - ignore_attr = TRUE, - tolerance = 10 - ) - if(nchar(Sys.getenv("MAPBOX_API_KEY")) < 8) - skip(message = "Skipping test, MAPBOX token in .Renviron needed") - r3d2 = elevation_add(r) - expect_equal( - sf::st_z_range(r3d2$geom), - c(86, 92), - ignore_attr = TRUE, - tolerance = 10 - ) + expect_error(slopes:::stopifnotsf(1)) expect_error(slopes:::stop_is_not_linestring(1)) }) -test_that("elevation_add() works with POINT geometries", { - dem <- dem_lisbon() - pts <- sf::st_cast(sf::st_geometry(lisbon_road_network[204, ]), "POINT") - pts <- sf::st_sf(id = seq_along(pts), geometry = pts) - - # Default: add_z = TRUE, no elevation column - pts_z <- elevation_add(pts, dem) - expect_s3_class(pts_z, "sf") - expect_equal(as.character(sf::st_geometry_type(pts_z)[1]), "POINT") - expect_true(!is.na(sf::st_z_range(pts_z)[1])) # has Z coordinate - expect_false("elevation" %in% names(pts_z)) - expect_equal( - sf::st_z_range(pts_z), - c(86, 92), - ignore_attr = TRUE, - tolerance = 10 - ) - - # add_column = TRUE: Z in geometry AND elevation column - pts_z_col <- elevation_add(pts, dem, add_column = TRUE) - expect_true("elevation" %in% names(pts_z_col)) - expect_true(!is.na(sf::st_z_range(pts_z_col)[1])) # still has Z - expect_equal( - round(pts_z_col$elevation[1:3], 2), - c(92.31, 91.93, 91.60) - ) - - # add_z = FALSE, add_column = TRUE: XY geometry with elevation column - pts_col <- elevation_add(pts, dem, add_z = FALSE, add_column = TRUE) - expect_true("elevation" %in% names(pts_col)) - expect_true(is.null(sf::st_z_range(pts_col))) # XY: no Z coordinate - expect_equal( - round(pts_col$elevation[1:3], 2), - c(92.31, 91.93, 91.60) - ) -}) diff --git a/tests/testthat/test-z.R b/tests/testthat/test-z.R deleted file mode 100644 index 1c1f4e8..0000000 --- a/tests/testthat/test-z.R +++ /dev/null @@ -1,16 +0,0 @@ -test_that("Functions on Z values work", { - x = slopes::lisbon_route_3d - expect_equal(class(z_value(x)), "numeric") - x = slopes::lisbon_route - expect_error(z_value(x)) - expect_error(z_start(x)) - expect_error(z_end(x)) - expect_error(z_mean(x)) - expect_error(z_min(x)) - expect_error(z_max(x)) - expect_error(z_elevation_change_start_end(x)) - expect_error(z_direction(x)) - expect_error(z_cumulative_difference(x)) - # x_slope = - # expect_is() -}) From 01cf9d3aa706923c0705f7b4103a438df1c5a9c8 Mon Sep 17 00:00:00 2001 From: RFlx Date: Mon, 22 Jun 2026 14:02:49 +0100 Subject: [PATCH 4/4] update documentation for elevation_add with points #21 --- man/elevation_add.Rd | 2 +- man/elevation_extract.Rd | 2 +- man/elevation_get.Rd | 2 +- man/sf_mid_ext_lonlat.Rd | 5 +++-- vignettes/slopes.Rmd | 10 +++++++--- 5 files changed, 13 insertions(+), 8 deletions(-) diff --git a/man/elevation_add.Rd b/man/elevation_add.Rd index b2efeea..4b7e0bf 100644 --- a/man/elevation_add.Rd +++ b/man/elevation_add.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/slopes.R +% Please edit documentation in R/elevation.R \name{elevation_add} \alias{elevation_add} \title{Add elevation data to routes linestrings or points} diff --git a/man/elevation_extract.Rd b/man/elevation_extract.Rd index cbf8367..17b637b 100644 --- a/man/elevation_extract.Rd +++ b/man/elevation_extract.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/slopes.R +% Please edit documentation in R/elevation.R \name{elevation_extract} \alias{elevation_extract} \title{Extract elevation values from coordinates} diff --git a/man/elevation_get.Rd b/man/elevation_get.Rd index 072fb58..b0d353c 100644 --- a/man/elevation_get.Rd +++ b/man/elevation_get.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/slope_get.R +% Please edit documentation in R/elevation.R \name{elevation_get} \alias{elevation_get} \title{Get elevation data for routes} diff --git a/man/sf_mid_ext_lonlat.Rd b/man/sf_mid_ext_lonlat.Rd index 1b813cd..c5362d5 100644 --- a/man/sf_mid_ext_lonlat.Rd +++ b/man/sf_mid_ext_lonlat.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/slope_get.R +% Please edit documentation in R/elevation.R \name{sf_mid_ext_lonlat} \alias{sf_mid_ext_lonlat} \title{Extract midpoint and extent from routes in lonlat} @@ -13,5 +13,6 @@ sf_mid_ext_lonlat(routes) A list with midpoint coordinates and width/height dimensions } \description{ -Internal helper function to get midpoint and extent of routes in lon/lat coordinates. +Internal helper function to get midpoint and extent of routes in lon/lat +coordinates. } diff --git a/vignettes/slopes.Rmd b/vignettes/slopes.Rmd index d0dae5e..11f9438 100644 --- a/vignettes/slopes.Rmd +++ b/vignettes/slopes.Rmd @@ -101,7 +101,7 @@ data(lisbon_route) dem_lisbon <- dem_lisbon() ``` -## Add elevation to a linestring +## Add elevation to linestrings or points If you have a 2D linestring and a DEM, you can add elevation data to the linestring using `elevation_add()`: @@ -118,6 +118,8 @@ If you don't have a local DEM, `elevation_add()` can download elevation data (th # head(sf::st_coordinates(sf_linestring_xyz_mapbox)) ``` +You can also use `elevation_add()` with points. + ## Calculate slope Once you have a 3D linestring (with XYZ coordinates), you can calculate its average slope using `slope_xyz()`: @@ -157,7 +159,8 @@ summary(lisbon_route_segments_xyz$slope) # Segments are coloured by steepness (absolute slope), regardless of direction # (uphill or downhill). slope_breaks are in proportions, matching slope_xyz() output. col_idx <- cut(abs(lisbon_route_segments_xyz$slope), - breaks = slope_breaks, labels = FALSE, include.lowest = TRUE) + breaks = slope_breaks, labels = FALSE, include.lowest = TRUE +) plot(st_geometry(lisbon_route_segments_xyz), col = slope_colors[col_idx], lwd = 3, main = "Slope by vertex segments" @@ -178,7 +181,8 @@ summary(lisbon_route_100m_xyz$slope) ```{r} col_idx <- cut(abs(lisbon_route_100m_xyz$slope), - breaks = slope_breaks, labels = FALSE, include.lowest = TRUE) + breaks = slope_breaks, labels = FALSE, include.lowest = TRUE +) plot(st_geometry(lisbon_route_100m_xyz), col = slope_colors[col_idx], lwd = 3, main = "Slope by 100 m segments"