合并并绘制多个等时线

Merge and plot multiple isochrones

我喜欢做什么

我喜欢在地图上绘制多个位置的等时线,这样我就可以直观地找到从任意城镇到最近位置的旅行时间。它应该看起来像核密度二维图:

library(purrr)
library(ggmap)

locations <- tibble::tribble(
  ~city,     ~lon,      ~lat,
  "Hamburg",  9.992246, 53.550354,
  "Berlin",  13.408163, 52.518527,
  "Rostock", 12.140776, 54.088581
)

data <- map2_dfr(locations$lon, locations$lat, ~ data.frame(lon = rnorm(10000, .x, 0.8),
                                                            lat = rnorm(10000, .y, 0.7)))
ger <- c(left = min(locations$lon) - 1,  bottom = min(locations$lat) - 1, 
         right = max(locations$lon) + 1, top = max(locations$lat) + 1)
get_stamenmap(ger, zoom = 7, maptype = "toner-lite") %>% 
  ggmap() + 
  stat_density_2d(data = data, aes(x= lon, y = lat, fill = ..level.., alpha = ..level..), 
                  geom = "polygon") +
  scale_fill_distiller(palette = "Blues", direction = 1, guide = FALSE) +
  scale_alpha_continuous(range = c(0.1,0.3), guide = FALSE)

我试过的

您可以通过 osrm 轻松获取等时线并使用 leaflet 绘制它们。然而,这些等时线彼此独立。当我绘制它们时,它们相互重叠。

library(osrm)
library(leaflet)
library(purrr)
library(ggmap)

locations <- tibble::tribble(
  ~city,     ~lon,      ~lat,
  "Hamburg",  9.992246, 53.550354,
  "Berlin",  13.408163, 52.518527,
  "Rostock", 12.140776, 54.088581
)


isochrone <- map2(locations$lon, locations$lat, 
                  ~ osrmIsochrone(loc = c(.x, .y),
                                  breaks = seq(0, 120, 30))) %>%
  do.call(what = rbind)

isochrone@data$drive_times <- factor(paste(isochrone@data$min, "bis", 
                                           isochrone@data$max, "Minuten"))

factpal <- colorFactor("Blues", isochrone@data$drive_times, reverse = TRUE)

leaflet() %>% 
  setView(mean(locations$lon), mean(locations$lat), zoom = 7) %>%
  addProviderTiles("Stamen.TonerLite") %>%
  addPolygons(fill = TRUE, stroke = TRUE, color = "black",
              fillColor = ~factpal(isochrone@data$drive_times),
              weight = 0.5, fillOpacity = 0.6,
              data = isochrone, popup = isochrone@data$drive_times,
              group = "Drive Time") %>% 
  addLegend("bottomright", pal = factpal, values = isochrone@data$drive_time,   
            title = "Fahrtzeit")

如何合并这些等时线以使它们不重叠?

非常酷的问题。您要做的是按 ID 合并形状,因此所有 0-30 分钟区域都是一个形状,所有 30-60 分钟区域都是另一个形状,依此类推。有一些方法可以用其他空间包来做到这一点,但似乎 well-suited 到 sf,它使用 dplyr 风格的函数。

创建isochrone后,可以将其转换为sf对象,制作同类型的距离标签,按ID分组,调用summarise。当您汇总 sf 个对象时,默认情况下只是一个空间联合,因此您不需要在那里提供函数。

library(sf)
library(dplyr)

iso_sf <- st_as_sf(isochrone)

iso_union <- iso_sf %>%
  mutate(label = paste(min, max, sep = "-")) %>%
  group_by(id, label) %>%
  summarise()

我没有 leaflet 方便,所以这里只是默认打印方法:

plot(iso_union["label"], pal = RColorBrewer::brewer.pal(4, "Blues"))

我不确定具有突然垂直边缘的区域是怎么回事,但它们也在您的绘图中。

我很难使用您使用的 map2 方法,因为它既执行并集,又执行另一个集合论,例如创建特定区间的函数。相反,我会建议创建一个栅格层您创建的图层并将一个不透明度应用于该栅格,就像 ggmap 示例所做的那样。有一个很棒的博客 post,我从 here(以及 user:camille)窃取了很多代码。

它使用不同的 API 需要 mapbox 但它是免费的。另一个限制是它不会 return 等长线是您喜欢的大小,但我在另一个位置重新创建了它,其中三个点靠得更近以证明该方法。

我也没有费心去矢量化创建 isocrone 网络请求的过程,所以我把它留给了更聪明的人。

# First be sure to get your mapbox token

library(fasterize)
library(sf)
library(mapboxapi)
library(leaflet)
#mapboxapi::mb_access_token("Go get the token and put it here",
#                           install = TRUE, overwrite = TRUE)

isos1 <- mb_isochrone(
  location = c("-149.883234, 61.185765"),
  profile = "driving",
  time = c(5,10,15),
)

isos2 <- mb_isochrone(
  location = c("-149.928200, 61.191227"),
  profile = "driving",
  time = c(5,10,15),
)
isos3 <- mb_isochrone(
  location = c("-149.939484, 61.160192"),
  profile = "driving",
  time = c(5,10,15),
)

library(sf)
library(dplyr)

isocrones <- rbind(isos1,isos2,isos3)

iso_sf <- st_as_sf(isocrones)

iso_union <- iso_sf %>%
  group_by(time) %>%
  summarise()

isos_proj <- st_transform(iso_sf, 32615)

template <- raster(isos_proj, resolution = 100)

iso_surface <- fasterize(isos_proj, template, field = "time", fun = "min")

pal <- colorNumeric("viridis", isos_proj$time, na.color = "transparent")
leaflet() %>%
  addTiles() %>%
  addRasterImage(iso_surface, colors = pal, opacity = 0.5) %>%
  addLegend(values = isos_proj$time, pal = pal,
            title = "Minutes of Travel") %>% 
  addMarkers(lat = c(61.185765, 61.191227, 61.160192), lng = c(-149.883234, -149.928200, -149.939484))