---
title: "The amenity basket"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{The amenity basket}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  eval = nzchar(Sys.getenv("CLOSECITY_KEY"))
)
library(closecity)
library(sf)
close <- closecity::close_client(api_key = Sys.getenv("CLOSECITY_KEY"))
```

A city planner wants every resident to be able to walk to a basket of six everyday
amenities: a supermarket, a library, a park, a frequent-transit stop, a restaurant,
and a cafe. This tutorial measures how many residents already have that, and shows
where the gaps are. The idea follows
[this analysis](https://nathenry.com/writing/2023-02-07-seattle-walkability.html),
here applied to Richmond, Virginia.

*Running this tutorial uses about 3,600 tokens.*

## Set up

Turn the city name into a centre point and a boundary, and build the basket as a
small table: each amenity paired with its destination-type id.

```r
library(closecity)
library(sf)
close <- closecity::close_client(api_key = "ck_live_your_key")   # use your own key here
```

```{r}
amenity_types <- close$destination_types()
ids <- setNames(amenity_types$dest_type_id, amenity_types$label)

basket <- data.frame(
  amenity = c("supermarket", "library", "park", "transit", "restaurant", "cafe"),
  dest_type_id = unname(ids[c("grocery_stores", "libraries", "parks",
                              "frequent_transit", "restaurants", "cafes")])
)

city <- close$places(q = "Richmond")[1, ]
city_boundary <- close$place_boundary(geoid = city$geoid)
```

## Pull the blocks, with population

`$blocks_query()` reads every page: the walk time from every block to each of the
six categories, plus each block's population. To keep this tutorial cheap we take
the central blocks within a radius; `$place_blocks(geoid = city$geoid)` pulls
**every** block in the city the same way (at a higher token cost).

```{r}
blocks <- close$blocks_query(
  center = list(lon = city$lon, lat = city$lat),
  radius_m = 2500,
  mode = "walk",
  type = basket$dest_type_id,
  include_population = TRUE
)
```

Reshape to one row per block: keep just the GEOID and geometry, add the population,
and add a walk-time column for each amenity (`NA` when it is more than 30 minutes
away). The result is one clean table, with no leftover census-block metadata.

```{r}
one_per_block <- blocks[!duplicated(blocks$geoid), "geoid"]
one_per_block$population <-
  setNames(blocks$population, blocks$geoid)[one_per_block$geoid]
for (i in seq_len(nrow(basket))) {
  sub <- blocks[blocks$dest_type_id == basket$dest_type_id[i], ]
  one_per_block[[paste0(basket$amenity[i], "_min")]] <-
    setNames(sub$travel_time, sub$geoid)[one_per_block$geoid]
}
total_pop <- sum(one_per_block$population)

# The six walk-time columns, as a plain matrix, for the tallies below.
mins <- as.matrix(sf::st_drop_geometry(one_per_block)[, paste0(basket$amenity, "_min")])
within_15 <- mins <= 15
within_15[is.na(within_15)] <- FALSE
```

## Coverage, one amenity at a time

For each amenity, a block counts as covered when it is within a 15-minute walk. Add
up the population of the covered blocks.

```{r}
for (i in seq_len(nrow(basket))) {
  pop <- sum(one_per_block$population[within_15[, i]])
  cat(sprintf("%-11s %3.0f%%\n", basket$amenity[i], 100 * pop / total_pop))
}
```

Parks and restaurants tend to be everywhere; supermarkets and libraries are usually
the hardest to reach. Map the library coverage: every block shown, the covered ones
highlighted, the library locations as points, and the city boundary behind.

```{r}
one_per_block$has_library <- within_15[, basket$amenity == "library"]
libraries <- close$pois_search(
  lat = city$lat,
  lon = city$lon,
  radius_m = 2500,
  type = basket$dest_type_id[basket$amenity == "library"]
)
closecity::close_map(
  x = one_per_block,
  highlight = "has_library",
  color = "#058040",
  points = libraries,
  boundary = city_boundary
)
```

## The 15-minute-city score

Count, for each block, how many of the six amenities are within a 15-minute walk.
That score, from 0 to 6, is the map planners reach for; blue marks the best-served
blocks. It reuses the data you already pulled, so it costs nothing more.

```{r}
one_per_block$score <- rowSums(within_15)
closecity::close_map(
  x = one_per_block,
  fill = "score",
  reverse = TRUE,
  boundary = city_boundary
)
```

## Who can reach all six

A block is fully covered only when all six amenities are within 15 minutes.

```{r}
one_per_block$full_basket <- one_per_block$score == 6
basket_pop <- sum(one_per_block$population[one_per_block$full_basket])
cat(sprintf("All six amenities: %.0f%% of residents\n", 100 * basket_pop / total_pop))

closecity::close_map(
  x = one_per_block,
  highlight = "full_basket",
  color = "#f36e21",
  boundary = city_boundary
)
```

## Which amenity to add first

Look at the residents who are not yet fully covered, and count how many of them are
missing each amenity. The amenity that the most people lack is the one to add first.

```{r}
for (i in seq_len(nrow(basket))) {
  lacking <- !one_per_block$full_basket & !within_15[, i]
  cat(sprintf("%-11s %6.0f residents would gain access\n", basket$amenity[i],
              sum(one_per_block$population[lacking])))
}
```

## Who is one or two amenities away

The residents worth targeting first are the ones *almost* there: a block that has
five of the six is a much easier win than one that has none. Count how many amenities
each block is missing, and for the blocks short by just one or two, break down which
amenities are the gap.

```{r}
n_missing <- 6 - one_per_block$score
gap <- apply(within_15, 1, function(row) paste(basket$amenity[!row], collapse = " + "))

almost <- n_missing %in% c(1, 2)
almost_pop <- sum(one_per_block$population[almost])
cat(sprintf("%.0f%% of residents are one or two amenities short of the full basket.\n\n",
            100 * almost_pop / total_pop))

pop_by_gap <- tapply(one_per_block$population[almost], gap[almost], sum)
pop_by_gap <- sort(pop_by_gap, decreasing = TRUE)
for (g in names(pop_by_gap)) {
  cat(sprintf("  missing %-22s %3.0f%%\n", g, 100 * pop_by_gap[g] / almost_pop))
}
```

## Site a new supermarket

The counts above say *which* amenity to add; the next question is *where*. Pick the
candidate site automatically: the uncovered block (no supermarket within a 15-minute
walk) nearest the middle of the study area. A `direction = "to"` isochrone then gives
the blocks that could reach the site on foot in 15 minutes, and the ones not already
served are the population this site would newly reach.

```{r}
uncovered <- one_per_block[!within_15[, basket$amenity == "supermarket"] &
                             one_per_block$population > 0, ]
uncovered_points <- sf::st_point_on_surface(sf::st_geometry(uncovered))
middle <- sf::st_centroid(sf::st_union(sf::st_geometry(one_per_block)))
site <- sf::st_coordinates(
  uncovered_points[which.min(sf::st_distance(uncovered_points, middle))]
)

reachable <- close$isochrone(lon = site[1], lat = site[2], mode = "walk",
                             direction = "to", minutes = 15, format = "blocks")
near_supermarket <- one_per_block$geoid[within_15[, basket$amenity == "supermarket"]]
newly_served <- setdiff(intersect(reachable$geoid, one_per_block$geoid), near_supermarket)
cat(sprintf("A supermarket here would newly serve %.0f residents\n",
            sum(one_per_block$population[one_per_block$geoid %in% newly_served])))
```

Map the whole city and highlight the blocks that would newly gain access, with the
candidate site marked by an X.

```{r}
one_per_block$newly_served <- one_per_block$geoid %in% newly_served
closecity::close_map(
  x = one_per_block,
  highlight = "newly_served",
  color = "#e8590c",
  mark = c(site[1], site[2]),
  boundary = city_boundary
)
```