使用R sf包计算多国家多几何点间距离及相对高程的技术问询

Hey there! Let's work through your problems step by step—first fixing that distance calculation, then simplifying your linestring code, and finally tackling elevation estimation.

1. 简化线串（Linestring）的创建 #

Your current approach works, but we can make it much cleaner using dplyr's grouping with sf functions directly:

library(sf)
library(dplyr)

# 按id和country分组，将每组的点转换为线串
aus_linestrings <- aus %>%
  group_by(id, country) %>%
  arrange(point_id) %>%  # 确保点按point_id顺序排列
  summarise(
    geometry = st_sfc(st_linestring(as.matrix(select(., lon, lat)))),
    .groups = "drop"
  ) %>%
  st_set_crs(4326)

mapview(aus_linestrings, zcol = "id")

This way, you don't need to manually split into multiple groups—summarise handles creating one linestring per group automatically.

2. 计算连续点之间的距离 #

The error in your attempt comes because st_point expects a single point (1 row of coordinates), but you were passing two points. Instead, we can create lagged coordinates for each group, then calculate the distance between each point and its predecessor:

aus_with_distances <- aus %>%
  group_by(id, country) %>%
  arrange(point_id) %>%
  # 创建前一个点的经纬度列
  mutate(
    lon_prev = lag(lon),
    lat_prev = lag(lat),
    # 计算当前点与前一个点的距离（第一个点距离为0）
    distance = ifelse(is.na(lon_prev), 0, 
                      st_distance(
                        st_sfc(st_point(c(lon, lat)), crs = 4326),
                        st_sfc(st_point(c(lon_prev, lat_prev)), crs = 4326),
                        by_element = TRUE
                      ))
  ) %>%
  ungroup() %>%
  # 移除临时列（可选）
  select(-lon_prev, -lat_prev)

# 查看结果
head(aus_with_distances)

This gives you a distance value for every row: 0 for the first point in each group, and the actual distance from the previous point for all others. If you want the total length of each line, you can group by id and sum the distance column.

3. 相对高程估算 #

For elevation data in R, the elevatr package is a great starting point—it pulls elevation data from various sources (like AWS Terrain Tiles). Here's how to get elevation for your points and calculate relative elevation differences:

First, install and load the package:

install.packages("elevatr")
library(elevatr)

Then, add elevation to your dataset and compute relative differences:

# 将数据转换为sf点对象
aus_sf <- aus %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

# 获取高程数据（单位：米）
aus_sf <- get_elev_point(aus_sf, prj = "+proj=longlat +datum=WGS84", src = "aws")

# 计算连续点之间的高程差
aus_sf <- aus_sf %>%
  group_by(id, country) %>%
  arrange(point_id) %>%
  mutate(
    elevation_prev = lag(elevation),
    relative_elevation = elevation - elevation_prev
  ) %>%
  ungroup()

# 查看结果
head(aus_sf)

Note: The relative_elevation will be NA for the first point in each group—you can set that to 0 if needed, similar to the distance calculation.

4. 解释do.call中的"c" #

In the post you referenced, do.call("c", ...) is used to concatenate a list of sfc (simple feature geometry column) objects into a single sfc vector.

When you use lapply to create a list of individual st_sfc elements (each containing one linestring), do.call("c", list) takes all those elements and combines them into one continuous sfc column—this is necessary because an sf object's geometry column needs to be a single sfc vector, not a list of separate ones. It's equivalent to using c(list[[1]], list[[2]], ...) but works dynamically for any number of elements.

备注：内容来源于stack exchange，提问作者dmoyaec