由于您正在使用meuse数据,因此使用空间对象和空间计算似乎是很自然的:
library(sf)
library(sp)
data(meuse)
newline = cbind(x = seq(from = 178000, to = 181000, length.out = 100),
y = seq(from = 330000, to = 334000, length.out = 100))
meuse <- st_as_sf(meuse, coords = c("x", "y"))
newline <- st_linestring(newline)
st_distance(meuse, newline) [1:10]
plot_sf(meuse)
plot(meuse, add = TRUE)
plot(newline, add = TRUE)
您可以通过在仅有2个坐标的直线上运行相同的代码来证明这些是相对于您的直线的垂直距离。
但请注意,这是到直线的最短距离。因此,对于靠近线段端点或超出直线范围的点,您将无法得到垂直距离(只有到端点的最短距离)。
您必须有足够长的线以避免这种情况...
newline = cbind(x = c(178000, 181000),
y = c(330000, 334000))
meuse <- st_as_sf(meuse, coords = c("x", "y"), crs = 31370)
newline <- st_linestring(newline)
st_distance(meuse, newline) [1:10]
plot_sf(meuse)
plot(meuse, add = TRUE)
plot(newline, add = TRUE)
如果直线的斜率为1,则可以使用勾股定理(x的差异= y的差异=到直线的距离/ sqrt(2))计算点在直线上的投影。但是这里不适用(红色线段与直线不垂直),因为斜率不是1,因此y坐标的差异不等于x坐标的差异。勾股定理在这里无法解决。(但是由st_distance计算出的距离是垂直于直线的。)
distances <- st_distance(meuse, newline)
newline = data.frame(x = c(178000, 181000),
y = c(330000, 334000))
segments <- as.data.frame(st_coordinates(meuse))
segments <- data.frame(
segments,
X2 = segments$X - distances/sqrt(2),
Y2 = segments$Y + distances/sqrt(2)
)
library(ggplot2)
ggplot() +
geom_point(data = segments, aes(X,Y)) +
geom_line(data = newline, aes(x,y)) +
geom_segment(data = segments, aes(x = X, y = Y, xend = X2, yend = Y2),
color = "orangered", alpha = 0.5) +
coord_equal() + xlim(c(177000, 182000))
![enter image description here](https://istack.dev59.com/6qWgX.webp)
如果您只想绘制它,可以使用
rgeos
的
gProject
函数(目前在sf中没有等价函数)来获取点在线上投影的坐标。 您需要将点和线以s
sp
格式而不是
sf
格式提供,并在 sp 的矩阵和 ggplot 的数据框之间进行转换。
library(sp)
library(rgeos)
newline = cbind(x = c(178000, 181000),
y = c(330000, 334000))
spline <- as(st_as_sfc(st_as_text(st_linestring(newline))), "Spatial")
position <- gProject(spline, as(meuse, "Spatial"))
position <- coordinates(gInterpolate(spline, position))
colnames(position) <- c("X2", "Y2")
segments <- data.frame(st_coordinates(meuse), position)
library(ggplot2)
ggplot() +
geom_point(data = segments, aes(X,Y)) +
geom_line(data = as.data.frame(newline), aes(x,y)) +
geom_segment(data = segments, aes(x = X, y = Y, xend = X2, yend = Y2),
color = "orangered", alpha = 0.5) +
coord_equal() + xlim(c(177000, 182000))
![enter image description here](https://istack.dev59.com/wkCzI.webp)