New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
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>
|
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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: