使用 R 高效排序随机生成的横断面
Use R to Efficiently Order Randomly Generated Transects
问题
我正在寻找一种方法来有效地对随机选择的采样横断面进行排序，这些横断面发生在固定对象周围。这些横断面一旦生成，就需要以一种在空间上有意义的方式进行排序，以使行进的距离最小化。这将通过确保当前横断面的终点尽可能接近下一条横断面的起点来完成。此外，none 的横断面可以重复。
因为要订购数以千计的横断面，手动完成这是一项非常繁琐的任务，我正在尝试使用 R 来自动执行此过程。我已经生成了横断面，每个横断面都有一个起点和终点，其位置使用 360 度系统表示（例如，0 是北，90 是东，180 是南，270 是西）。我还生成了一些代码，似乎指示下一个样带的起点和 ID，但此代码存在一些问题：(1) 它可能会根据所考虑的起点和终点产生错误，(2 ) 它没有实现我最终需要它实现的目标，并且 (3) 代码本身似乎过于复杂，我不禁想知道是否有更直接的方法来做到这一点。
理想情况下，代码会导致横断面被重新排序，以便它们匹配它们应该飞行的顺序而不是它们最初输入的顺序。
数据
为简单起见，我们假设只有 10 个样带需要订购。
# Transect ID for the start point
StID <- c(seq(1, 10, 1))
# Location of transect start point, based on a 360-degree circle
StPt <- c(342.1, 189.3, 116.5, 67.9, 72, 208.4, 173.2, 97.8, 168.7, 138.2)
# Transect ID for the end point
EndID <- c(seq(1, 10, 1))
# Location of transect start point, based on a 360-degree circle
EndPt <- c(122.3, 313.9, 198.7, 160.4, 166, 26.7, 312.7, 273.7, 288.8, 287.5)
# Dataframe
df <- cbind.data.frame(StPt, StID, EndPt, EndID)
我试过的
请忽略这段代码，必须有更好的方法，但它并没有真正达到预期的结果。现在我正在使用一个嵌套的 for 循环，它很难直观地遵循但代表了我迄今为止的最佳尝试。
# Create two new columns that will be populated using a loop
df$StPt_Next <- NA
df$ID_Next <- NA
# Also create a list to be populated as end and start points are matched
used <- c(df$StPt[1]) #puts the start point of transect #1 into the used vector since we will start with 1 and do not want to have it used again
# Then, for every row in the dataframe...
for (i in seq(1,length(df$EndPt)-1, 1)){ # Selects all rows except the last one as the last transect should have no "next" transect
# generate some print statements to indicate that the script is indeed running while you wait....
print(paste("######## ENDPOINT", i, ":", df$EndPt[i], " ########"))
print(paste("searching for a start point that fits criteria to follow this endpoint",sep=""))
# sequentially select each end point
valueEndPt <- df[i,1]
# and order the index by taking the absolute difference of end and start points and, if this value is greater than 180, also subtract from 360 so all differences are less than 180, then order differences from smallest to largest
orderx <- order(ifelse(360-abs(df$StPt-valueEndPt) > 180,
abs(df$StPt-valueEndPt),
360-abs(df$StPt-valueEndPt)))
tmp <- as.data.frame(orderx)
# specify index value
index=1
# for as long as there is an "NA" present in the StPt_Next created before for loop...
while (is.na(df$StPt_Next[i])) {
#select the value of the ordered index in sequential order
j=orderx[index]
# if the start point associated with a given index is present in the list of used values...
if (df$StPt[j] %in% used){
# then have R print a statement indicate this is the case...
print(paste("passing ",df$StPt[j], " as it has already been used",sep=""))
# and move onto the next index
index=index+1
# break statement intended to skip the remainder of the code for values that have already been used
next
# if the start point associated with a given index is not present in the list of used values...
} else {
# then identify the start point value associated with that index ID...
valueStPt <- df$StPt[j]
# and have R print a statement indicating an attempt is being made to use the next value
print(paste("trying ",df$StPt[j],sep=""))
# if the end transect number is different from the start end transect number...
if (df$EndID[i] != df$StID[j]) {
# then put the start point in the new column...
df$StPt_Next[i] <- df$StPt[j]
# note which record this start point came from for ease of reference/troubleshooting...
df$ID_Next[i] <- j
# have R print a statement that indicates a value for the new column has beed selected...
print(paste("using ",df$StPt[j],sep=""))
# and add that start point to the list of used ones
used <- c(used,df$StPt[j])
# otherwise, if the end transect number matches the start end transect number...
} else {
# keep NA in this column and try again
df$StPt_Next[i] <- NA
# and indicate that this particular matched pair can not be used
print(paste("cant use ",valueStPt," as the column EndID (related to index in EndPt) and StID (related to index in StPt) values are matching",sep=""))
}# end if else statement to ensure that start and end points come from different transects
# and move onto the next index
index=index+1
}# end if else statement to determine if a given start point still needs to be used
}# end while loop to identify if there are still NA's in the new column
}# end for loop
输出
当代码没有产生显式错误时，例如对于提供的示例数据，输出如下：
StPt StID EndPt EndID StPt_Next ID_Next
1 342.1 1 122.3 1 67.9 4
2 189.3 2 313.9 2 173.2 7
3 116.5 3 198.7 3 97.8 8
4 67.9 4 160.4 4 72.0 5
5 72.0 5 166.0 5 116.5 3
6 208.4 6 26.7 6 189.3 2
7 173.2 7 312.7 7 168.7 9
8 97.8 8 273.7 8 138.2 10
9 168.7 9 288.8 9 208.4 6
10 138.2 10 287.5 10 NA NA
最后两列由代码生成并添加到原始数据框中。 StPt_Next 具有下一个最近起点的位置，ID_Next 指示与下一个起点位置关联的 transectID。 ID_Next 列表示横断面的飞行顺序如下 1,4,5,3,8,10,NA（又名结束），2,7,9,6 形成自己的循环可以追溯到 2。
有两个具体问题我解决不了：
(1) 存在形成一个连续序列链的问题。我认为这与 1 是起始横断面而 10 是最后横断面有关，无论如何，但不知道如何在代码中指示倒数第二个横断面必须与 10 匹配，以便序列包括所有 10 个横断面在代表最终终点的 "NA" 处终止之前。
(2) 为了真正使这个过程自动化，在修复了由于 "NA" 作为 ID_next 过早引入而导致的序列提前终止之后，将创建一个新列，它将允许根据最有效的进展而不是它们 EndID/StartID 的原始顺序对横断面进行重新排序。
预期结果
如果我们假装我们只有 6 个样带需要订购，而忽略由于过早引入 "NA" 而无法订购的 4 个样带，这将是预期的结果：
StPt StID EndPt EndID StPt_Next ID_Next TransNum
1 342.1 1 122.3 1 67.9 4 1
4 67.9 4 160.4 4 72.0 5 2
5 72.0 5 166.0 5 116.5 3 3
3 116.5 3 198.7 3 97.8 8 4
8 97.8 8 273.7 8 138.2 10 5
10 138.2 10 287.5 10 NA NA 6
编辑：关于代码明确生成的错误消息的说明
如前所述，代码有一些缺陷。另一个缺陷是，在尝试订购大量横断面时，它通常会产生错误。我不完全确定错误是在过程中的哪一步产生的，但我猜测这与无法匹配最后几个横断面有关，可能是由于不符合 "orderx" 规定的标准.打印语句说 "trying NA" 而不是数据库中的起点，这会导致此错误："Error in if (df$EndID[i] != df$StID[j]) { : missing value where TRUE/FALSE needed"。我猜还需要另一个 if-else 语句以某种方式指示 "if the remaining points do not meet the orderx criteria, then just force them to match up with whatever transect remains so that everything is assigned a StPt_Next and ID_Next".
这是一个会产生错误的更大的数据集：
EndPt <- c(158.7,245.1,187.1,298.2,346.8,317.2,74.5,274.2,153.4,246.7,193.6,302.3,6.8,359.1,235.4,134.5,111.2,240.5,359.2,121.3,224.5,212.6,155.1,353.1,181.7,334,249.3,43.9,38.5,75.7,344.3,45.1,285.7,155.5,183.8,60.6,301,132.1,75.9,112,342.1,302.1,288.1,47.4,331.3,3.4,185.3,62,323.7,188,313.1,171.6,187.6,291.4,19.2,210.3,93.3,24.8,83.1,193.8,112.7,204.3,223.3,210.7,201.2,41.3,79.7,175.4,260.7,279.5,82.4,200.2,254.2,228.9,1.4,299.9,102.7,123.7,172.9,23.2,207.3,320.1,344.6,39.9,223.8,106.6,156.6,45.7,236.3,98.1,337.2,296.1,194,307.1,86.6,65.5,86.6,296.4,94.7,279.9)
StPt <- c(56.3,158.1,82.4,185.5,243.9,195.6,335,167,39.4,151.7,99.8,177.2,246.8,266.1,118.2,358.6,357.9,99.6,209.9,342.8,106.5,86.4,35.7,200.6,65.6,212.5,159.1,297,285.9,300.9,177,245.2,153.1,8.1,76.5,322.4,190.8,35.2,342.6,8.8,244.6,202,176.2,308.3,184.2,267.2,26.6,293.8,167.3,30.5,176,74.3,96.9,186.7,288.2,62.6,331.4,254.7,324.1,73.4,16.4,64,110.9,74.4,69.8,298.8,336.6,58.8,170.1,173.2,330.8,92.6,129.2,124.7,262.3,140.4,321.2,34,79.5,263,66.4,172.8,205.5,288,98.5,335.2,38.7,289.7,112.7,350.7,243.2,185.4,63.9,170.3,326.3,322.9,320.6,199.2,287.1,158.1)
EndID <- c(seq(1, 100, 1))
StID <- c(seq(1, 100, 1))
df <- cbind.data.frame(StPt, StID, EndPt, EndID)
如有任何建议，我们将不胜感激！
正如@chinsoon12 指出隐藏在您的问题中的，您有一个（不对称的）旅行商问题。出现不对称是因为你的transecs的起点和终点不同。
ATSP 是著名的 NP 完全问题。因此，即使对于中等规模的问题，精确的解决方案也非常困难（有关更多信息，请参阅 wikipedia）。因此，在大多数情况下，我们能做的最好的就是近似或启发式。正如您提到的，有数千条横断面，这至少是一个中等规模的问题。
不是从一开始就编写一个 ATSP 近似算法，而是有一个现有的 R 的 TSP 库。这包括几个近似算法。参考文档是 here.
以下是我使用的TSP包应用于你的问题。从设置开始（假设我有 运行 StPt、StID、EndPt 和 EndID，如您的问题。
install.packages("TSP")
library(TSP)
library(dplyr)
# Dataframe
df <- cbind.data.frame(StPt, StID, EndPt, EndID)
# filter to 6 example nodes for requested comparison
df = df %>% filter(StID %in% c(1,3,4,5,8,10))
我们将从距离矩阵使用 ATSP。矩阵中的位置 [row,col] 是从样带 row（结束）到样带 col（开始）的 cost/distance。此代码创建整个距离矩阵。
# distance calculation
transec_distance = function(end,start){
abs_dist = abs(start-end)
ifelse(360-abs_dist > 180, abs_dist, 360-abs_dist)
}
# distance matrix
matrix_distance = matrix(data = NA, nrow = nrow(df), ncol = nrow(df))
for(start_id in 1:nrow(df)){
start_point = df[start_id,'StPt']
for(end_id in 1:nrow(df)){
end_point = df[end_id,'EndPt']
matrix_distance[end_id,start_id] = transec_distance(end_point, start_point)
}
}
请注意，有更有效的方法来构建距离矩阵。但是，出于清晰度的考虑，我选择了这种方法。根据您的计算机和横断面的确切数量，此代码可能 运行 非常慢。
另外，请注意该矩阵的大小是横断面数量的二次方。所以对于大量的样带，你会发现内存不够。
求解很平淡。距离矩阵变成一个 ATSP 对象，ATSP 对象被传递给求解器。然后我们继续将 ordering/traveling 信息添加到原始 df.
answer = solve_TSP(as.ATSP(matrix_distance))
# get length of cycle
print(answer)
# sort df to same order as solution
df_w_answer = df[as.numeric(answer),]
# add info about next transect to each transect
df_w_answer = df_w_answer %>%
mutate(visit_order = 1:nrow(df_w_answer)) %>%
mutate(next_StID = lead(StID, order_by = visit_order),
next_StPt = lead(StPt, order_by = visit_order))
# add info about next transect to each transect (for final transect)
df_w_answer[df_w_answer$visit_order == nrow(df_w_answer),'next_StID'] =
df_w_answer[df_w_answer$visit_order == 1,'StID']
df_w_answer[df_w_answer$visit_order == nrow(df_w_answer),'next_StPt'] =
df_w_answer[df_w_answer$visit_order == 1,'StPt']
# compute distance between end of each transect and start of next
df_w_answer = df_w_answer %>% mutate(dist_between = transec_distance(EndPt, next_StPt))
此时我们有一个循环。您可以选择任何节点作为起点，按照 df 中给出的顺序：从 EndID 到 next_StID，您将覆盖最小距离（一个很好的近似值）中的每个横断面。
然而，在您的 'intended outcome' 中，您有一个路径解决方案（例如，从横断面 1 开始到横断面 10 结束）。我们可以通过排除单个最昂贵的转换将循环变成路径：
# as path (without returning to start)
min_distance = sum(df_w_answer$dist_between) - max(df_w_answer$dist_between)
path_start = df_w_answer[df_w_answer$dist_between == max(df_w_answer$dist_between), 'next_StID']
path_end = df_w_answer[df_w_answer$dist_between == max(df_w_answer$dist_between), 'EndID']
print(sprintf("minimum cost path = %.2f, starting at node %d, ending at node %d",
min_distance, path_start, path_end))
运行 以上所有内容为我提供了一个不同但更好的答案来满足您的预期结果。我得到以下顺序：1 --> 5 --> 8 --> 4 --> 3 --> 10 --> 1.
- 你从横断面 1 到横断面 10 的总距离为 428，如果我们也从横断面 10 返回到横断面 1，使这成为一个循环，则总距离为 483。
- 使用 R 中的 TSP 包，我们得到一条从 1 到 10 的路径，总距离为 377，循环为 431。
- 如果我们改为从节点 4 开始并在节点 8 结束，我们得到的总距离为 277。
一些额外的节点：
- 并非所有的 TSP 求解器都是确定性的，因此如果您再次 运行 或 运行 以不同的顺序输入行，您的答案可能会有所不同。
- TSP 是一个比您描述的横断面问题更普遍的问题。有可能你的问题有足够的 additional/special 特征，这意味着它可以在合理的时间内完美解决。但这会将您的问题转移到数学领域。
- 如果您运行创建距离矩阵时内存不足，请查看 TSP 包的文档。它包含几个使用地理坐标而不是距离矩阵作为输入的示例。这是一个小得多的输入大小（大概是包计算动态距离）所以如果你将起点和终点转换为坐标并指定欧几里德（或其他一些常见的距离函数）你可以绕过（一些）计算机内存限制.
另一个使用TSP包的版本...
这是设置。
library(TSP)
planeDim = 15
nTransects = 26
# generate some random transect beginning points in a plane, the size of which
# is defined by planeDim
b = cbind(runif(nTransects)*planeDim, runif(nTransects)*planeDim)
# generate some random transect ending points that are a distance of 1 from each
# beginning point
e = t(
apply(
b,
1,
function(x) {
bearing = runif(1)*2*pi
x + c(cos(bearing), sin(bearing))
}
)
)
为了好玩，我们可以可视化横断面：
# make an empty graph space
plot(1,1, xlim = c(-1, planeDim + 1), ylim = c(-1, planeDim + 1), ty = "n")
# plot the beginning of each transect as a green point, the end as a red point,
# with a thick grey line representing the transect
for(i in 1:nrow(e)) {
xs = c(b[i,1], e[i,1])
ys = c(b[i,2], e[i,2])
lines(xs, ys, col = "light grey", lwd = 4)
points(xs, ys, col = c("green", "red"), pch = 20, cex = 1.5)
text(mean(xs), mean(ys), letters[i])
}
所以给定一个 x,y 对矩阵 ("b") 作为起始点和一个 x,y 矩阵
对 ("e") 每个横断面的端点，解决方案是...
# a function that calculates the distance from all endpoints in the ePts matrix
# to the single beginning point in bPt
dist = function(ePts, bPt) {
# apply pythagorean theorem to calculate the distance between every end point
# in the matrix ePts to the point bPt
apply(ePts, 1, function(p) sum((p - bPt)^2)^0.5)
}
# apply the "dist" function to all begining points to create the distance
# matrix. since the distance from the end of transect "foo" to the beginning of
# "bar" is not the same as from the end of "bar" to the beginning of "foo," we
# have an asymmetric travelling sales person problem. Therefore, distance
# matrix is directional. The distances at any position in the matrix must be
# the distance from the transect shown in the row label and to the transect
# shown in the column label.
distMatrix = apply(b, 1, FUN = dist, ePts = e)
# for covenience, we'll labels the trasects a to z
dimnames(distMatrix) = list(letters, letters)
# set the distance between the beginning and end of each transect to zero so
# that there is no "cost" to walking the transect
diag(distMatrix) = 0
这里是距离矩阵的左上角：
> distMatrix[1:6, 1:6]
a b c d e f
a 0.00000 15.4287270 12.637979 12.269356 15.666710 12.3919715
b 13.58821 0.0000000 5.356411 13.840444 1.238677 12.6512352
c 12.48161 6.3086852 0.000000 8.427033 6.382304 7.1387840
d 10.69748 13.5936114 7.708183 0.000000 13.718517 0.9836146
e 14.00920 0.7736654 5.980220 14.470826 0.000000 13.2809601
f 12.24503 12.8987043 6.984763 2.182829 12.993283 0.0000000
现在 TSP 包中的三行代码解决了这个问题。
atsp = as.ATSP(distMatrix)
tour = solve_TSP(atsp)
# assume we want to start our circuit at transect "a".
path = cut_tour(tour, "a", exclude_cut = F)
变量 path 显示了您访问横断面的顺序：
> path
a w x q i o l d f s h y g v t k c m e b p u z j r n
1 23 24 17 9 15 12 4 6 19 8 25 7 22 20 11 3 13 5 2 16 21 26 10 18 14
我们可以将路径添加到可视化：
for(i in 1:(length(path)-1)) {
lines(c(e[path[i],1], b[path[i+1],1]), c(e[path[i],2], b[path[i+1], 2]), lty = "dotted")
}
感谢大家的建议，@Simon 的解决方案最适合我的 OP。 @Geoffrey 使用 x,y 坐标的实际方法很棒，因为它允许绘制横断面和旅行顺序。因此，我发布了一个混合解决方案，该解决方案是使用他们两人的代码生成的，以及其他评论和代码来详细说明该过程并获得我想要的实际最终结果。我不确定这是否会在未来帮助任何人，但是，由于没有答案提供 100% 解决我的问题的解决方案，我想我会分享我的想法。
正如其他人所指出的，这是一种旅行推销员问题。它是不对称的，因为从样线终点 "t" 到样线起点 "t+1" 的距离与从终点样线 "t+1" 到样线起点 [=22] 的距离不同=].此外，它是 "path" 解决方案而不是 "cycle" 解决方案。
#=========================================
# Packages
#=========================================
library(TSP)
library(useful)
library(dplyr)
#=========================================
# Full dataset for testing
#=========================================
EndPt <- c(158.7,245.1,187.1,298.2,346.8,317.2,74.5,274.2,153.4,246.7,193.6,302.3,6.8,359.1,235.4,134.5,111.2,240.5,359.2,121.3,224.5,212.6,155.1,353.1,181.7,334,249.3,43.9,38.5,75.7,344.3,45.1,285.7,155.5,183.8,60.6,301,132.1,75.9,112,342.1,302.1,288.1,47.4,331.3,3.4,185.3,62,323.7,188,313.1,171.6,187.6,291.4,19.2,210.3,93.3,24.8,83.1,193.8,112.7,204.3,223.3,210.7,201.2,41.3,79.7,175.4,260.7,279.5,82.4,200.2,254.2,228.9,1.4,299.9,102.7,123.7,172.9,23.2,207.3,320.1,344.6,39.9,223.8,106.6,156.6,45.7,236.3,98.1,337.2,296.1,194,307.1,86.6,65.5,86.6,296.4,94.7,279.9)
StPt <- c(56.3,158.1,82.4,185.5,243.9,195.6,335,167,39.4,151.7,99.8,177.2,246.8,266.1,118.2,358.6,357.9,99.6,209.9,342.8,106.5,86.4,35.7,200.6,65.6,212.5,159.1,297,285.9,300.9,177,245.2,153.1,8.1,76.5,322.4,190.8,35.2,342.6,8.8,244.6,202,176.2,308.3,184.2,267.2,26.6,293.8,167.3,30.5,176,74.3,96.9,186.7,288.2,62.6,331.4,254.7,324.1,73.4,16.4,64,110.9,74.4,69.8,298.8,336.6,58.8,170.1,173.2,330.8,92.6,129.2,124.7,262.3,140.4,321.2,34,79.5,263,66.4,172.8,205.5,288,98.5,335.2,38.7,289.7,112.7,350.7,243.2,185.4,63.9,170.3,326.3,322.9,320.6,199.2,287.1,158.1)
EndID <- c(seq(1, 100, 1))
StID <- c(seq(1, 100, 1))
df <- cbind.data.frame(StPt, StID, EndPt, EndID)
#=========================================
# Convert polar coordinates to cartesian x,y data
#=========================================
# Area that the transect occupy in space only used for graphing
planeDim <- 1
# Number of transects
nTransects <- 100
# Convert 360-degree polar coordinates to x,y cartesian coordinates to facilitate calculating a distance matrix based on the Pythagorean theorem
EndX <- as.matrix(pol2cart(planeDim, EndPt, degrees = TRUE)["x"])
EndY <- as.matrix(pol2cart(planeDim, EndPt, degrees = TRUE)["y"])
StX <- as.matrix(pol2cart(planeDim, StPt, degrees = TRUE)["x"])
StY <- as.matrix(pol2cart(planeDim, StPt, degrees = TRUE)["y"])
# Matrix of x,y pairs for the beginning ("b") and end ("e") points of each transect
b <- cbind(c(StX), c(StY))
e <- cbind(c(EndX), c(EndY))
#=========================================
# Function to calculate the distance from all endpoints in the ePts matrix to a single beginning point in bPt
#=========================================
dist <- function(ePts, bPt) {
# Use the Pythagorean theorem to calculate the hypotenuse (i.e., distance) between every end point in the matrix ePts to the point bPt
apply(ePts, 1, function(p) sum((p - bPt)^2)^0.5)
}
#=========================================
# Distance matrix
#=========================================
# Apply the "dist" function to all beginning points to create a matrix that has the distance between every start and endpoint
## Note: because this is an asymmetric traveling salesperson problem, the distance matrix is directional, thus, the distances at any position in the matrix must be the distance from the transect shown in the row label and to the transect shown in the column label
distMatrix <- apply(b, 1, FUN = dist, ePts = e)
## Set the distance between the beginning and end of each transect to zero so that there is no "cost" to walking the transect
diag(distMatrix) <- 0
#=========================================
# Solve asymmetric TSP
#=========================================
# This creates an instance of the asymmetric traveling salesperson (ASTP)
atsp <- as.ATSP(distMatrix)
# This creates an object of Class Tour that travels to all of the points
## In this case, the repetitive_nn produces the smallest overall and transect-to-transect
tour <- solve_TSP(atsp, method = "repetitive_nn")
#=========================================
# Create a path by cutting the tour at the most "expensive" transition
#=========================================
# Sort the original data frame to match the order of the solution
dfTour = df[as.numeric(tour),]
# Add the following columns to the original dataframe:
dfTour = dfTour %>%
# Assign visit order (1 to 100, ascending)
mutate(visit_order = 1:nrow(dfTour)) %>%
# The ID of the next transect to move to
mutate(next_StID = lead(StID, order_by = visit_order),
# The angle of the start point for the next transect
next_StPt = lead(StPt, order_by = visit_order))
# lead() generates the NA's in the last record for next_StID, next_StPt, replace these by adding that information
dfTour[dfTour$visit_order == nrow(dfTour),'next_StID'] <-
dfTour[dfTour$visit_order == 1,'StID']
dfTour[dfTour$visit_order == nrow(dfTour),'next_StPt'] <-
dfTour[dfTour$visit_order == 1,'StPt']
# Function to calculate distance for 360 degrees rather than x,y coordinates
transect_distance <- function(end,start){
abs_dist = abs(start-end)
ifelse(360-abs_dist > 180, abs_dist, 360-abs_dist)
}
# Compute distance between end of each transect and start of next using polar coordinates
dfTour = dfTour %>% mutate(dist_between = transect_distance(EndPt, next_StPt))
# Identify the longest transition point for breaking the cycle
min_distance <- sum(dfTour$dist_between) - max(dfTour$dist_between)
path_start <- dfTour[dfTour$dist_between == max(dfTour$dist_between), 'next_StID']
path_end <- dfTour[dfTour$dist_between == max(dfTour$dist_between), 'EndID']
# Make a statement about the least cost path
print(sprintf("minimum cost path = %.2f, starting at node %d, ending at node %d",
min_distance, path_start, path_end))
# The variable path shows the order in which you should visit the transects
path <- cut_tour(tour, path_start, exclude_cut = F)
# Arrange df from smallest to largest travel distance
tmp1 <- dfTour %>%
arrange(dist_between)
# Change dist_between and visit_order to NA for transect with the largest distance to break cycle
# (e.g., we will not travel this distance, this represents the path endpoint)
tmp1[length(dfTour$dist_between):length(dfTour$dist_between),8] <- NA
tmp1[length(dfTour$dist_between):length(dfTour$dist_between),5] <- NA
# Set df order back to ascending by visit order
tmp2 <- tmp1 %>%
arrange(visit_order)
# Detect the break in a sequence of visit_order introduced by the NA (e.g., 1,2,3....5,6) and mark groups before the break with 0 and after the break with 1 in the "cont_per" column
tmp2$cont_per <- cumsum(!c(TRUE, diff(tmp2$visit_order)==1))
# Sort "cont_per" such that the records following the break become the beginning of the path and the ones following the break represent the middle orders and the point with the NA being assigned the last visit order, and assign a new visit order
tmp3 <- tmp2%>%
arrange(desc(cont_per))%>%
mutate(visit_order_FINAL=seq(1, length(tmp2$visit_order), 1))
# Datframe ordered by progression of transects
trans_order <- cbind.data.frame(tmp3[2], tmp3[1], tmp3[4], tmp3[3], tmp3[6], tmp3[7], tmp3[8], tmp3[10])
# Insert NAs for "next" info for final transect
trans_order[nrow(trans_order),'next_StPt'] <- NA
trans_order[nrow(trans_order), 'next_StID'] <- NA
#=========================================
# View data
#=========================================
head(trans_order)
#=========================================
# Plot
#=========================================
#For fun, we can visualize the transects:
# make an empty graph space
plot(1,1, xlim = c(-planeDim-0.1, planeDim+0.1), ylim = c(-planeDim-0.1, planeDim+0.1), ty = "n")
# plot the beginning of each transect as a green point, the end as a red point,
and a grey line representing the transect
for(i in 1:nrow(e)) {
xs = c(b[i,1], e[i,1])
ys = c(b[i,2], e[i,2])
lines(xs, ys, col = "light grey", lwd = 1, lty = 1)
points(xs, ys, col = c("green", "red"), pch = 1, cex = 1)
#text((xs), (ys), i)
}
# Add the path to the visualization
for(i in 1:(length(path)-1)) {
# This makes a line between the x coordinates for the end point of path i and beginning point of path i+1
lines(c(e[path[i],1], b[path[i+1],1]), c(e[path[i],2], b[path[i+1], 2]), lty = 1, lwd=1)
}
这就是最终结果的样子
问题
我正在寻找一种方法来有效地对随机选择的采样横断面进行排序，这些横断面发生在固定对象周围。这些横断面一旦生成，就需要以一种在空间上有意义的方式进行排序，以使行进的距离最小化。这将通过确保当前横断面的终点尽可能接近下一条横断面的起点来完成。此外，none 的横断面可以重复。
因为要订购数以千计的横断面，手动完成这是一项非常繁琐的任务，我正在尝试使用 R 来自动执行此过程。我已经生成了横断面，每个横断面都有一个起点和终点，其位置使用 360 度系统表示（例如，0 是北，90 是东，180 是南，270 是西）。我还生成了一些代码，似乎指示下一个样带的起点和 ID，但此代码存在一些问题：(1) 它可能会根据所考虑的起点和终点产生错误，(2 ) 它没有实现我最终需要它实现的目标，并且 (3) 代码本身似乎过于复杂，我不禁想知道是否有更直接的方法来做到这一点。
理想情况下，代码会导致横断面被重新排序，以便它们匹配它们应该飞行的顺序而不是它们最初输入的顺序。
数据
为简单起见，我们假设只有 10 个样带需要订购。
# Transect ID for the start point
StID <- c(seq(1, 10, 1))
# Location of transect start point, based on a 360-degree circle
StPt <- c(342.1, 189.3, 116.5, 67.9, 72, 208.4, 173.2, 97.8, 168.7, 138.2)
# Transect ID for the end point
EndID <- c(seq(1, 10, 1))
# Location of transect start point, based on a 360-degree circle
EndPt <- c(122.3, 313.9, 198.7, 160.4, 166, 26.7, 312.7, 273.7, 288.8, 287.5)
# Dataframe
df <- cbind.data.frame(StPt, StID, EndPt, EndID)
我试过的
请忽略这段代码，必须有更好的方法，但它并没有真正达到预期的结果。现在我正在使用一个嵌套的 for 循环，它很难直观地遵循但代表了我迄今为止的最佳尝试。
# Create two new columns that will be populated using a loop
df$StPt_Next <- NA
df$ID_Next <- NA
# Also create a list to be populated as end and start points are matched
used <- c(df$StPt[1]) #puts the start point of transect #1 into the used vector since we will start with 1 and do not want to have it used again
# Then, for every row in the dataframe...
for (i in seq(1,length(df$EndPt)-1, 1)){ # Selects all rows except the last one as the last transect should have no "next" transect
# generate some print statements to indicate that the script is indeed running while you wait....
print(paste("######## ENDPOINT", i, ":", df$EndPt[i], " ########"))
print(paste("searching for a start point that fits criteria to follow this endpoint",sep=""))
# sequentially select each end point
valueEndPt <- df[i,1]
# and order the index by taking the absolute difference of end and start points and, if this value is greater than 180, also subtract from 360 so all differences are less than 180, then order differences from smallest to largest
orderx <- order(ifelse(360-abs(df$StPt-valueEndPt) > 180,
abs(df$StPt-valueEndPt),
360-abs(df$StPt-valueEndPt)))
tmp <- as.data.frame(orderx)
# specify index value
index=1
# for as long as there is an "NA" present in the StPt_Next created before for loop...
while (is.na(df$StPt_Next[i])) {
#select the value of the ordered index in sequential order
j=orderx[index]
# if the start point associated with a given index is present in the list of used values...
if (df$StPt[j] %in% used){
# then have R print a statement indicate this is the case...
print(paste("passing ",df$StPt[j], " as it has already been used",sep=""))
# and move onto the next index
index=index+1
# break statement intended to skip the remainder of the code for values that have already been used
next
# if the start point associated with a given index is not present in the list of used values...
} else {
# then identify the start point value associated with that index ID...
valueStPt <- df$StPt[j]
# and have R print a statement indicating an attempt is being made to use the next value
print(paste("trying ",df$StPt[j],sep=""))
# if the end transect number is different from the start end transect number...
if (df$EndID[i] != df$StID[j]) {
# then put the start point in the new column...
df$StPt_Next[i] <- df$StPt[j]
# note which record this start point came from for ease of reference/troubleshooting...
df$ID_Next[i] <- j
# have R print a statement that indicates a value for the new column has beed selected...
print(paste("using ",df$StPt[j],sep=""))
# and add that start point to the list of used ones
used <- c(used,df$StPt[j])
# otherwise, if the end transect number matches the start end transect number...
} else {
# keep NA in this column and try again
df$StPt_Next[i] <- NA
# and indicate that this particular matched pair can not be used
print(paste("cant use ",valueStPt," as the column EndID (related to index in EndPt) and StID (related to index in StPt) values are matching",sep=""))
}# end if else statement to ensure that start and end points come from different transects
# and move onto the next index
index=index+1
}# end if else statement to determine if a given start point still needs to be used
}# end while loop to identify if there are still NA's in the new column
}# end for loop
输出
当代码没有产生显式错误时，例如对于提供的示例数据，输出如下：
StPt StID EndPt EndID StPt_Next ID_Next
1 342.1 1 122.3 1 67.9 4
2 189.3 2 313.9 2 173.2 7
3 116.5 3 198.7 3 97.8 8
4 67.9 4 160.4 4 72.0 5
5 72.0 5 166.0 5 116.5 3
6 208.4 6 26.7 6 189.3 2
7 173.2 7 312.7 7 168.7 9
8 97.8 8 273.7 8 138.2 10
9 168.7 9 288.8 9 208.4 6
10 138.2 10 287.5 10 NA NA
最后两列由代码生成并添加到原始数据框中。 StPt_Next 具有下一个最近起点的位置，ID_Next 指示与下一个起点位置关联的 transectID。 ID_Next 列表示横断面的飞行顺序如下 1,4,5,3,8,10,NA（又名结束），2,7,9,6 形成自己的循环可以追溯到 2。
有两个具体问题我解决不了：
(1) 存在形成一个连续序列链的问题。我认为这与 1 是起始横断面而 10 是最后横断面有关，无论如何，但不知道如何在代码中指示倒数第二个横断面必须与 10 匹配，以便序列包括所有 10 个横断面在代表最终终点的 "NA" 处终止之前。
(2) 为了真正使这个过程自动化，在修复了由于 "NA" 作为 ID_next 过早引入而导致的序列提前终止之后，将创建一个新列，它将允许根据最有效的进展而不是它们 EndID/StartID 的原始顺序对横断面进行重新排序。
预期结果
如果我们假装我们只有 6 个样带需要订购，而忽略由于过早引入 "NA" 而无法订购的 4 个样带，这将是预期的结果：
StPt StID EndPt EndID StPt_Next ID_Next TransNum
1 342.1 1 122.3 1 67.9 4 1
4 67.9 4 160.4 4 72.0 5 2
5 72.0 5 166.0 5 116.5 3 3
3 116.5 3 198.7 3 97.8 8 4
8 97.8 8 273.7 8 138.2 10 5
10 138.2 10 287.5 10 NA NA 6
编辑：关于代码明确生成的错误消息的说明
如前所述，代码有一些缺陷。另一个缺陷是，在尝试订购大量横断面时，它通常会产生错误。我不完全确定错误是在过程中的哪一步产生的，但我猜测这与无法匹配最后几个横断面有关，可能是由于不符合 "orderx" 规定的标准.打印语句说 "trying NA" 而不是数据库中的起点，这会导致此错误："Error in if (df$EndID[i] != df$StID[j]) { : missing value where TRUE/FALSE needed"。我猜还需要另一个 if-else 语句以某种方式指示 "if the remaining points do not meet the orderx criteria, then just force them to match up with whatever transect remains so that everything is assigned a StPt_Next and ID_Next".
这是一个会产生错误的更大的数据集：
EndPt <- c(158.7,245.1,187.1,298.2,346.8,317.2,74.5,274.2,153.4,246.7,193.6,302.3,6.8,359.1,235.4,134.5,111.2,240.5,359.2,121.3,224.5,212.6,155.1,353.1,181.7,334,249.3,43.9,38.5,75.7,344.3,45.1,285.7,155.5,183.8,60.6,301,132.1,75.9,112,342.1,302.1,288.1,47.4,331.3,3.4,185.3,62,323.7,188,313.1,171.6,187.6,291.4,19.2,210.3,93.3,24.8,83.1,193.8,112.7,204.3,223.3,210.7,201.2,41.3,79.7,175.4,260.7,279.5,82.4,200.2,254.2,228.9,1.4,299.9,102.7,123.7,172.9,23.2,207.3,320.1,344.6,39.9,223.8,106.6,156.6,45.7,236.3,98.1,337.2,296.1,194,307.1,86.6,65.5,86.6,296.4,94.7,279.9)
StPt <- c(56.3,158.1,82.4,185.5,243.9,195.6,335,167,39.4,151.7,99.8,177.2,246.8,266.1,118.2,358.6,357.9,99.6,209.9,342.8,106.5,86.4,35.7,200.6,65.6,212.5,159.1,297,285.9,300.9,177,245.2,153.1,8.1,76.5,322.4,190.8,35.2,342.6,8.8,244.6,202,176.2,308.3,184.2,267.2,26.6,293.8,167.3,30.5,176,74.3,96.9,186.7,288.2,62.6,331.4,254.7,324.1,73.4,16.4,64,110.9,74.4,69.8,298.8,336.6,58.8,170.1,173.2,330.8,92.6,129.2,124.7,262.3,140.4,321.2,34,79.5,263,66.4,172.8,205.5,288,98.5,335.2,38.7,289.7,112.7,350.7,243.2,185.4,63.9,170.3,326.3,322.9,320.6,199.2,287.1,158.1)
EndID <- c(seq(1, 100, 1))
StID <- c(seq(1, 100, 1))
df <- cbind.data.frame(StPt, StID, EndPt, EndID)
如有任何建议，我们将不胜感激！
正如@chinsoon12 指出隐藏在您的问题中的，您有一个（不对称的）旅行商问题。出现不对称是因为你的transecs的起点和终点不同。
ATSP 是著名的 NP 完全问题。因此，即使对于中等规模的问题，精确的解决方案也非常困难（有关更多信息，请参阅 wikipedia）。因此，在大多数情况下，我们能做的最好的就是近似或启发式。正如您提到的，有数千条横断面，这至少是一个中等规模的问题。
不是从一开始就编写一个 ATSP 近似算法，而是有一个现有的 R 的 TSP 库。这包括几个近似算法。参考文档是 here.
以下是我使用的TSP包应用于你的问题。从设置开始（假设我有 运行 StPt、StID、EndPt 和 EndID，如您的问题。
install.packages("TSP")
library(TSP)
library(dplyr)
# Dataframe
df <- cbind.data.frame(StPt, StID, EndPt, EndID)
# filter to 6 example nodes for requested comparison
df = df %>% filter(StID %in% c(1,3,4,5,8,10))
我们将从距离矩阵使用 ATSP。矩阵中的位置 [row,col] 是从样带 row（结束）到样带 col（开始）的 cost/distance。此代码创建整个距离矩阵。
# distance calculation
transec_distance = function(end,start){
abs_dist = abs(start-end)
ifelse(360-abs_dist > 180, abs_dist, 360-abs_dist)
}
# distance matrix
matrix_distance = matrix(data = NA, nrow = nrow(df), ncol = nrow(df))
for(start_id in 1:nrow(df)){
start_point = df[start_id,'StPt']
for(end_id in 1:nrow(df)){
end_point = df[end_id,'EndPt']
matrix_distance[end_id,start_id] = transec_distance(end_point, start_point)
}
}
请注意，有更有效的方法来构建距离矩阵。但是，出于清晰度的考虑，我选择了这种方法。根据您的计算机和横断面的确切数量，此代码可能 运行 非常慢。
另外，请注意该矩阵的大小是横断面数量的二次方。所以对于大量的样带，你会发现内存不够。
求解很平淡。距离矩阵变成一个 ATSP 对象，ATSP 对象被传递给求解器。然后我们继续将 ordering/traveling 信息添加到原始 df.
answer = solve_TSP(as.ATSP(matrix_distance))
# get length of cycle
print(answer)
# sort df to same order as solution
df_w_answer = df[as.numeric(answer),]
# add info about next transect to each transect
df_w_answer = df_w_answer %>%
mutate(visit_order = 1:nrow(df_w_answer)) %>%
mutate(next_StID = lead(StID, order_by = visit_order),
next_StPt = lead(StPt, order_by = visit_order))
# add info about next transect to each transect (for final transect)
df_w_answer[df_w_answer$visit_order == nrow(df_w_answer),'next_StID'] =
df_w_answer[df_w_answer$visit_order == 1,'StID']
df_w_answer[df_w_answer$visit_order == nrow(df_w_answer),'next_StPt'] =
df_w_answer[df_w_answer$visit_order == 1,'StPt']
# compute distance between end of each transect and start of next
df_w_answer = df_w_answer %>% mutate(dist_between = transec_distance(EndPt, next_StPt))
此时我们有一个循环。您可以选择任何节点作为起点，按照 df 中给出的顺序：从 EndID 到 next_StID，您将覆盖最小距离（一个很好的近似值）中的每个横断面。
然而，在您的 'intended outcome' 中，您有一个路径解决方案（例如，从横断面 1 开始到横断面 10 结束）。我们可以通过排除单个最昂贵的转换将循环变成路径：
# as path (without returning to start)
min_distance = sum(df_w_answer$dist_between) - max(df_w_answer$dist_between)
path_start = df_w_answer[df_w_answer$dist_between == max(df_w_answer$dist_between), 'next_StID']
path_end = df_w_answer[df_w_answer$dist_between == max(df_w_answer$dist_between), 'EndID']
print(sprintf("minimum cost path = %.2f, starting at node %d, ending at node %d",
min_distance, path_start, path_end))
运行 以上所有内容为我提供了一个不同但更好的答案来满足您的预期结果。我得到以下顺序：1 --> 5 --> 8 --> 4 --> 3 --> 10 --> 1.
- 你从横断面 1 到横断面 10 的总距离为 428，如果我们也从横断面 10 返回到横断面 1，使这成为一个循环，则总距离为 483。
- 使用 R 中的 TSP 包，我们得到一条从 1 到 10 的路径，总距离为 377，循环为 431。
- 如果我们改为从节点 4 开始并在节点 8 结束，我们得到的总距离为 277。
一些额外的节点：
- 并非所有的 TSP 求解器都是确定性的，因此如果您再次 运行 或 运行 以不同的顺序输入行，您的答案可能会有所不同。
- TSP 是一个比您描述的横断面问题更普遍的问题。有可能你的问题有足够的 additional/special 特征，这意味着它可以在合理的时间内完美解决。但这会将您的问题转移到数学领域。
- 如果您运行创建距离矩阵时内存不足，请查看 TSP 包的文档。它包含几个使用地理坐标而不是距离矩阵作为输入的示例。这是一个小得多的输入大小（大概是包计算动态距离）所以如果你将起点和终点转换为坐标并指定欧几里德（或其他一些常见的距离函数）你可以绕过（一些）计算机内存限制.
另一个使用TSP包的版本...
这是设置。
library(TSP)
planeDim = 15
nTransects = 26
# generate some random transect beginning points in a plane, the size of which
# is defined by planeDim
b = cbind(runif(nTransects)*planeDim, runif(nTransects)*planeDim)
# generate some random transect ending points that are a distance of 1 from each
# beginning point
e = t(
apply(
b,
1,
function(x) {
bearing = runif(1)*2*pi
x + c(cos(bearing), sin(bearing))
}
)
)
为了好玩，我们可以可视化横断面：
# make an empty graph space
plot(1,1, xlim = c(-1, planeDim + 1), ylim = c(-1, planeDim + 1), ty = "n")
# plot the beginning of each transect as a green point, the end as a red point,
# with a thick grey line representing the transect
for(i in 1:nrow(e)) {
xs = c(b[i,1], e[i,1])
ys = c(b[i,2], e[i,2])
lines(xs, ys, col = "light grey", lwd = 4)
points(xs, ys, col = c("green", "red"), pch = 20, cex = 1.5)
text(mean(xs), mean(ys), letters[i])
}
所以给定一个 x,y 对矩阵 ("b") 作为起始点和一个 x,y 矩阵 对 ("e") 每个横断面的端点，解决方案是...
# a function that calculates the distance from all endpoints in the ePts matrix
# to the single beginning point in bPt
dist = function(ePts, bPt) {
# apply pythagorean theorem to calculate the distance between every end point
# in the matrix ePts to the point bPt
apply(ePts, 1, function(p) sum((p - bPt)^2)^0.5)
}
# apply the "dist" function to all begining points to create the distance
# matrix. since the distance from the end of transect "foo" to the beginning of
# "bar" is not the same as from the end of "bar" to the beginning of "foo," we
# have an asymmetric travelling sales person problem. Therefore, distance
# matrix is directional. The distances at any position in the matrix must be
# the distance from the transect shown in the row label and to the transect
# shown in the column label.
distMatrix = apply(b, 1, FUN = dist, ePts = e)
# for covenience, we'll labels the trasects a to z
dimnames(distMatrix) = list(letters, letters)
# set the distance between the beginning and end of each transect to zero so
# that there is no "cost" to walking the transect
diag(distMatrix) = 0
这里是距离矩阵的左上角：
> distMatrix[1:6, 1:6]
a b c d e f
a 0.00000 15.4287270 12.637979 12.269356 15.666710 12.3919715
b 13.58821 0.0000000 5.356411 13.840444 1.238677 12.6512352
c 12.48161 6.3086852 0.000000 8.427033 6.382304 7.1387840
d 10.69748 13.5936114 7.708183 0.000000 13.718517 0.9836146
e 14.00920 0.7736654 5.980220 14.470826 0.000000 13.2809601
f 12.24503 12.8987043 6.984763 2.182829 12.993283 0.0000000
现在 TSP 包中的三行代码解决了这个问题。
atsp = as.ATSP(distMatrix)
tour = solve_TSP(atsp)
# assume we want to start our circuit at transect "a".
path = cut_tour(tour, "a", exclude_cut = F)
变量 path 显示了您访问横断面的顺序：
> path
a w x q i o l d f s h y g v t k c m e b p u z j r n
1 23 24 17 9 15 12 4 6 19 8 25 7 22 20 11 3 13 5 2 16 21 26 10 18 14
我们可以将路径添加到可视化：
for(i in 1:(length(path)-1)) {
lines(c(e[path[i],1], b[path[i+1],1]), c(e[path[i],2], b[path[i+1], 2]), lty = "dotted")
}
感谢大家的建议，@Simon 的解决方案最适合我的 OP。 @Geoffrey 使用 x,y 坐标的实际方法很棒，因为它允许绘制横断面和旅行顺序。因此，我发布了一个混合解决方案，该解决方案是使用他们两人的代码生成的，以及其他评论和代码来详细说明该过程并获得我想要的实际最终结果。我不确定这是否会在未来帮助任何人，但是，由于没有答案提供 100% 解决我的问题的解决方案，我想我会分享我的想法。
正如其他人所指出的，这是一种旅行推销员问题。它是不对称的，因为从样线终点 "t" 到样线起点 "t+1" 的距离与从终点样线 "t+1" 到样线起点 [=22] 的距离不同=].此外，它是 "path" 解决方案而不是 "cycle" 解决方案。
#=========================================
# Packages
#=========================================
library(TSP)
library(useful)
library(dplyr)
#=========================================
# Full dataset for testing
#=========================================
EndPt <- c(158.7,245.1,187.1,298.2,346.8,317.2,74.5,274.2,153.4,246.7,193.6,302.3,6.8,359.1,235.4,134.5,111.2,240.5,359.2,121.3,224.5,212.6,155.1,353.1,181.7,334,249.3,43.9,38.5,75.7,344.3,45.1,285.7,155.5,183.8,60.6,301,132.1,75.9,112,342.1,302.1,288.1,47.4,331.3,3.4,185.3,62,323.7,188,313.1,171.6,187.6,291.4,19.2,210.3,93.3,24.8,83.1,193.8,112.7,204.3,223.3,210.7,201.2,41.3,79.7,175.4,260.7,279.5,82.4,200.2,254.2,228.9,1.4,299.9,102.7,123.7,172.9,23.2,207.3,320.1,344.6,39.9,223.8,106.6,156.6,45.7,236.3,98.1,337.2,296.1,194,307.1,86.6,65.5,86.6,296.4,94.7,279.9)
StPt <- c(56.3,158.1,82.4,185.5,243.9,195.6,335,167,39.4,151.7,99.8,177.2,246.8,266.1,118.2,358.6,357.9,99.6,209.9,342.8,106.5,86.4,35.7,200.6,65.6,212.5,159.1,297,285.9,300.9,177,245.2,153.1,8.1,76.5,322.4,190.8,35.2,342.6,8.8,244.6,202,176.2,308.3,184.2,267.2,26.6,293.8,167.3,30.5,176,74.3,96.9,186.7,288.2,62.6,331.4,254.7,324.1,73.4,16.4,64,110.9,74.4,69.8,298.8,336.6,58.8,170.1,173.2,330.8,92.6,129.2,124.7,262.3,140.4,321.2,34,79.5,263,66.4,172.8,205.5,288,98.5,335.2,38.7,289.7,112.7,350.7,243.2,185.4,63.9,170.3,326.3,322.9,320.6,199.2,287.1,158.1)
EndID <- c(seq(1, 100, 1))
StID <- c(seq(1, 100, 1))
df <- cbind.data.frame(StPt, StID, EndPt, EndID)
#=========================================
# Convert polar coordinates to cartesian x,y data
#=========================================
# Area that the transect occupy in space only used for graphing
planeDim <- 1
# Number of transects
nTransects <- 100
# Convert 360-degree polar coordinates to x,y cartesian coordinates to facilitate calculating a distance matrix based on the Pythagorean theorem
EndX <- as.matrix(pol2cart(planeDim, EndPt, degrees = TRUE)["x"])
EndY <- as.matrix(pol2cart(planeDim, EndPt, degrees = TRUE)["y"])
StX <- as.matrix(pol2cart(planeDim, StPt, degrees = TRUE)["x"])
StY <- as.matrix(pol2cart(planeDim, StPt, degrees = TRUE)["y"])
# Matrix of x,y pairs for the beginning ("b") and end ("e") points of each transect
b <- cbind(c(StX), c(StY))
e <- cbind(c(EndX), c(EndY))
#=========================================
# Function to calculate the distance from all endpoints in the ePts matrix to a single beginning point in bPt
#=========================================
dist <- function(ePts, bPt) {
# Use the Pythagorean theorem to calculate the hypotenuse (i.e., distance) between every end point in the matrix ePts to the point bPt
apply(ePts, 1, function(p) sum((p - bPt)^2)^0.5)
}
#=========================================
# Distance matrix
#=========================================
# Apply the "dist" function to all beginning points to create a matrix that has the distance between every start and endpoint
## Note: because this is an asymmetric traveling salesperson problem, the distance matrix is directional, thus, the distances at any position in the matrix must be the distance from the transect shown in the row label and to the transect shown in the column label
distMatrix <- apply(b, 1, FUN = dist, ePts = e)
## Set the distance between the beginning and end of each transect to zero so that there is no "cost" to walking the transect
diag(distMatrix) <- 0
#=========================================
# Solve asymmetric TSP
#=========================================
# This creates an instance of the asymmetric traveling salesperson (ASTP)
atsp <- as.ATSP(distMatrix)
# This creates an object of Class Tour that travels to all of the points
## In this case, the repetitive_nn produces the smallest overall and transect-to-transect
tour <- solve_TSP(atsp, method = "repetitive_nn")
#=========================================
# Create a path by cutting the tour at the most "expensive" transition
#=========================================
# Sort the original data frame to match the order of the solution
dfTour = df[as.numeric(tour),]
# Add the following columns to the original dataframe:
dfTour = dfTour %>%
# Assign visit order (1 to 100, ascending)
mutate(visit_order = 1:nrow(dfTour)) %>%
# The ID of the next transect to move to
mutate(next_StID = lead(StID, order_by = visit_order),
# The angle of the start point for the next transect
next_StPt = lead(StPt, order_by = visit_order))
# lead() generates the NA's in the last record for next_StID, next_StPt, replace these by adding that information
dfTour[dfTour$visit_order == nrow(dfTour),'next_StID'] <-
dfTour[dfTour$visit_order == 1,'StID']
dfTour[dfTour$visit_order == nrow(dfTour),'next_StPt'] <-
dfTour[dfTour$visit_order == 1,'StPt']
# Function to calculate distance for 360 degrees rather than x,y coordinates
transect_distance <- function(end,start){
abs_dist = abs(start-end)
ifelse(360-abs_dist > 180, abs_dist, 360-abs_dist)
}
# Compute distance between end of each transect and start of next using polar coordinates
dfTour = dfTour %>% mutate(dist_between = transect_distance(EndPt, next_StPt))
# Identify the longest transition point for breaking the cycle
min_distance <- sum(dfTour$dist_between) - max(dfTour$dist_between)
path_start <- dfTour[dfTour$dist_between == max(dfTour$dist_between), 'next_StID']
path_end <- dfTour[dfTour$dist_between == max(dfTour$dist_between), 'EndID']
# Make a statement about the least cost path
print(sprintf("minimum cost path = %.2f, starting at node %d, ending at node %d",
min_distance, path_start, path_end))
# The variable path shows the order in which you should visit the transects
path <- cut_tour(tour, path_start, exclude_cut = F)
# Arrange df from smallest to largest travel distance
tmp1 <- dfTour %>%
arrange(dist_between)
# Change dist_between and visit_order to NA for transect with the largest distance to break cycle
# (e.g., we will not travel this distance, this represents the path endpoint)
tmp1[length(dfTour$dist_between):length(dfTour$dist_between),8] <- NA
tmp1[length(dfTour$dist_between):length(dfTour$dist_between),5] <- NA
# Set df order back to ascending by visit order
tmp2 <- tmp1 %>%
arrange(visit_order)
# Detect the break in a sequence of visit_order introduced by the NA (e.g., 1,2,3....5,6) and mark groups before the break with 0 and after the break with 1 in the "cont_per" column
tmp2$cont_per <- cumsum(!c(TRUE, diff(tmp2$visit_order)==1))
# Sort "cont_per" such that the records following the break become the beginning of the path and the ones following the break represent the middle orders and the point with the NA being assigned the last visit order, and assign a new visit order
tmp3 <- tmp2%>%
arrange(desc(cont_per))%>%
mutate(visit_order_FINAL=seq(1, length(tmp2$visit_order), 1))
# Datframe ordered by progression of transects
trans_order <- cbind.data.frame(tmp3[2], tmp3[1], tmp3[4], tmp3[3], tmp3[6], tmp3[7], tmp3[8], tmp3[10])
# Insert NAs for "next" info for final transect
trans_order[nrow(trans_order),'next_StPt'] <- NA
trans_order[nrow(trans_order), 'next_StID'] <- NA
#=========================================
# View data
#=========================================
head(trans_order)
#=========================================
# Plot
#=========================================
#For fun, we can visualize the transects:
# make an empty graph space
plot(1,1, xlim = c(-planeDim-0.1, planeDim+0.1), ylim = c(-planeDim-0.1, planeDim+0.1), ty = "n")
# plot the beginning of each transect as a green point, the end as a red point,
and a grey line representing the transect
for(i in 1:nrow(e)) {
xs = c(b[i,1], e[i,1])
ys = c(b[i,2], e[i,2])
lines(xs, ys, col = "light grey", lwd = 1, lty = 1)
points(xs, ys, col = c("green", "red"), pch = 1, cex = 1)
#text((xs), (ys), i)
}
# Add the path to the visualization
for(i in 1:(length(path)-1)) {
# This makes a line between the x coordinates for the end point of path i and beginning point of path i+1
lines(c(e[path[i],1], b[path[i+1],1]), c(e[path[i],2], b[path[i+1], 2]), lty = 1, lwd=1)
}
这就是最终结果的样子