如何从一个点画一条线到多边形边缘，然后在 R 中以 sf 形式获得线的长度？

Question

还有一些其他帖子相关与此相关，例如：Post 1、Post 2、Post 3。但是，它们都没有提供我所希望的。我想要的是能够在特定方向（“正南”又名向下）从特定点（采样位置）到完全围绕该点（湖泊边界）的多边形边缘绘制一条线段。然后我想测量采样点和多边形边缘之间的线段的长度（真的，这只是我想要的距离，所以如果我们可以在不绘制线段的情况下获得距离，那就更好了！）。不幸的是，sf 包中似乎不存在执行此操作的功能：See closed issue here。

我怀疑，不过，这可以通过修改此处提供的解决方案来实现：See copy-pasted code below, modified by me。但是，我对sf 中的工具非常不满意--我只制作了从点本身到多边形南部范围的线段，在某个点与多边形相交：

library(sf)
library(dplyr)

df = data.frame(
  lon = c(119.4, 119.4, 119.4, 119.5, 119.5),
  lat = c(-5.192,-5.192,-5.167,-5.167,-5.191)
)

polygon <- df %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("POLYGON")

plot(polygon)

df2 <- data.frame(lon = c(119.45, 119.49, 119.47),
                  lat = c(-5.172,-5.190,-5.183))

points <- df2 %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326) %>%
  summarise(geometry = st_combine(geometry)) %>%
  st_cast("MULTIPOINT")

plot(points, add = TRUE, col = "red")

# Solution via a loop

xmin <- min(df$lat)

m = list()
# Iterate and create lines
for (i in 1:3) {
  m[[i]] = st_linestring(matrix(
    c(df2[i, "lon"],
      df2[i, "lat"],
      df2[i, "lon"],
      xmin),
    nrow = 2,
    byrow = TRUE
  ))
}
test = st_multilinestring(m)

# Result is line MULTILINESTRING object

plot(test, col = "green", add = TRUE)

但现在我不知道如何使用st_intersection 或任何此类函数来确定交点的位置。我认为，大部分问题在于我创建的不是sf 对象，我不知道如何让它成为一个对象。我假设，如果我能找出线段与多边形相交的位置（或者理想情况下它们相交的最北时间），我可以使用类似st_distance 的函数以某种方式从交点到采样点进行测量。但是，由于湖泊多边形通常真的很复杂，因此一个线段可能会多次与多边形相交（例如，如果在给定点以南有一个半岛），在这种情况下，我想我可以找到每个采样点的“最北”交点并使用该交点，否则为每个采样点取最小的此类距离。

无论如何，如果有人可以告诉我我缺少的几个步骤，那就太好了！我觉得我很近，但又很远......

Answer 1

考虑一下这种方法，大致灵感来自我之前关于 lines from points 的帖子

为了使其更具可重复性，我使用了众所周知且深受喜爱的北卡罗来纳州 shapefile，该文件随 {sf} 和三个半随机 NC 城市的数据框一起提供。

代码的作用是：

• 通过 for 循环遍历城市的数据框
• 创建一条从每个城市开始（“观察”）并在南极结束的线
• 与解散的北卡罗来纳州相交
• 将交叉点炸成单独的线串
• 选择在原点 1 米范围内通过的线串
• 通过sf::st_lenghth()计算长度
• 将结果保存为名为@9​​87654328@（结果的缩写：）的{sf}数据框

我在最终对象中包含了实际的行以使结果更清晰，但您可以选择省略它。

library(sf)
library(dplyr)
library(ggplot2)

shape <- st_read(system.file("shape/nc.shp", package="sf")) %>%  # included with sf package
  summarise() %>% 
  st_transform(4326) # to align CRS with cities

cities <- data.frame(name = c("Raleigh", "Greensboro", "Plymouth"),
                     x = c(-78.633333, -79.819444, -76.747778),
                     y = c(35.766667, 36.08, 35.859722)) %>% 
  st_as_sf(coords = c("x", "y"), crs = 4326)

# a quick overview
ggplot() +
   geom_sf(data = shape) + # polygon of North Carolina
   geom_sf(data = cities, color = "red") # 3 cities  

# now here's the action!!!
for (i in seq_along(cities$name)) {
  
  # create a working linestring object
  wrk_line <- st_coordinates(cities[i, ]) %>% 
    rbind(c(0, -90)) %>% 
    st_linestring() %>% 
    st_sfc(crs = 4326) %>% 
    st_intersection(shape) %>% 
    st_cast("LINESTRING") # separate individual segments of multilines
  
  first_segment <- unlist(st_is_within_distance(cities[i, ], wrk_line, dist = 1))  
  
  # a single observation
  line_data <- data.frame(
    name = cities$name[i],
    length = st_length(wrk_line[first_segment]),
    geometry = wrk_line[first_segment]
  )
  
  # bind results rows to a single object
  if (i == 1) {
    res <- line_data
    
  } else {
    res <- dplyr::bind_rows(res, line_data)
    
  } # /if - saving results
  
  
} # /for

# finalize results
res <- sf::st_as_sf(res, crs = 4326) 
  
# result object    
res
# Simple feature collection with 3 features and 2 fields
# Geometry type: LINESTRING
# Dimension:     XY
# Bounding box:  xmin: -79.81944 ymin: 33.92945 xmax: -76.74778 ymax: 36.08
# Geodetic CRS:  WGS 84
#         name        length                       geometry
# 1    Raleigh 204289.21 [m] LINESTRING (-78.63333 35.76...
# 2 Greensboro 141552.67 [m] LINESTRING (-79.81944 36.08...
# 3   Plymouth  48114.32 [m] LINESTRING (-76.74778 35.85...

# a quick overview of the lines
ggplot() +
  geom_sf(data = shape) + # polygon of North Carolina
  geom_sf(data = res, color = "red") # 3 lines