如何从点到多边形边缘画一条线，然后在 R 中以 sf 获取线的长度？
How can I draw a line from a point to a polygon edge and then get the line's length in sf in R?
还有一些其他帖子与此相关，例如：Post 1, Post 2, . However, none of them deliver what I am hoping for. What I want is to be able to draw a line segment from a specific point (a sampling location) to the edge of a polygon fully surrounding that point (a lake border) in a specific direction ("due south" aka downward). I then want to measure the length of that line segment in between the sampling point and the polygon edge (really, it's only the distance I want, so if we can get the distance without drawing the line segment, so much the better!). Unfortunately, it doesn't seem like functionality to do this already exists within the sf package: 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 之类的函数以某种方式测量从交点到采样点的距离.但是，由于湖泊多边形通常 确实 复杂，因此一个线段可能会多次与多边形相交（例如给定点以南有一个半岛），在这种情况下我认为我可以为每个采样点找到“最北端”的交点并使用它，或者为每个采样点取最小这样的距离。
无论如何，如果有人能告诉我我遗漏的几个步骤，那就太好了！我觉得我很近又很远...
考虑一下这种方法，大致上受到我之前 post 关于 lines from points
的启发
为了使其更具可重现性，我使用了众所周知且备受喜爱的北卡罗来纳州 shapefile，它随 {sf} 和三个 semi-random NC 城市的数据框一起提供。
代码的作用是：
- 循环遍历城市数据框
- 创建一条线，从每个城市（“观测点”）开始，到南极结束
- 与已解散的北卡罗来纳州相交
- 将交叉点爆破为单独的线串
- 选择原点 1 米以内的线串
- 通过
sf::st_lenghth() 计算长度
- 将结果保存为名为
res（结果的缩写 :) 的 {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
还有一些其他帖子与此相关，例如：Post 1, Post 2, sf package: 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 之类的函数以某种方式测量从交点到采样点的距离.但是，由于湖泊多边形通常 确实 复杂，因此一个线段可能会多次与多边形相交（例如给定点以南有一个半岛），在这种情况下我认为我可以为每个采样点找到“最北端”的交点并使用它，或者为每个采样点取最小这样的距离。
无论如何，如果有人能告诉我我遗漏的几个步骤，那就太好了！我觉得我很近又很远...
考虑一下这种方法，大致上受到我之前 post 关于 lines from points
的启发为了使其更具可重现性，我使用了众所周知且备受喜爱的北卡罗来纳州 shapefile，它随 {sf} 和三个 semi-random NC 城市的数据框一起提供。
代码的作用是：
- 循环遍历城市数据框
- 创建一条线，从每个城市（“观测点”）开始，到南极结束
- 与已解散的北卡罗来纳州相交
- 将交叉点爆破为单独的线串
- 选择原点 1 米以内的线串
- 通过
sf::st_lenghth()计算长度
- 将结果保存为名为
res（结果的缩写 :) 的 {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