又一个WordPress站点

忆江南三首古诗番外篇(2):用R语言规划旅行的第一次尝试-老宅不宅

番外篇(2):用R语言规划旅行的第一次尝试-老宅不宅
上个《番外篇》虽然以失败告终,不过老宅干劲满满,又开始准备新的挑战。还记得上回结尾时说过的吗?老宅的终极目标是路线规划。所谓“人有多大胆,地有多大产”,老宅这不就摩拳擦掌,跃跃欲试了吗?其实最初的动机是因为宅妈下周就从老宅身边离开回国了,老宅想给宅妈作个汇报,展示一下最近老宅的学习成果,让宅妈高兴高兴。走之前再用AmEx白金卡给的Uber额度请宅妈吃一顿免费大餐,让宅妈的美国之旅有个完满的结局。正文开始之前,大家先猜猜,这次老宅成功了吗?前期的数据准备
先说说老宅的计划:
1月18号,乘坐AS769航班,晚上5:30从波士顿直飞圣地亚哥,当地时间9:10到达,入住圣地亚哥市中心的圣地亚哥共和国,奥途格精选;
1月19号玩一天以后,入住五星级的格兰特将军,豪华精选(这家酒店不仅是万豪在圣地亚哥的头牌,而且AmEx有375-75刀的offer,AmEx FHR更有100刀的酒店消费额度);
1月20号再玩一天,晚上入住圣地亚哥机场万怡酒店(这家酒店不仅离机场近,适合最后一天早上的航班,AmEx还有200-40刀的offer);
1月21号早上7:50乘坐AS798直飞波士顿。
在这两天时间里,老宅打算玩的地方有:巴波亚公园(半天)孙菲菲被打 ,中途岛号航空母舰博物馆(2小时),科罗纳多岛(半天),卡波里奥国家公园(2小时),拉霍亚海湾(2小时),圣地亚哥老城(半天)。假如老宅打算玩A、B和C三个地方,因为在这三个地方的时间是一定的,那么要找的就是花最少时间在三个地方之间转移的路径。比如B->A->C花的时间最少,这条路径就是老宅要的。因为每个地方的游览时间涉及到下面的旅行方案的计算(游览时间参考了穷游的时间推荐),于是先制作了下面的表。
> place
# A tibble: 10 x 5
namehour type short abbr
<chr> <dbl> <chr> <chr> <chr>
1 Balboa Park, san diego, ca4 place Balboa Park bal
2 USS Midway,李思晓 san diego, ca2 place USS Midwaymidway
3 Coronado, san diego, ca 4 place Coronado coro
4 Cabrillo National Monument, san diego, ca 2 place Cabrillo cab
5 La Jolla Cove, san diego, ca 2 place La Jolla ljc
6 4002 Wallace St, San Diego, CA 92110 4 place Old Town ot
7 Hotel Republic San Diego, Autograph Collection NA hotel Autograph Collection NA
8 THE US GRANT, a Luxury Collection Hotel, San Diego NA hotel Luxury Collection NA
9 2592 Laning Road, San Diego, California 92106 NA hotel CourtyardNA
10 san diego airport NA airport airport NA 绘制坐标地图
上文说过,老宅试图在地图上标记地点而未成功。现在老宅知道如何做了!而且画的图还非常漂亮。这都要归功于 ggrepel包!
library(ggmap)
register_google(key = "my_google_key_here")
library(ggrepel)
geoc <- geocode(place$name)
place <- bind_cols(place, geoc)
ggmap(get_map(location = 'san diego', zoom = 11)) +
geom_label_repel(data = place, aes(label = short, color = type), force = 5, nudge_x = 1) +
geom_point(data = place, aes(lon, lat, color = type))
#https://cran.r-project.org/web/packages/ggrepel/vignettes/ggrepel.html

