测量 R 中一条直线上所有点的距离（sf 中的线串）

Question

我有一条 linestring 的平滑线（海岸线的简化抽象），我想沿着它每隔一段时间测量线的长度。我可以创建平滑线，并测量它的长度：

library(raster)
library(sf)
library(tidyverse)
library(rnaturalearth)
library(smoothr)

# create bounding box for the line 
xmin=-80
xmax=-66
ymin=24
ymax=45
bbox <- extent(xmin, xmax, ymin, ymax)

# get coarse outline of the US 
usamap <- rnaturalearth::ne_countries(scale = "small", country = "united states of america", returnclass = "sf")[1] %>% 
  st_cast("MULTILINESTRING")

# crop to extent of the bbox and get rid of a border line that isn't coastal 
bbox2 <- st_set_crs(st_as_sf(as(raster::extent(-80, -74, 42, 45.5), "SpatialPolygons")), st_crs(usamap))
cropmap <- usamap %>% 
  st_crop(bbox) %>% 
  st_difference(bbox2)

# smooth the line 
smoothmap <- cropmap %>% 
  smoothr::smooth(method="ksmooth", smoothness=8)

# measure the line length
st_length(smoothmap) # I get 1855956m

我将成为 "snapping" 沿着这条线的点的示例站点，我需要知道它们沿着海岸线有多远。但是，我不知道如何沿着海岸线每隔一段时间测量线的长度。间隔大小并不是非常重要，只要它具有相对精细的分辨率，也许每 1 公里或每 0.01 度。

我想要生成的数据框包含 x,y 列，其中包含沿线的 lat/lon 个点，"length" 列包含沿线的距离（从原点至此）。以下是我尝试过的一些方法：

迭代边界框。 我尝试用越来越小的边界框裁剪线（在 lat/lon 中使用规则间隔），但是因为线弯曲 "backward" 在 -76、38 左右，一些裁剪框没有包含我预期的完整线段。这种方法适用于该行的右上半部分，但不适用于左下半部分——只有 returns 长度为零。

在实现平滑器之前裁剪范围，然后测量线。因为平滑器函数不会产生相同的形状，如果它只测量原始线的一段，这实际上并不测量沿同一条线的距离。

用st_coordinates得到linestring的坐标，剪掉一行（线上一个点），将剩余的坐标重铸为linestring. 这种方法不会产生单个 linestring 而是产生一系列点（因为 st_cast 不知道如何再次连接它们），因此无法测量通常情况下。

"edit"smoothmap的几何最好一次删除一行坐标，重复测线，写出终点坐标和线长一个数据框。但是，我不确定是否可以在不将 sf 对象转换为数据框的情况下编辑它的坐标。

Answer 1

如果我理解你的问题，我想你可以这样做：

x <- as(smoothmap, "Spatial")
g <- geom(x)
d <- pointDistance(g[-nrow(g), c("x", "y")], g[-1, c("x", "y")], lonlat=TRUE)
gg <- data.frame(g[, c('x','y')], seglength=c(d, 0))

gg$lengthfromhere <- rev(cumsum(rev(gg[,"seglength"])))

head(gg)

#          x        y seglength lengthfromhere
#1 -67.06494 45.00000 70850.765        1855956
#2 -67.74832 44.58805  2405.180        1785105
#3 -67.77221 44.57474  2490.175        1782700
#4 -67.79692 44.56095  2577.254        1780210
#5 -67.82248 44.54667  2666.340        1777633
#6 -67.84890 44.53189  2757.336        1774967

tail(gg)
#            x        y  seglength lengthfromhere
#539 -79.09383 33.34224   2580.481       111543.5
#540 -79.11531 33.32753   2542.557       108963.0
#541 -79.13648 33.31306   2512.564       106420.5
#542 -79.15739 33.29874   2479.949       103907.9
#543 -79.17802 33.28460 101427.939       101427.9
#544 -80.00000 32.68751      0.000            0.0

Answer 2

我相信你需要sf::st_line_sample:

# Transform to metric
smoothmap_utm <- st_transform(smoothmap, 3857)

# Get samples at every kilometer
smoothmap_samples <- st_line_sample(smoothmap_utm, density = 1/1000)

# Transform back to a sf data.frame
smoothmaps_points <- map(smoothmap_samples, function(x) data.frame(geometry = st_geometry(x))) %>%
  map_df(as.data.frame) %>% st_sf() %>%
  st_cast("POINT") %>%
  st_set_crs(3857) %>%
  st_transform(4326)

mapview(smoothmaps_points) + mapview(smoothmap)

获得您想要的输出：

# Function to transform sf to lon,lat

sfc_as_cols <- function(x, names = c("lon","lat")) {
  stopifnot(inherits(x,"sf") && inherits(sf::st_geometry(x),"sfc_POINT"))
  ret <- sf::st_coordinates(x)
  ret <- tibble::as_tibble(ret)
  stopifnot(length(names) == ncol(ret))
  x <- x[ , !names(x) %in% names]
  ret <- setNames(ret,names)
  ui <- dplyr::bind_cols(x,ret)
  st_set_geometry(ui, NULL)
}

smoothmaps_points_xy <- sfc_as_cols(smoothmaps_points) %>%
  mutate(dist = cumsum(c(0, rep(1000, times = n() - 1))))

smoothmaps_points_xy

        lon      lat dist
1 -67.06836 44.99794    0
2 -67.07521 44.99383 1000
3 -67.08206 44.98972 2000
4 -67.08891 44.98560 3000
5 -67.09575 44.98149 4000
6 -67.10260 44.97737 5000

重要

但如果您的最终目标是获取路径中的点距离，我建议您检查 rgeos::gProject。

测量 R 中一条直线上所有点的距离（sf 中的线串）

Measuring the distance of all points along a line in R (linestring in sf)

r

spatial

sf

重要