Enhance elevation profile accuracy with cumulative path distance

[Aditya-Dahiya]

This commit adds an improved method for computing cumulative distances along a non-linear path when extracting elevation profiles from rasters. It replaces the straight-line (Euclidean) distance calculation with a stepwise sum of distances between consecutive points, ensuring more accurate results for curved transects. The new code enhances real-world applicability, especially for hiking trails and routes with turns or elevation changes.

Original Code from the Book: This computes the straight-line (Euclidean) distance from the first point to every other point using st_distance(). This works well only if the transect path is a straight line.

zion_transect <- zion_transect |> 
  group_by(id) |> 
  mutate(dist = st_distance(geometry)[, 1])

Proposed Improved Code: This code accounts for the real path distance, even when the route is not a straight line.

# This extracts coordinates of all points along the transect,
# grouped by ID to support multiple transects if present.

zion_transect_coordinates <- zion_transect |>  
  group_by(id) |>  
  st_coordinates()

# Calculates the cumulative distance along the actual (curved) path
# by computing Euclidean distance between each pair of consecutive points.

point_distances <- c(  
  0,   
  cumsum(sqrt(rowSums(diff(zion_transect_coordinates)^2)))  
)

# Adds the computed cumulative distances back to the main spatial object.

zion_transect <- zion_transect |>  
  mutate(dist = point_distances)

[Aditya-Dahiya]

Minor improvement in proposed code.

[Robinlovelace]

Minor changes suggested.

[Aditya-Dahiya]

Assignment operator changed from "<-" to "="

[Nowosad]

Thanks @Aditya-Dahiya -- I agree that the improvement is needed. One thing that still bothers me, though, is the generalization of the proposed approach. It will only work on projected CRSs -- if I understand correctly, the outcome for the data in geographical CRS would be in degrees (not meters), which is not good.

See an example and let me know what you think.

library(sf)
library(terra)
library(dplyr)

zion_transect = cbind(c(-113.2, -112.9), c(37.45, 37.2)) |>
  st_linestring() |> 
  st_sfc(crs = "EPSG:4326") |>
  st_sf(geometry = _)

zion_transect$id = 1:nrow(zion_transect)
zion_transect = st_segmentize(zion_transect, dfMaxLength = 250)
zion_transect = st_cast(zion_transect, "POINT")

zion_transect_coordinates = zion_transect |> group_by(id) |> st_coordinates()

# ~true
zion_tc = zion_transect_coordinates[c(1, 240), ]
zion_tc
#>              X        Y
#> [1,] -113.2000 37.45000
#> [2,] -112.9199 37.21662

zion_tc_sf = zion_tc |>
  as.data.frame() |>
  st_as_sf(coords = c("X", "Y"), crs = "EPSG:4326")


st_distance(zion_tc_sf)
#> Units: [m]
#>          [,1]     [,2]
#> [1,]     0.00 35872.97
#> [2,] 35872.97     0.00

# incorrect?
point_distances0 = c(
  0, 
  cumsum(sqrt(rowSums(diff(zion_tc)^2)))
)
point_distances0 # this result is in "degrees"
#> [1] 0.0000000 0.3646125