还不错吧?核心算法(一)
接下来是关键:算法的设计。
老宅的设计是:先把所有可能的排列组合找出来( 存在plan里),然后按照已经计划好的游览时间去除每天超过10小时或不足8小时的组合,剩下的组合留着下一步的计算。这个算法的核心是 gtools包中的 permutations函数穆雷桑 ,它可以找出一个数组的所有可能排列组合。
library(gtools)
plan <- permutations(2, 6, v=c(1, 2), repeats.allowed = T)
得到的前六种排列组合是这样子的:1和2分别代表在第一天或第二天游览,day1和 day2列就是这些组合下第一天和第二天分别需要的时间。组合一共有64个。再去除时间不符合安排(大于10小时或小于8小时)的组合:
plan <- plan %>%
filter(day1 >= 8 & day1 <= 10 & day2 >= 8 & day2 <= 10)
最后剩下24个组合。
> head(plan)
bal midway coro cab ljc ot day1 day2
1 11 1 2 2 2 10 8
2 11 2 1 1 2 10 8
3 11 2 1 2 2 8 10
4 11 2 2 1 2 8 10
5 11 2 2 2 1 10 8
6 12 1 1 2 2 10 8
> nrow(plan)
[1] 24核心算法(二)
上面得到的仅仅是每天游览地点的分组,下一步是要得到每个分组的每天内游玩顺序的排列组合。这里再一次用到了 permutations函数。因为不同组合中第一天、第二天游玩景点的数量不同,所以这里用了一个for循环来计算。这里需要先把1和2用真实地面替换,以便下一步用谷歌地图提取交通时间。老宅先把 paln里等于1和2的地方分别提取了出来(1和2用地名缩写代替)安泽熙 ,然后用新学到的 case_when函数,把排列组合中的地名缩写用真实名称替换。
day_1 <- list()
day_2 <- list()
for (i in 1:nrow(plan)) {
j <- sum(plan[i, ] == 1)
day_1[[i]] <- data.frame(permutations(j, j,
v = names(plan)[which(plan[i, ] == 1, arr.ind=T)[, "col"]]), stringsAsFactors = FALSE)
#https://stackoverflow.com/questions/36960010/get-column-name-that-matches-specific-row-value-in-dataframe
k <- sum(plan[i, ] == 2)
day_2[[i]] <- 同上,略
}
#https://stackoverflow.com/questions/53625126/how-to-create-data-frames-not-just-one-at-once-in-r
for (i in 1:length(day_1)) {
for (j in 1:length(day_1[[i]])) {
for (k in 1:length(day_1[[i]][,j])) {
day_1[[i]][, j][k] <- case_when(day_1[[i]][, j][k] == "bal" ~ "Balboa Park, san diego, ca",
day_1[[i]][, j][k] == "midway" ~ "USS Midway, san diego, ca",
day_1[[i]][, j][k] == "coro" ~ "Coronado, san diego, ca",
day_1[[i]][, j][k] == "cab" ~ "Cabrillo National Monument, san diego, ca",
day_1[[i]][, j][k] == "ljc" ~ "La Jolla Cove, san diego, ca",
day_1[[i]][, j][k] == "ot" ~ "4002 Wallace St, San Diego, CA 92110")
}
}
}
for (i in 1:length(day_2)) {
同上,略
}
因为这是组合中的组合,所以新生成的组合没法保存在组合中。这个问题困扰了老宅好久。在Stack Overflow上一位的网友建议下,把新的组合保存在列表中。看看第一天的可能游览方案:
> day_1[[1]]
X1 X2 X3
1 Balboa Park, san diego, ca Coronado, san diego, ca USS Midway, san diego, ca
2 Balboa Park, san diego军婚娇妻撩人, ca USS Midway鬼城凶梦 , san diego, ca Coronado, san diego, ca
3 Coronado, san diego, ca Balboa Park, san diego, ca USS Midway, san diego, ca
4 Coronado, san diego, ca USS Midway, san diego, ca Balboa Park, san diego, ca
5 USS Midway, san diego, ca Balboa Park, san diego, ca Coronado, san diego, ca
6 USS Midway, san diego, ca Coronado, san diego, ca Balboa Park, san diego, ca核心算法(三)
有了所有的排列组合,下一步就是利用谷歌地图API提取每个组合中的不同地点之间的交通时间之和,耗时最短的一个组合最佳。然后比较所有组合的耗时,取最少的一组。这里老宅遇到了一个大麻烦,原本老宅设计了一个多层嵌套的循环来提取列表中的矩阵中的耗时,但是运行的时候有错误无法运行,老宅又找不到错误在哪。老宅把这个问题又一次发到了Stack Overflow上。一位网友写了一个功能,成功解决了这个问题。这个功能经过老宅的修改以后,是下面的样子:
start <- place$name[7]
end1 <- place$name[8]
end2 <- place$name[9]
compute_route_time = function(row张咏琪 , start, end) {
# get a vector of stops that you can loop through
stops = unlist(row)
# distance from start to first stop
time = mapdist(from = start, to = stops[1])$seconds
# iterate through the rest of the route
n_stops = length(stops)
for (i in 1:(n_stops-1)) {
time = time + mapdist(from = stops[i], to = stops[i+1])$seconds
}
# finish rout
time = time + mapdist(from = stops[n_stops], to = end)$seconds
# collect row as a data frame to preserve old columns
row = as.data.frame(row)
# add time column
row$time = time
return(row)
}
for (i in 1:length(day_1)) {
day_1[[i]] <- day_1[[i]] %>%
rowwise() %>%
do(compute_route_time(., start = start, end = end1))
}
for (i in 1:length(day_2)) {
day_2[[i]] <- day_2[[i]] %>%
rowwise() %>%
do(compute_route_time(., start = end1, end = end2))
}
#https://stackoverflow.com/questions/53662127/non-numeric-argument-to-binary-operator-in-ggmap
为了防止服务器过载,ggmap人为设置了从谷歌地图上请求数据的间隔,速度非常慢。读取完毕后,来看看 plan里第一个组合的耗时情况(以秒为单位),可以看到第6种路径耗时最少:
# A tibble: 6 x 4
X1 X2 X3 time
* <chr> <chr> <chr> <dbl>
1 Balboa Park, san diego, ca Coronado, san diego, ca USS Midway, san diego, ca 2674
2 Balboa Park, san diego, ca USS Midway, san diego, ca Coronado, san diego, ca2906
3 Coronado, san diego, ca Balboa Park, san diego, ca USS Midway, san diego, ca 2608
4 Coronado, san diego, ca USS Midway, san diego, ca Balboa Park, san diego, ca 2931
5 USS Midway, san diego, ca Balboa Park, san diego, ca Coronado, san diego, ca2621
6 USS Midwayg5308w , san diego, ca Coronado, san diego, ca Balboa Park, san diego, ca 2458
然后把每一种组合中第一天、第二天的最小路上耗时相加,得到的既是该组合中的最佳方案。
route_1 <- list()
route_2 <- list()
total_time <- NULL
for (i in 1:length(day_1)) {
route_1[[i]] <- day_1[[i]][which.min(day_1[[i]]$time),]
route_2[[i]] <- day_2[[i]][which.min(day_2[[i]]$time),]
total_time <- c(total_time, route_1[[i]]$time + route_2[[i]]$time)
}
再比较各个最佳方案,选取总的最佳方案。1号方案最终胜出。
> (max <- which.min(total_time))
[1] 1制作路线图
有了最佳路线,下面就是要把最佳路线重新做成一个行程表格,以便下面的路线图绘制希斯堡惨案。以第一天为例:
plan_1 <- data.frame(matrix(ncol = 2, nrow = length(route_1[[max]]) + 1))
colnames(plan_1) <- c("id", "location")
for (i in 1:(length(route_1[[max]]) + 1)) {
if(i == 1){
plan_1$id[1] <- 1
plan_1$location[1] <- "Hotel Republic San Diego, Autograph Collection"
}
if(i > 1 & i < length(route_1[[max]]) + 1){
plan_1$id[i] <- i
plan_1$location[i] <- route_1[[max]][i - 1]
}
if(i == length(route_1[[max]]) + 1){
plan_1$id[i] <- length(route_1[[max]]) + 1
plan_1$location[i] <- "THE US GRANT, a Luxury Collection Hotel, San Diego"
}
}
再一次以第一天的行程为例,新的表格看起来如下,id代表了到达顺序:
> plan_1
id location
1 1Hotel Republic San Diego, Autograph Collection
2 2 USS Midway, san diego, ca
3 3 Coronado, san diego, ca
4 4 Balboa Park, san diego, ca
5 5 THE US GRANT, a Luxury Collection Hotel, San Diego
有了行程表格,然后就是利用谷歌地图API提取路线了:
legs_df_1 <- list()
for (i in 1:(nrow(plan_1) - 1)) {
legs_df_1[[i]] <- route(from = as.character(plan_1$location[i]), to = as.character(plan_1$location[i + 1]), alternatives = FALSE, mode = "driving", structure = "leg")
legs_df_1[[i]]$id <- i
}
for (i in 1:length(legs_df_1)) {
if(i == 1){
final_1 <- legs_df_1[[1]][which.min(legs_df_1[[1]]$seconds),]
}
if(i != 1){
final_1 <- bind_rows(final_1, legs_df_1[[i]][which.min(legs_df_1[[i]]$seconds),])
}
}
然后用这个路线信息绘制路线图:
ggmap(get_map(location = 'san diego', zoom = 13)) +
geom_label_repel(data = place, aes(label = short, color = type), force = 5, nudge_x = 1) +
geom_point(data = place, aes(lon, lat, color = type)) +
geom_leg(aes(x = startLon, y = startLat, xend = endLon, yend = endLat),
alpha = 3/4, size = 2, data = legs_df_1[[1]], color = "red") +
geom_leg(aes(x = startLon, y = startLat, xend = endLon, yend = endLat),
alpha = 3/4, size = 2, data = legs_df_1[[2]]第一等傻女, color = "blue") +
geom_leg(aes(x = startLon, y = startLat, xend = endLon, yend = endLat),
alpha = 3/4, size = 2, data = legs_df_1[[3]], color = "green") +
geom_leg(aes(x = startLon, y = startLat, xend = endLon, yend = endLat),
alpha = 3/4, size = 2, data = legs_df_1[[4]], color = "yellow") +
geom_label_repel(aes(x = startLon, y = startLat), data = final_1, label = final_1$id) +
ggtitle("Day 1") +
theme(plot.title = element_text(colour = "red", size = 20, face = "bold"))
这是最后的结果(谷歌地图提取坐标还是有问题):


