I’ve wanted to merge two ideas for a while. I did so a month or so back and wanted to share via a diary page (which is my sole programming blog outlet at this time).

I want to show a webmap that has overlapping isochrones where it shows the shortest distance between multiple zones. Now, this was sort of possible for me before with vectorized rings but I want it as a raster so I can show the spectrum of time it would take. It makes it more detailed and I think more pretty. Below is an example using the mapbox api for doing so. I assume like most web services it will be deprecated in the future but who knows.

library(mapboxapi)
library(tidyverse)
library(mapdeck)
library(osmdata)
library(leaflet)
library(sf)
library(raster)
library(rgeos)
library(fasterize)
library(rgdal)
library(htmlwidgets)

#mapboxapi::mb_access_token("put_your_own_key_here",
#                           install = TRUE, overwrite = TRUE)


x <- getbb("Anchorage, AK") %>% opq() %>%
  add_osm_feature("shop", "supermarket") %>%
  osmdata_sf()
#y <- getbb("Anchorage, AK") %>% opq() %>%
#  add_osm_feature("shop", "wholesale") %>%
#  osmdata_sf()
#x <- bind_rows(x$osm_polygons, y$osm_polygons)

x <- x$osm_polygons

all_latlon <- c()
for (i in x$osm_id) {
  
  latlon <-
    paste0(x$geometry[[i]][1][[1]][1,1], ", ",
           x$geometry[[i]][1][[1]][1,2])
  
  all_latlon <- c(latlon, all_latlon)
}

data.frame(name = x$name, latlon = all_latlon)



all_isos <- data.frame()

for(j in all_latlon) {
  isos <- mb_isochrone(
    location = j,
    profile = "walking",
    time = 1:20,
  )  
  print(j)
  all_isos <- rbind(all_isos, isos)
  Sys.sleep(1)
}


library(sf)
library(dplyr)

iso_sf <- st_as_sf(all_isos)

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")

(map <- leaflet() %>%
  addTiles() %>%
  addRasterImage(iso_surface, colors = pal, opacity = 0.5) %>%
  addLegend(values = isos_proj$time, pal = pal,
            title = "Minutes of Travel") %>%
  addPolygons(data = x, popup = ~name,
              opacity = 1.0, fillOpacity = 0.5, 
              color = "orange"))

writeRaster(iso_surface, filename = "supermarkets.tif")


saveWidget(map, file="grocery_stores.html")