Feature suggestion: add `get_profile_elev()` function for returning elevation along a line

[elipousson]

Thanks for the handy package! I was wondering if you may be interested in a get_profile_elev() function that uses sf::st_line_sample() to get points along a line before calling elevatr::get_elev_point() and optionally get distance between points on the line. I put together a prototype that I think works pretty well – let me know if you're interested and I'd be happy to open a pull request.

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}

[pnogas67]

I think this is a very interesting feature. Thx for sharing, PN

On Fri, Nov 3, 2023 at 4:19 AM Eli Pousson @.***> wrote:

Thanks for the handy package! I was wondering if you may be interested in a get_profile_elev() function that uses sf::st_line_sample() to get points along a line before calling elevatr::get_elev_point() and optionally get distance between points on the line. I put together a prototype that I think works pretty well – let me know if you're interested and I'd be happy to open a pull request.

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.)" )

https://camo.githubusercontent.com/dee70e678d68ca07c4ff57336364b849f878c796d639ef2d938a526b7723aaf4/68747470733a2f2f692e696d6775722e636f6d2f6b65686958486b2e706e67

Created on 2023-11-03 with reprex v2.0.2 https://reprex.tidyverse.org

— Reply to this email directly, view it on GitHub https://github.com/jhollist/elevatr/issues/93, or unsubscribe https://github.com/notifications/unsubscribe-auth/AB7E6RMOXTWLHRHLHNRZPZLYCRWGNAVCNFSM6AAAAAA633V7RCVHI2DSMVQWIX3LMV43ASLTON2WKOZRHE3TKNBTGI4DSNA . You are receiving this because you are subscribed to this thread.Message ID: @.***>

--

Pedro Nogueira (GMAIL)

Home: http://evunix.uevora.pt/~pmn Ph: +351 266 745301 Economic Geology (PhD) Departamento de Geociências da Universidade de Évora

http://goog_46339544

[jhollist]

@elipousson So sorry for not responding to your initial ping about this. I am a one man shop on elevatr and my time is short to spend on it.

That being said, this is a really nice idea and thanks for submitting the PR. I am out of the office for the next 10 days or so, but I'll take a look at the PR when I get back and will comment then.

I am not opposed at all to merging this I just want to think about it some. It might make more sense as a stand alone package of elevation tools like this, but then again maybe not. Need to mull!

[elipousson]

No apologies needed. Totally appreciate your work as a volunteer on this project.

I can imagine making the case for either approach. A third option could be to export the helper functions that are making the draft version of get_elev_profile() work and expand the documentation to show how you can get this result without the need to maintain a dedicated function.

FWIW – the PR has a couple of other changes, particularly the option to customize the coordinate column names when using an input data frame and the option to customize the output units for get_elev_point() that I think could be really helpful beyond this specific issue. If you wanted me to open some separate issues for those changes for good open-source housekeeping, I'd be happy to do so.

Regardless, take your time on the review – nothing urgent on my end!