我有多个国家/地区的多个点（经度、纬度）。这些点通过 ID 相关联，连接后会创建一条线。
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
我在这篇文章之后成功创建了几何图形 [1]，但是，我认为我的尝试可以进一步自动化。我尝试单独使用 group by ，如下所示：
# 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
我分享这一部分，以防有改进 group_by 部分的建议
我的实际问题是关于点之间的距离估计以及这些点之间的相对海拔。
对于距离估计，我已经尝试过：
aus_2 <- aus_1 %>% 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
(as.matrix(rbind(b[i, ], g[i, ]), by_element = TRUE)
)
),
crs = 4326
)
}))
Error in st_point(as.matrix(rbind(b[i, ], g[i, ]), by_element = TRUE)) :
nrow(x) == 1 is not TRUE
我认为该错误可能与矩阵创建有关。我很感谢这里的任何帮助。
对于相对高程估计，我还没有找到以前使用 R 进行的工作。欢迎提出关于从哪里开始进行高程估计的建议。
另外，1中的“c”是什么意思？它在函数中是如何工作的？
对于海拔，您需要一些数据源/提供者。
elevatr在这里，我假设“相对高程估计”是指高程范围，尽管对于线串来说它可以是多个指标。使用
sfst_length()group()summarise()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_sf <-
aus |>
# 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")
aus_sf
#> Simple feature collection with 17 features and 5 fields
#> Geometry type: POINT
#> Dimension: XY
#> Bounding box: xmin: 129.251 ymin: -20.06129 xmax: 137.2865 ymax: -18.30723
#> Geodetic CRS: WGS 84
#> First 10 features:
#> id country point_id geometry elevation elev_units
#> 1 1 Australia 0 POINT (130.1491 -19.5752) 382 meters
#> 2 1 Australia 1 POINT (129.9958 -19.4876) 425 meters
#> 3 1 Australia 2 POINT (129.7156 -19.25788) 444 meters
#> 4 1 Australia 3 POINT (129.7104 -19.20223) 441 meters
#> 5 2 Australia 0 POINT (129.251 -18.59016) 376 meters
#> 6 2 Australia 1 POINT (129.5436 -18.30723) 398 meters
#> 7 3 Australia 0 POINT (137.284 -20.06129) 229 meters
#> 8 3 Australia 1 POINT (137.2865 -20.04308) 234 meters
#> 9 3 Australia 2 POINT (137.1915 -20.00782) 237 meters
#> 10 3 Australia 3 POINT (137.122 -19.97166) 234 meters
从有序点构建线串：
aus_lines <-
aus_sf |>
# group by `country` to keep it as attribute and by `id`
# to create 4 multipoints, one per id
group_by(country, id) |>
# points to multipoints,
# do not use do_union as it will likely change point order,
# add elevation range
summarise(elev_range = max(elevation) - min(elevation), do_union = FALSE, .groups = "drop") |>
# multipoints to linestrings
st_cast("LINESTRING") |>
# add length column
mutate(length = st_length(geometry))
生成的具有高程范围和长度的线串：
aus_lines
#> Simple feature collection with 4 features and 4 fields
#> Geometry type: LINESTRING
#> Dimension: XY
#> Bounding box: xmin: 129.251 ymin: -20.06129 xmax: 137.2865 ymax: -18.30723
#> Geodetic CRS: WGS 84
#> # A tibble: 4 × 5
#> country id elev_range geometry length
#> * <chr> <int> <dbl> <LINESTRING [°]> [m]
#> 1 Australia 1 62 (130.1491 -19.5752, 129.9958 -19.4876, 129.… 63941.
#> 2 Australia 2 22 (129.251 -18.59016, 129.5436 -18.30723) 44072.
#> 3 Australia 3 8 (137.284 -20.06129, 137.2865 -20.04308, 137… 48462.
#> 4 Australia 4 4 (136.8961 -19.85932, 136.8791 -19.88669, 13… 10190.
数据集示例：
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")