基于运输时间的热图/轮廓（反向等时轮廓）
Posted
技术标签:
【中文标题】基于运输时间的热图/轮廓（反向等时轮廓）【英文标题】：Heatmap/Contours based on Transportation Time (Reverse Isochronic Contours) 【发布时间】：2018-12-05 11:03:10 【问题描述】：注意： 需要r、python、java 或c++ 或c# 中的解决方案。
我正在尝试根据运输时间绘制轮廓。更清楚地说，我想将具有相似旅行时间（比如说 10 分钟间隔）的点聚集到特定点（目的地），并将它们映射为等高线或热图。
现在，我唯一的想法是使用 R 包 gmapsdistance 来查找不同来源的旅行时间，然后将它们聚类并绘制在地图。但是，正如您所知，这绝不是一个可靠的解决方案。
GIS 社区上的 thread 和 python 上的 this one 说明了类似的问题，但在特定时间到达目的地的起点。我想找到可以在一定时间内到达目的地的起点。
现在，下面的代码展示了我的基本想法（使用 R）：
library(gmapsdistance)
set.api.key("YOUR.API.KEY")
mdestination <- "40.7+-73"
morigin1 <- "40.6+-74.2"
morigin2 <- "40+-74"
gmapsdistance(origin = morigin1,
destination = mdestination,
mode = "transit")
gmapsdistance(origin = morigin2,
destination = mdestination,
mode = "transit")
这张地图也可能有助于理解这个问题：
使用这个 answer 我可以得到我可以从原点到达的点，但我需要反转它并找到有旅行的点到我的目的地的时间小于一定时间；
library(httr)
library(googleway)
library(jsonlite)
appId <- "TravelTime_APP_ID"
apiKey <- "TravelTime_API_KEY"
mapKey <- "GOOGLE_MAPS_API_KEY"
location <- c(40, -73)
CommuteTime <- (5 / 6) * 60 * 60
url <- "http://api.traveltimeapp.com/v4/time-map"
requestBody <- paste0('
"departure_searches" : [
"id" : "test",
"coords": "lat":', location[1], ', "lng":', location[2],' ,
"transportation" : "type" : "driving" ,
"travel_time" : ', CommuteTime, ',
"departure_time" : "2017-05-03T07:20:00z"
]
')
res <- httr::POST(url = url,
httr::add_headers('Content-Type' = 'application/json'),
httr::add_headers('Accept' = 'application/json'),
httr::add_headers('X-Application-Id' = appId),
httr::add_headers('X-Api-Key' = apiKey),
body = requestBody,
encode = "json")
res <- jsonlite::fromJSON(as.character(res))
pl <- lapply(res$results$shapes[[1]]$shell, function(x)
googleway::encode_pl(lat = x[['lat']], lon = x[['lng']])
)
df <- data.frame(polyline = unlist(pl))
df_marker <- data.frame(lat = location[1], lon = location[2])
google_map(key = mapKey) %>%
add_markers(data = df_marker) %>%
add_polylines(data = df, polyline = "polyline")
此外，Documentation of Travel Time Map Platform 谈到 Multi Origins with Arrival time，这正是我想做的事情。但是我需要为公共交通和驾驶（对于通勤时间不到一小时的地方）这样做，我认为由于公共交通很棘手（基于您靠近的车站），也许热图是比等高线更好的选择。
【问题讨论】：
看这里：***.com/questions/40489162/… @NicolásVelásquez 我看了看；它正在回答一个问题，即找到您可以在特定时间到达的地方，而不是您可以开始并在所需时间到达特定目的地的地方。 （您可能对制作等时地图的mapumental.com 感兴趣。） @Rich 看来他们目前只覆盖英国。 mapumental.com/about 谢谢。 【参考方案1】：此答案基于获得（大致）等距点的网格之间的起点-终点矩阵。这是一项计算机密集型操作，不仅因为它需要对地图服务进行大量 API 调用，还因为服务器必须为每次调用计算一个矩阵。所需调用的数量随着网格中的点数呈指数增长。
要解决这个问题，我建议您考虑在本地计算机或本地服务器上运行映射服务器。 Project OSRM 提供了一个相对简单、免费和开源的解决方案，使您能够在 Linux docker (https://github.com/Project-OSRM/osrm-backend) 中运行 OpenStreetMap 服务器。拥有自己的本地映射服务器将允许您进行任意数量的 API 调用。 R 的 osrm 包允许您与 OpenStreetMaps 的 API 进行交互，包括那些放置到本地服务器的 API。
library(raster) # Optional
library(sp)
library(ggmap)
library(tidyverse)
library(osrm)
devtools::install_github("cmartin/ggConvexHull") # Needed to quickly draw the contours
library(ggConvexHull)
我在布鲁塞尔（比利时）大都市圈周围创建了一个由 96 个大致相等距离的点组成的网格。 这个网格没有考虑地球的曲率，在城市距离的层面上可以忽略不计。
为方便起见，我使用 raster 包下载比利时的 ShapeFile 并提取布鲁塞尔市的节点。
BE <- raster::getData("GADM", country = "BEL", level = 1)
Bruxelles <- BE[BE$NAME_1 == "Bruxelles", ]
df_grid <- makegrid(Bruxelles, cellsize = 0.02) %>%
SpatialPoints() %>%
## I convert the SpatialPoints object into a simple data.frame
as.data.frame() %>%
## create a unique id for each point in the data.frame
rownames_to_column() %>%
## rename variables of the data.frame with more explanatory names.
rename(id = rowname, lat = x2, lon = x1)
## I point osrm.server to the OpenStreet docker running in my Linux machine. ...
### ... Do not run this if you are getting your data from OpenStreet public servers.
options(osrm.server = "http://127.0.0.1:5000/")
## I obtain a list with distances (Origin Destination Matrix in ...
### ... minutes, origins and destinations)
Distance_Tables <- osrmTable(loc = df_grid)
OD_Matrix <- Distance_Tables$durations %>% ## subset the previous list
## convert the Origin Destination Matrix into a tibble
as_data_frame() %>%
rownames_to_column() %>%
## make sure we have an id column for the OD tibble
rename(origin_id = rowname) %>%
## transform the tibble into long/tidy format
gather(key = destination_id, value = distance_time, -origin_id) %>%
left_join(df_grid, by = c("origin_id" = "id")) %>%
## set origin coordinates
rename(origin_lon = lon, origin_lat = lat) %>%
left_join(df_grid, by = c("destination_id" = "id")) %>%
## set destination coordinates
rename(destination_lat = lat, destination_lon = lon)
## Obtain a nice looking road map of Brussels
Brux_map <- get_map(location = "bruxelles, belgique",
zoom = 11,
source = "google",
maptype = "roadmap")
ggmap(Brux_map) +
geom_point(aes(x = origin_lon, y = origin_lat),
data = OD_Matrix %>%
## Here I selected point_id 42 as the desired target, ...
## ... just because it is not far from the City Center.
filter(destination_id == 42),
size = 0.5) +
## Draw a diamond around point_id 42
geom_point(aes(x = origin_lon, y = origin_lat),
data = OD_Matrix %>%
filter(destination_id == 42, origin_id == 42),
shape = 5, size = 3) +
## Countour marking a distance of up to 8 minutes
geom_convexhull(alpha = 0.2,
fill = "blue",
colour = "blue",
data = OD_Matrix %>%
filter(destination_id == 42,
distance_time <= 8),
aes(x = origin_lon, y = origin_lat)) +
## Countour marking a distance of up to 16 minutes
geom_convexhull(alpha = 0.2,
fill = "red",
colour = "red",
data = OD_Matrix %>%
filter(destination_id == 42,
distance_time <= 15),
aes(x = origin_lon, y = origin_lat))
结果
蓝色等高线表示到市中心最多 8 分钟的距离。 红色轮廓代表最长 15 分钟的距离。
【讨论】：
好吧，如果计算成本是关键，而且这对您来说是一个非常重要的项目，您可能会使用其他工具。如果您有道路和公共交通网络的多边形形状文件，以及用于计算速度的模型或每种道路/交通模式的速度值数据，您可以使用 ArcGis 的 Network Analisys 工具，或者通过一些额外的编码， Qgis 的路线图。您需要将代表道路的每条线分成非常小的段并计算所有段之间的 OD 矩阵。不过，我不知道 R 如何在 R 中做到这一点。 感谢 Nicolas 的想法。通过结合一些想法（主要是您的答案和您指出的主题），我能够获得预期的结果。澄清一下，你所说的关于 GIS 软件的内容在 R 中也是可行的。SymbolixAU 在开发人员模式下有几个包，可以下载公共交通时间表，你也可以获得道路网络。但是，仍然存在巨大的内存/处理问题。我认为我在回答中的建议（从等时线的交点中挑选一些点）会让我们足够接近。比这更接近，需要更多的计算能力。干杯-M
太棒了！很高兴能提供帮助，这是一个非常有趣的问题。祝项目顺利！【参考方案2】：
我想出了一种与进行大量 api 调用相比更适用的方法。
这个想法是找到你可以在特定时间到达的地方（看看这个thread）。可以通过将时间从早上更改为晚上来模拟交通。你最终会得到一个重叠的区域，你可以从两个地方到达。
然后您可以使用Nicolas answer 并在该重叠区域内绘制一些点，并为您拥有的目的地绘制热图。这样，您将有更少的区域（点）需要覆盖，因此您将进行更少的 api 调用（请记住为此使用适当的时间）。
在下面，我试图证明我所说的这些的意思，并让您明白您可以制作另一个答案中提到的网格，以使您的估计更加稳健。
这显示了如何绘制相交区域。
library(httr)
library(googleway)
library(jsonlite)
appId <- "Travel.Time.ID"
apiKey <- "Travel.Time.API"
mapKey <- "Google.Map.ID"
locationK <- c(40, -73) #K
locationM <- c(40, -74) #M
CommuteTimeK <- (3 / 4) * 60 * 60
CommuteTimeM <- (0.55) * 60 * 60
url <- "http://api.traveltimeapp.com/v4/time-map"
requestBodyK <- paste0('
"departure_searches" : [
"id" : "test",
"coords": "lat":', locationK[1], ', "lng":', locationK[2],' ,
"transportation" : "type" : "public_transport" ,
"travel_time" : ', CommuteTimeK, ',
"departure_time" : "2018-06-27T13:00:00z"
]
')
requestBodyM <- paste0('
"departure_searches" : [
"id" : "test",
"coords": "lat":', locationM[1], ', "lng":', locationM[2],' ,
"transportation" : "type" : "driving" ,
"travel_time" : ', CommuteTimeM, ',
"departure_time" : "2018-06-27T13:00:00z"
]
')
resKi <- httr::POST(url = url,
httr::add_headers('Content-Type' = 'application/json'),
httr::add_headers('Accept' = 'application/json'),
httr::add_headers('X-Application-Id' = appId),
httr::add_headers('X-Api-Key' = apiKey),
body = requestBodyK,
encode = "json")
resMi <- httr::POST(url = url,
httr::add_headers('Content-Type' = 'application/json'),
httr::add_headers('Accept' = 'application/json'),
httr::add_headers('X-Application-Id' = appId),
httr::add_headers('X-Api-Key' = apiKey),
body = requestBodyM,
encode = "json")
resK <- jsonlite::fromJSON(as.character(resKi))
resM <- jsonlite::fromJSON(as.character(resMi))
plK <- lapply(resK$results$shapes[[1]]$shell, function(x)
googleway::encode_pl(lat = x[['lat']], lon = x[['lng']])
)
plM <- lapply(resM$results$shapes[[1]]$shell, function(x)
googleway::encode_pl(lat = x[['lat']], lon = x[['lng']])
)
dfK <- data.frame(polyline = unlist(plK))
dfM <- data.frame(polyline = unlist(plM))
df_markerK <- data.frame(lat = locationK[1], lon = locationK[2], colour = "#green")
df_markerM <- data.frame(lat = locationM[1], lon = locationM[2], colour = "#lavender")
iconK <- "red"
df_markerK$icon <- iconK
iconM <- "blue"
df_markerM$icon <- iconM
google_map(key = mapKey) %>%
add_markers(data = df_markerK,
lat = "lat", lon = "lon",colour = "icon",
mouse_over = "K_K") %>%
add_markers(data = df_markerM,
lat = "lat", lon = "lon", colour = "icon",
mouse_over = "M_M") %>%
add_polygons(data = dfM, polyline = "polyline", stroke_colour = '#461B7E',
fill_colour = '#461B7E', fill_opacity = 0.6) %>%
add_polygons(data = dfK, polyline = "polyline",
stroke_colour = '#F70D1A',
fill_colour = '#FF2400', fill_opacity = 0.4)
你可以像这样提取相交区域：
# install.packages(c("rgdal", "sp", "raster","rgeos","maptools"))
library(rgdal)
library(sp)
library(raster)
library(rgeos)
library(maptools)
Kdata <- resK$results$shapes[[1]]$shell
Mdata <- resM$results$shapes[[1]]$shell
xyfunc <- function(mydf)
xy <- mydf[,c(2,1)]
return(xy)
spdf <- function(xy, mydf)
sp::SpatialPointsDataFrame(
coords = xy, data = mydf,
proj4string = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"))
for (i in (1:length(Kdata))) Kdata[[i]] <- xyfunc(Kdata[[i]])
for (i in (1:length(Mdata))) Mdata[[i]] <- xyfunc(Mdata[[i]])
Kshp <- list(); for (i in (1:length(Kdata))) Kshp[i] <- spdf(Kdata[[i]],Kdata[[i]])
Mshp <- list(); for (i in (1:length(Mdata))) Mshp[i] <- spdf(Mdata[[i]],Mdata[[i]])
Kbind <- do.call(bind, Kshp)
Mbind <- do.call(bind, Mshp)
#plot(Kbind);plot(Mbind)
x <- intersect(Kbind,Mbind)
#plot(x)
xdf <- data.frame(x)
xdf$icon <- "https://i.stack.imgur.com/z7NnE.png"
google_map(key = mapKey,
location = c(mean(latmax,latmin), mean(lngmax,lngmin)), zoom = 8) %>%
add_markers(data = xdf, lat = "lat", lon = "lng", marker_icon = "icon")
这只是相交区域的示意图。
现在，您可以从 xdf 数据框获取坐标并围绕这些点构建您的网格，最终得出热图。为了尊重提出该想法/答案的其他用户，我没有将其包含在我的内容中，而只是在引用它。
Nicolás Velásquez - Obtaining an Origin-Destination Matrix between a Grid of (Roughly) Equally Distant Points
【讨论】：
以上是关于基于运输时间的热图/轮廓（反向等时轮廓）的主要内容，如果未能解决你的问题，请参考以下文章
ML之kmeans：通过数据预处理(分布图箱线图热图/文本转数字/构造特征/编码/PCA)利用kmeans实现汽车产品聚类分析(SSE-平均轮廓系数图/聚类三维图/雷达图/饼图柱形图)/竞品分析之详细
PSO基于PSO粒子群优化的物料点货物运输成本最低值计算matlab仿真，包括运输费用代理人转换费用运输方式转化费用和时间惩罚费用