Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions NAMESPACE
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,7 @@ export(elevation_extract)
export(elevation_get)
export(plot_dz)
export(plot_slope)
export(route_to_segments)
export(sequential_dist)
export(slope_distance)
export(slope_distance_mean)
Expand Down
25 changes: 24 additions & 1 deletion R/slopes.R
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,7 @@ elevation_add <- function(routes, dem = NULL, method = "bilinear", terra = NULL)
m_xyz <- cbind(m[, 1:2], z)
}
n <- nrow(routes)
linestrings <- lapply(seq(n), function(i) sf::st_linestring(m_xyz[m[, 3] == i, ]))
linestrings <- lapply(seq(n), function(i) sf::st_linestring(m_xyz[m[, "L1"] == i, ]))

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Would this be more specific and clear, seeing as we want a vector output?

Suggested change
linestrings <- lapply(seq(n), function(i) sf::st_linestring(m_xyz[m[, "L1"] == i, ]))
linestrings <- lapply(seq(n), function(i) sf::st_linestring(m_xyz[m[["L1"]] == i, ]))

Not tested but imagine it should work..

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, for a data.frame. but this is a matrix (from st_coordinates()), right? not sure if works the same way...

rgeom3d_sfc <- sf::st_sfc(linestrings, crs = sf::st_crs(routes))
sf::st_geometry(routes) <- rgeom3d_sfc
routes
Expand All @@ -287,3 +287,26 @@ has_terra <- function() requireNamespace("terra", quietly = TRUE)
is_linestring <- function(x) unique(sf::st_geometry_type(x)) == "LINESTRING"
stop_is_not_linestring <- function(x) if (!is_linestring(x)) stop("Only works with LINESTRINGs. Convert with sf::st_cast()")
stopifnotsf <- function(x, arg_name = "routes") if (!methods::is(x, "sf")) stop(arg_name, " is not an sf object. Try again with an sf object.")

#' Split a route into vertex-to-vertex segments
#'
#' Splits a linestring with XYZ coordinates into individual 2-point segments,
#' one per consecutive vertex pair. Useful for computing per-segment slopes
#' with [slope_xyz()].
#'
#' @param route_xyz An sf object with a single LINESTRING geometry with Z coordinates,
#' as returned by [elevation_add()].
#' @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)
#' summary(segs$slope)
route_to_segments <- function(route_xyz) {
coords <- sf::st_coordinates(route_xyz)
n <- nrow(coords)
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)))
}

26 changes: 26 additions & 0 deletions man/route_to_segments.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

57 changes: 34 additions & 23 deletions vignettes/slopes.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -98,15 +98,15 @@ library(sf)

# Load example data
data(lisbon_route)
dem_lisbon = dem_lisbon()
dem_lisbon <- dem_lisbon()

Copy link
Copy Markdown
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

these operator changes are annoying. please ignore.

```

## Add elevation to a linestring

If you have a 2D linestring and a DEM, you can add elevation data to the linestring using `elevation_add()`:

```{r}
sf_linestring_xyz_local = elevation_add(lisbon_route, dem = dem_lisbon)
sf_linestring_xyz_local <- elevation_add(lisbon_route, dem = dem_lisbon)
head(sf::st_coordinates(sf_linestring_xyz_local))
```

Expand All @@ -123,7 +123,7 @@ If you don't have a local DEM, `elevation_add()` can download elevation data (th
Once you have a 3D linestring (with XYZ coordinates), you can calculate its average slope using `slope_xyz()`:

```{r}
slope = slope_xyz(sf_linestring_xyz_local)
slope <- slope_xyz(sf_linestring_xyz_local)
slope
```

Expand All @@ -141,37 +141,48 @@ plot_slope(sf_linestring_xyz_local, pal = pal, brks = brks)

## Working with segments

The `slopes` package can also work with individual segments of a linestring.
First, let's segment the `lisbon_route`:
The `slopes` package can also work with individual segments of a linestring.
There are two ways to split a route into segments:

```{r}
lisbon_route_segments = sf::st_segmentize(lisbon_route, dfMaxLength = 100) # Arbitrary length
lisbon_route_segments = sf::st_cast(lisbon_route_segments, "LINESTRING")
# Add elevation to segments
lisbon_route_segments_xyz = elevation_add(lisbon_route_segments, dem = dem_lisbon)
```
### By vertex (native segments of the linestring)

Now calculate the slope for each segment:
`route_to_segments()` splits the route at every existing vertex, producing one 2-point segment per coordinate pair:

```{r}
lisbon_route_segments_xyz$slope = slope_xyz(lisbon_route_segments_xyz)
lisbon_route_xyz <- elevation_add(lisbon_route, dem = dem_lisbon())
lisbon_route_segments_xyz <- route_to_segments(lisbon_route_xyz)
lisbon_route_segments_xyz$slope <- slope_xyz(lisbon_route_segments_xyz)
summary(lisbon_route_segments_xyz$slope)
```

You can plot these segments, for example, colored by their slope. Here we use `tmap` for a more advanced plot (requires `tmap` package).

```{r, eval=FALSE}
# Requires tmap package
# library(tmap)
# qtm(lisbon_route_segments_xyz, lines.col = "slope", lines.lwd = 3)
```{r}
plot(st_geometry(lisbon_route_segments_xyz),
col = heat.colors(length(lisbon_route_segments_xyz$slope))[rank(lisbon_route_segments_xyz$slope)],
lwd = 3, main = "Slope by vertex segment"
)
```

Alternatively, using base R graphics:
### By fixed length (using stplanr)

`stplanr::line_segment()` splits the route into segments of a given length (e.g. 100 m).
Elevation must be added after segmenting, since the new endpoints won't have Z coordinates yet:

```{r}
plot(st_geometry(lisbon_route_segments_xyz), col = heat.colors(length(lisbon_route_segments_xyz$slope))[rank(lisbon_route_segments_xyz$slope)], lwd = 3)
if (requireNamespace("stplanr", quietly = TRUE)) {
lisbon_route_100m <- stplanr::line_segment(lisbon_route, segment_length = 100)
lisbon_route_100m_xyz <- elevation_add(lisbon_route_100m, dem = dem_lisbon())
lisbon_route_100m_xyz$slope <- slope_xyz(lisbon_route_100m_xyz)
summary(lisbon_route_100m_xyz$slope)
}
```

This vignette provides a basic overview. For more detailed information and advanced use cases, please refer to the other vignettes and the function documentation.

```{r}
if (requireNamespace("stplanr", quietly = TRUE)) {
plot(st_geometry(lisbon_route_100m_xyz),
col = heat.colors(length(lisbon_route_100m_xyz$slope))[rank(lisbon_route_100m_xyz$slope)],
lwd = 3, main = "Slope by 100 m segment"
)
}
```

This vignette provides a basic overview. For more detailed information and advanced use cases, please refer to the other vignettes and the function documentation.
Loading