Consider adding a get_profile_elev() and get_location_elev() function wrapping {elevatr} package

[elipousson]

I have an open GitHub issue in the {elevatr} package proposing the addition of this function but I may add it to getdata regardless:

library(sf)
#> Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(elevatr)
#> elevatr v0.99.0 NOTE: Version 0.99.0 of 'elevatr' uses 'sf' and 'terra'.  Use 
#> of the 'sp', 'raster', and underlying 'rgdal' packages by 'elevatr' is being 
#> deprecated; however, get_elev_raster continues to return a RasterLayer.  This 
#> will be dropped in future versions, so please plan accordingly.
library(ggplot2)

get_profile_elev <- function(locations,
                             units = NULL,
                             dist = FALSE,
                             n = NULL,
                             density = NULL,
                             type = "regular",
                             sample = NULL,
                             ...,
                             keep_units = FALSE,
                             cumulative = FALSE) {
  if (!inherits(locations, "sf")) {
    locations <- sf::st_as_sf(locations)
  }

  if (is.numeric(sample) || is.numeric(n) || is.numeric(density)) {
    locations <- sf::st_line_sample(
      locations,
      n = n, density = density,
      sample = sample, type = type
    )

    locations <- sf::st_cast(locations, to = "POINT")
  } else if (!sf::st_is(locations, "POINT")) {
    locations <- suppressWarnings(sf::st_cast(locations, to = "POINT"))
  }

  location_coords <- as.data.frame(sf::st_coordinates(locations))
  location_coords <- setNames(location_coords, tolower(names(location_coords)))
  locations_crs <- sf::st_crs(locations)

  elev_point <- elevatr::get_elev_point(location_coords, prj = locations_crs)

  if (nrow(elev_point) < 2) {
    dist <- FALSE
  }

  if (dist) {
    dist_point <- 0
    elev_point_geom <- vctrs::vec_chop(sf::st_geometry(elev_point))

    for (i in (seq(length(elev_point_geom) - 1))) {
      dist_add <- sf::st_distance(
        elev_point_geom[[i]],
        elev_point_geom[[i + 1]]
      )

      dist_point <- c(dist_point, dist_add)
    }

    if (cumulative) {
      dist_point <- cumsum(dist_point)
    }

    elev_point[["distance"]] <- units::set_units(
      dist_point, locations_crs$units_gdal,
      mode = "standard"
    )
  }

  elev_units <- unique(elev_point[["elev_units"]])

  if (!is.null(units)) {
    mode <- "standard"
    elev_units_label <- units

    if (!is.character(units)) {
      elev_units_label <- unique(as.character(base::units(x)[["numerator"]]))
      mode <- units::units_options("set_units_mode")
    }

    elev_point[["elevation"]] <- units::as_units(
      elev_point[["elevation"]],
      elev_units
    )

    elev_point[["elevation"]] <- units::set_units(
      elev_point[["elevation"]],
      value = units, mode = mode
    )

    elev_point[["elev_units"]] <- elev_units_label

    if (dist) {
      elev_point[["distance"]] <- units::set_units(
        elev_point[["distance"]],
        value = units, mode = mode
      )
    }
  } else if (dist) {
    elev_point[["distance"]] <- units::set_units(
      elev_point[["distance"]], elev_units
    )
  }

  if (!keep_units) {
    elev_point[["elevation"]] <- units::drop_units(elev_point[["elevation"]])

    if (dist) {
      elev_point[["distance"]] <- units::drop_units(elev_point[["distance"]])
    }
  }

  elev_point
}

nc <- st_read(system.file("shape/nc.shp", package = "sf")) |>
  st_transform(3857)
#> Reading layer `nc' from data source 
#>   `/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/sf/shape/nc.shp' 
#>   using driver `ESRI Shapefile'
#> Simple feature collection with 100 features and 14 fields
#> Geometry type: MULTIPOLYGON
#> Dimension:     XY
#> Bounding box:  xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
#> Geodetic CRS:  NAD27

nc_line <- suppressWarnings(
  st_cast(
    st_union(
      st_centroid(nc[1, ]),
      st_centroid(nc[2, ])
    ),
    to = "LINESTRING"
  )
)

elev_profile <- get_profile_elev(
  nc_line,
  units = "ft",
  dist = TRUE,
  cumulative = TRUE,
  n = 15
)
#> Downloading point elevations:
#> Note: Elevation units are in meters

elev_profile |>
  ggplot(aes(distance, elevation)) +
  geom_area() +
  scale_x_continuous(labels = scales::label_number()) +
  labs(
    x = "Distance (ft.)",
    y = "Elevation (ft.)"
  )

^{Created on 2023-11-03 with reprex v2.0.2}