I have multiple points (lon, lat) for a range of countries. These points are related by an ID that when connected creates a line.
data exmaple:
aus:
id country point_id lon lat
1 Australia 0 130.1491 -19.57520
1 Australia 1 129.9958 -19.48760
1 Australia 2 129.7156 -19.25788
1 Australia 3 129.7104 -19.20223
2 Australia 0 129.2510 -18.59016
2 Australia 1 129.5436 -18.30723
3 Australia 0 137.2840 -20.06129
3 Australia 1 137.2865 -20.04308
3 Australia 2 137.1915 -20.00782
3 Australia 3 137.1220 -19.97166
3 Australia 4 137.0650 -19.91363
3 Australia 5 136.8961 -19.85932
4 Australia 0 136.8961 -19.85932
4 Australia 1 136.8791 -19.88669
4 Australia 2 136.8594 -19.91227
4 Australia 3 136.8454 -19.92507
4 Australia 4 136.8360 -19.92976
I managed to create the geometry following this post [1], however, I think my attempt can be automated further. I attempted using group by individually, as follow:
attempt 1 to get linestrings:
# same logic as in [1] but by group
aus_group <- aus %>% group_by(id, country)
b_group <- aus %>% group_by(id, country) %>% select(lon, lat)
e_group <- aus %>% group_by(id, country) %>% filter(row_number()!=1) %>% select(lon, lat)
f_group <- e_group %>% group_by(id, country) %>% summarise_all(last) %>% select(country, lon, lat)
g_group <- e_group %>% group_by(id, country) %>% rbind(f_group) %>% arrange(id)
aus_group$geometry = do.call(
"c",
lapply(seq(nrow(b_group)), function(i) {
st_sfc(
st_linestring(
as.matrix(
rbind(b_group[i, c("lon", "lat")], g_group[i, c("lon", "lat")])
)
),
crs = 4326
)
}))
dat_g_sf = st_as_sf(aus_group)
mapview(dat_g_sf, zcol = "id")
# to compare with approach in [1]
dat_g_sf %>% filter(id==1) %>% mapview(zcol = "point_id") # got same plot
# so it works
I share this part in case there are suggestions to improve the group_by part
question:
My actual question is regarding distance estimation between points and the relative elevation between those points.
For the distance estimation, I have tried this:
attempt 2 to get distances:
# edited aus_1 becomes aus
aus_2 <- aus %>% select(id, country, point_id, lon, lat) #to avoid error with geometry
str(aus_2)
aus_2$distance = do.call(
"c",
lapply(seq(nrow(b)), function(i) {
st_distance(
st_sfc(st_point
# edited b becomes b_group, and g becomes g_group
(as.matrix(rbind(b_group[i, c("lon", "lat")], g_group[i, c("lon", "lat")], by_element = TRUE)
)
),
crs = 4326
)
}))
Error in st_point(as.matrix(rbind(b_group[i, c("lon", "lat")], g_group[i, c("lon", "lat")]), by_element = TRUE)) :
nrow(x) == 1 is not TRUE
I think the error might be related with the matrix creation. I appreciate any help here.
expected results: edited.
As can be seen in my attempt 2, I tried to get the length between points (by point_id), e.g., distance between point_id: 0 & 1, in id: 1; and so on; using st_distance. (I guess the sum of all distances by point_id would result on the total length estimation).
In my understanding, the distance for point_id == 0 should be zero because it is the distance between point_id: 0 and 0, while the distance between point_id: 0 and 1 is the actual distance for point_id: 1. This is the reason why in my attempt 1 and 2, I tried to create two groups of c("lon", "lat"): b_group and g_group, and then I try applying st_distance between point b and point g. Unfortunately, this did not work.
Overall: it is expected "an sf object of point geometries with the same number of rows as an input dataset, with distance to lagging or leading point (between point_id) within the same group (id) + elevation (between by point_id)" from margusl comment.
no edited:
For the relative elevation estimation, I haven’t found previous work on this using R. Suggestions regarding where to start about elevation estimation are welcome.
Also, what is the meaning of "c" in 1? How does it work in the function?
For elevation you'd need some data source / provider. elevatr package potentially can provide you values for point locations, for example through Amazon Web Service Terrain Tiles, but you'd probably need to evaluate those results and work out acceptable resolution / zoom levels, elevatr docs should cover this part.
To get distances between every current and previous point, you can use element-wise distance with lag:
library(elevatr)
library(dplyr)
library(sf)
#> Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
aus |>
# tibble gets us more compact printout
tibble() |>
# make sure points are correctyl sorted
arrange(country, id, point_id) |>
# convert to sf object
st_as_sf(coords = c("lon", "lat"), crs = "WGS84") |>
# get elevation from Amazon Web Service Terrain Tiles
get_elev_point(src = "aws") |>
# distance from previuos point to the current one, within the same id group
mutate(dist = st_distance(geometry, lag(geometry), by_element = TRUE),
# or switch from lag() to lead() to get distance to the next point
# dist_to_next = st_distance(geometry, lead(geometry), by_element = TRUE),
.by = id)
#> Mosaicing & Projecting
#> Note: Elevation units are in meters
#> Simple feature collection with 17 features and 6 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 129.251 ymin: -20.06129 xmax: 137.2865 ymax: -18.30723
#> Geodetic CRS: WGS 84
#> # A tibble: 17 × 7
#> id country point_id geometry elevation elev_units dist
#> * <int> <chr> <int> <POINT [°]> <dbl> <chr> [m]
#> 1 1 Austral… 0 (130.1491 -19.5752) 382 meters NA
#> 2 1 Austral… 1 (129.9958 -19.4876) 425 meters 18788.
#> 3 1 Austral… 2 (129.7156 -19.25788) 444 meters 38941.
#> 4 1 Austral… 3 (129.7104 -19.20223) 441 meters 6212.
#> 5 2 Austral… 0 (129.251 -18.59016) 376 meters NA
#> 6 2 Austral… 1 (129.5436 -18.30723) 398 meters 44072.
#> 7 3 Austral… 0 (137.284 -20.06129) 229 meters NA
#> 8 3 Austral… 1 (137.2865 -20.04308) 234 meters 2042.
#> 9 3 Austral… 2 (137.1915 -20.00782) 237 meters 10671.
#> 10 3 Austral… 3 (137.122 -19.97166) 234 meters 8301.
#> 11 3 Austral… 4 (137.065 -19.91363) 235 meters 8783.
#> 12 3 Austral… 5 (136.8961 -19.85932) 235 meters 18665.
#> 13 4 Austral… 0 (136.8961 -19.85932) 235 meters NA
#> 14 4 Austral… 1 (136.8791 -19.88669) 235 meters 3525.
#> 15 4 Austral… 2 (136.8594 -19.91227) 233 meters 3512.
#> 16 4 Austral… 3 (136.8454 -19.92507) 233 meters 2042.
#> 17 4 Austral… 4 (136.836 -19.92976) 237 meters 1112.
Dataset example:
aus <- read.table(header = TRUE, text =
"id country point_id lon lat
1 Australia 0 130.1491 -19.57520
1 Australia 1 129.9958 -19.48760
1 Australia 2 129.7156 -19.25788
1 Australia 3 129.7104 -19.20223
2 Australia 0 129.2510 -18.59016
2 Australia 1 129.5436 -18.30723
3 Australia 0 137.2840 -20.06129
3 Australia 1 137.2865 -20.04308
3 Australia 2 137.1915 -20.00782
3 Australia 3 137.1220 -19.97166
3 Australia 4 137.0650 -19.91363
3 Australia 5 136.8961 -19.85932
4 Australia 0 136.8961 -19.85932
4 Australia 1 136.8791 -19.88669
4 Australia 2 136.8594 -19.91227
4 Australia 3 136.8454 -19.92507
4 Australia 4 136.8360 -19.92976")
Original answer - https://stackoverflow.com/revisions/78251444/1