## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = identical(tolower(Sys.getenv("NOT_CRAN")), "true"),
  out.width = "100%"
)

## ----message = FALSE----------------------------------------------------------
options(java.parameters = "-Xmx2G")

library(r5r)
library(accessibility)
library(sf)
library(data.table)
library(ggplot2)
library(interp)
library(h3jsr)
library(dplyr)

## ----message = FALSE----------------------------------------------------------
# system.file returns the directory with example data inside the r5r package
# set data path to directory containing your own data if not running this example
data_path <- system.file("extdata/poa", package = "r5r")

r5r_network <- build_network(data_path)

## ----message = FALSE----------------------------------------------------------
# read all points in the city
points <- fread(file.path(data_path, "poa_hexgrid.csv"))

# routing inputs
mode <- c("WALK", "TRANSIT")
max_walk_time <- 30      # in minutes
travel_time_cutoff <- 20 # in minutes
time_window <- 60        # in minutes
departure_datetime <- as.POSIXct("13-05-2019 14:00:00",
                                 format = "%d-%m-%Y %H:%M:%S")

# calculate accessibility
access1 <- r5r::accessibility(r5r_network,
                        origins = points,
                        destinations = points,
                        mode = mode,
                        opportunities_colnames = c("schools", "healthcare"),
                        decay_function = "step",
                        cutoffs = travel_time_cutoff,
                        departure_datetime = departure_datetime,
                        max_walk_time = max_walk_time,
                        time_window = time_window,
                        progress = FALSE)

head(access1)

## ----message = FALSE----------------------------------------------------------
# calculate travel time matrix
ttm <- r5r::travel_time_matrix(r5r_network,
                        origins = points,
                        destinations = points,
                        mode = mode,
                        departure_datetime = departure_datetime,
                        max_walk_time = max_walk_time,
                        time_window = time_window,
                        progress = FALSE)
head(ttm)

## ----message = FALSE----------------------------------------------------------
# calculate accessibility
access_edu <- accessibility::cumulative_cutoff(travel_matrix = ttm, 
                                        land_use_data = points,
                                        opportunity = 'schools',
                                        travel_cost = 'travel_time_p50',
                                        cutoff = 20)

access_health <- accessibility::cumulative_cutoff(travel_matrix = ttm, 
                                        land_use_data = points,
                                        opportunity = 'healthcare',
                                        travel_cost = 'travel_time_p50',
                                        cutoff = 20)
head(access_edu)
head(access_health)

## ----message = FALSE, out.width='100%'----------------------------------------
# retrieve polygons of H3 spatial grid
grid <- h3jsr::cell_to_polygon(points$id, simple = FALSE)

# merge accessibility estimates
access_sf <- left_join(grid, access1, by = c('h3_address'='id'))

# plot
ggplot() +
  geom_sf(data = access_sf, aes(fill = accessibility), color= NA) +
  scale_fill_viridis_c(direction = -1, option = 'B') +
  labs(fill = "Number of\nfacilities within\n20 minutes") +
  theme_minimal() +
  theme(axis.title = element_blank()) +
  facet_wrap(~opportunity) +
  theme_void()


## ----message = FALSE, out.width='100%'----------------------------------------
# interpolate estimates to get spatially smooth result
access_schools <- access1 %>% 
  filter(opportunity == "schools") %>%
  inner_join(points, by='id') %>%
  with(interp::interp(lon, lat, accessibility)) %>%
  with(cbind(acc=as.vector(z),  # Column-major order
             x=rep(x, times=length(y)),
             y=rep(y, each=length(x)))) %>% as.data.frame() %>% na.omit() %>%
  mutate(opportunity = "schools")

access_health <- access1 %>% 
  filter(opportunity == "healthcare") %>%
  inner_join(points, by='id') %>%
  with(interp::interp(lon, lat, accessibility)) %>%
  with(cbind(acc=as.vector(z),  # Column-major order
             x=rep(x, times=length(y)),
             y=rep(y, each=length(x)))) %>% as.data.frame() %>% na.omit() %>%
  mutate(opportunity = "healthcare")

access.interp <- rbind(access_schools, access_health)

# find results' bounding box to crop the map
bb_x <- c(min(access.interp$x), max(access.interp$x))
bb_y <- c(min(access.interp$y), max(access.interp$y))

# extract OSM network, to plot over map
street_net <- street_network_to_sf(r5r_network)

# plot
ggplot(na.omit(access.interp)) +
  geom_sf(data = street_net$edges, color = "gray55", size=0.01, alpha = 0.7) +
  geom_contour_filled(aes(x=x, y=y, z=acc), alpha=.7) +
  scale_fill_viridis_d(direction = -1, option = 'B') +
  scale_x_continuous(expand=c(0,0)) +
  scale_y_continuous(expand=c(0,0)) +
  coord_sf(xlim = bb_x, ylim = bb_y, datum = NA) + 
  labs(fill = "Number of\nfacilities within\n20 minutes") +
  theme_void() +
  facet_wrap(~opportunity)

## ----message = FALSE----------------------------------------------------------
r5r::stop_r5(r5r_network)
rJava::.jgc(R.gc = TRUE)