交点线与轮廓线的循环坐标- r

Question

基于这个答案如何得到一个轮廓线相交的坐标- R，我尝试使用下面的脚本运行一个循环。知道为什么我不能画出所有的交点和直线吗?形状与给出的答案不同

代码:

library(ggplot2)
library(sf)

t <- seq(0, 2*pi, by=0.1)
df <- data.frame(x = 13*sin(t)^3,
y = 4*cos(t)-2*cos(3*t)-5*cos(4*t)-cos(2*t))
df <- rbind(df, df[1,]) # close the polygon

meanX <- mean(df$x)
meanY <- mean(df$y)

# Transform your data.frame in a sf polygon (the first and last points
# must have the same coordinates)

#> Linking to GEOS 3.5.1, GDAL 2.1.3, proj.4 4.9.2
poly <- st_sf(st_sfc(st_polygon(list(as.matrix(df)))))
# Choose the angle (in degrees)
rotAngles <- 5
for(angle in seq(0,359,rotAngles)) {

# Find the minimum length for the line segment to be always 
# outside the cloud whatever the choosen angle
maxX <- max(abs(abs(df[,"x"]) - abs(meanX)))
maxY <- max(abs(abs(df[,"y"]) - abs(meanY)))
line_length = sqrt(maxX^2 + maxY^2) + 1
# Find the coordinates of the 2 points to draw a line with 
# the intended angle.
# This is the gray line on the graph below
line <- rbind(c(meanX,meanY), 
c(meanX + line_length * cos((pi/180)*angle), 
meanY + line_length * sin((pi/180)*angle)))
# Transform into a sf line object
line <- st_sf(st_sfc(st_linestring(line)))
# Intersect the polygon and line. The result is a two points line
# shown in black on the plot below
intersect_line <- st_intersection(poly, line)
# Extract only the second point of this line.
# This is the intersecting point
intersect_point <- st_coordinates(intersect_line)[2,c("X","Y")] 
# Visualise this with ggplot and without geom_sf
# you need first transform back the lines into data.frame
line <- as.data.frame(st_coordinates(line))[,1:2]
intersect_line <- as.data.frame(st_coordinates(intersect_line))[,1:2]

ggplot() + geom_path(data=df, aes(x = x, y = y)) + 
geom_line(data=line, aes(x = X, y = Y), color = "gray80", lwd = 3) + 
geom_line(data=intersect_line, aes(x = X, y = Y), color = "gray20", lwd = 1) + 
geom_point(aes(meanX, meanY), colour="orangered", size=2) +
geom_point(aes(intersect_point["X"], intersect_point["Y"]), 
colour="orangered", size=2) +
theme_bw()
}

Answer 1

首先，我们将回到@Gilles多边形形状，因为它更符合他的推理和演示:

# Generate a heart shape
t <- seq(0, 2*pi, by=0.1)
df <- data.frame(x = 16*sin(t)^3, 
y = 13*cos(t)-5*cos(2*t)-2*cos(3*t)-cos(4*t))
df <- rbind(df, df[1,]) # close the polygon
meanX <- mean(df$x)
meanY <- mean(df$y)
library(sf)
poly <- st_sf(st_sfc(st_polygon(list(as.matrix(df)))))

这些元素不会改变，也不需要在循环中多次计算:

maxX <- max(abs(abs(df[,"x"]) - abs(meanX)))
maxY <- max(abs(abs(df[,"y"]) - abs(meanY)))
line_length = sqrt(maxX^2 + maxY^2) + 1

然后你的rotAngle和angle:

rotAngle <- 5
angle <- seq(0, 359, rotAngle)

将注意力集中在for循环上，第一个line调用具有变化和不变化的元素。让我们创建一个空列表来保存我们的结果，在for循环之外创建，它将保存2x2个矩阵:

line_lst <- list()
for (j in 1:length(angle)) {
line_lst[[j]] <- matrix(nrow = 2, ncol=2)
line_lst[[j]][1,1] <- meanX
line_lst[[j]][1,2] <- meanY
line_lst[[j]][2,1] <- meanX + line_length * cos((pi/180)*angle[j])
line_lst[[j]][2,2] <- meanY + line_length * sin((pi/180)*angle[j])
}
line_lst[[1]]
[,1]       [,2]
[1,] 1.225402e-06 0.09131118
[2,] 2.425684e+01 0.09131118
line_lst[[72]]
[,1]        [,2]
[1,] 1.225402e-06  0.09131118
[2,] 2.416454e+01 -2.02281169

这些似乎是合理的，这主要是我想要展示的，解释LHS[j] <- RHS[j]，我们在1中的迭代:长度(角度)。然后是线串，交点和点，同样创建一个空接收器，循环通过:

# here we have mismatch of establishing an `i` counter then
# counting `j`, which look close enough to tired eyes
# this will result in NULL(s)
linestring_lst <- list()
for (i in 1:length(line_lst)) { # this causes future error
linestring_lst[[j]] <- st_sf(st_sfc(st_linestring(line_lst[[j]]))) 
}
# simply keeping our accounting right, using all `i` or all `j`,
# or staying away from things that look alike and using `k` here
for (k in 1:length(line_lst)) { 
linestring_lst[[k]] <- st_sf(st_sfc(st_linestring(line_lst[[k]]))) 
}
intersection_lst <- list()
for (j in 1:length(linestring_lst)) {
intersection_lst[[j]] <- st_intersection(poly, linestring_lst[[j]])
}
intersect_points <- list()
for (j in 1:length(intersection_lst)) {
intersect_points[[j]] <- st_coordinates(intersection_lst[[j]])[2,c('X','Y')]
}

这里要记住的事情与for循环有关，在循环外创建接收器对象，索引LHS[j]和RHS[j] ([用于向量接收器，[[用于列表)。在独立完成了这些之后，你可以把它们都放到一个for循环中。

最后一步，将列表放到data.frame中，以便在ggplot中使用。

intersect_pts_df <- as.data.frame(do.call('rbind', intersect_points))
head(intersect_pts_df, n = 3)
X          Y
1 14.96993 0.09131118
2 15.56797 1.45333163
3 15.87039 2.88968964