第一天从圣地亚哥共和国宾馆出发,依次游览中途岛号航空母舰博物馆、科罗纳多岛、巴波亚公园,晚上回到格兰特将军宾馆,结束第一天的行程;
第二天从格兰特将军宾馆出发,依次游览拉霍亚海湾、圣地亚哥老城、卡波里奥国家公园,晚上回到万怡宾馆,结束第二天的行程。
后记与感想
这次老宅的设计,至少在理论上,取得了圆满成功。不过这个设计有一个缺陷,没有考虑到中午吃饭的问题:假如在某一天的第一个景点呆两小时,第二个景点呆四小时,那么在第二个景点玩到一半的时候就饿了。不过这在技术上也不是问题,只要把2个整天分为4个半天就可以了。所以这次的设计仅供参考。老宅懒,不想再修改代码了,所以这次就这样吧忆江南三首古诗。
通过这个项目,老宅的收获非常大。在完成这个项目的过程中,老宅遇到了无数问题。有些问题老宅自己解决了,有些问题要到网上查资料、提问题璀璨王座 。这里老宅强烈安利一个网站:Stack Overflow。这是一个编程问答社区,有问题到这里提问,总会有高手回答的。老宅一共写了260行代码,以前从来没有在一个项目里写过这么多代码。当然,在老手眼里260行不算什么,但对老宅来说是个突破。以后回头看,肯定会觉得现在的算法逻辑、代码太稚嫩。不过什么事都有个过程。
本文中的完整代码可以在https://github.com/vincent507cpu/r_practice/blob/master/plan_to_san_diego.r上看到威尔惠顿 。
而且这次的公众号文章的写作让老宅对Markdown这个文本编辑语言有了初步了解,并且亲自实战体验了一把。可以说,这是这个项目的最大收获之一。
虽然这次的实践对这一次的旅游仅有指导意义,但是老宅计划三月去纽约玩,预计会入住几家顶级豪华酒店(每晚1000+的房费肯定不是老宅出啦),去几家米其林餐厅用餐。希望到时候,这次的技术储备能派上大用场。
老宅的职业目标是成为一名数据科学家。不过除了数据科学家,还有数据分析师这个职位。老宅不想成为数据分析师,不仅仅因为数据分析师的薪水远没有数据科学家高,还因为两者的工作其实有区别:数据分析师是从既有数据中得出结论,而数据科学家是用数据回答新问题。换句话说,数据科学家创造了新东西。老宅希望能够有所创造。从最近的两次尝试来看,第一次尝试趋向数据分析师的工作,而第二次尝试趋向数据科学家的工作。所以老宅对第二次的尝试尤其满意。
不说了,一月份圣地亚哥见!
如果你对小钱大玩的奢华旅行和最大化消费回报感兴趣,欢迎订阅我的公众号“老宅不宅”!
每周日更新,美国信用卡老司机带你上车!

关键词: