Feature suggestion: add get_profile_elev() function for returning elevation along a line
#93
Comments
|
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
<#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>
|
|
@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! |
|
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 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 Regardless, take your time on the review – nothing urgent on my end! |
Thanks for the handy package! I was wondering if you may be interested in a
get_profile_elev()function that usessf::st_line_sample()to get points along a line before callingelevatr::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.Created on 2023-11-03 with reprex v2.0.2
The text was updated successfully, but these errors were encountered: