Rewrite amenity-basket around a basket table, per-block travel-time columns, and an auto-picked site.

Author: Nathaniel Henry

docs/tutorials/amenity_basket.md

@@ -17,12 +17,13 @@ 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,500 tokens.*
+*Running this tutorial uses about 3,600 tokens.*
 
 ```{code-cell} python
 :tags: [remove-cell]
 import os
 import pandas as pd
+import numpy as np
 from closecity import Client, close_map
 import plotly.io as pio
 
@@ -33,8 +34,8 @@ close = Client(os.environ.get("CLOSECITY_KEY"))
 
 ## Set up
 
-Read the six category ids from the free catalog, and turn the city name into a centre
-point and a boundary.
+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.
 
 ```python
 from closecity import Client, close_map
@@ -45,14 +46,11 @@ close = Client("ck_live_your_key")   # use your own key here
 ```{code-cell} python
 amenity_types = close.destination_types()
 ids = dict(zip(amenity_types["label"], amenity_types["dest_type_id"]))
-basket = {
-    "supermarket": ids["grocery_stores"],
-    "library": ids["libraries"],
-    "park": ids["parks"],
-    "transit": ids["frequent_transit"],
-    "restaurant": ids["restaurants"],
-    "cafe": ids["cafes"],
-}
+basket = pd.DataFrame({
+    "amenity": ["supermarket", "library", "park", "transit", "restaurant", "cafe"],
+    "dest_type_id": [ids[k] for k in ["grocery_stores", "libraries", "parks",
+                                      "frequent_transit", "restaurants", "cafes"]],
+})
 
 city = close.places(q = "Richmond").iloc[0]
 city_boundary = close.place_boundary(geoid = city["geoid"])
@@ -67,43 +65,50 @@ same way (at a higher token cost).
 
 ```{code-cell} python
 blocks = close.blocks_query(
-    center = {"lon": city["lon"], "lat": city["lat"]},
-    radius_m = 2500,
-    mode = "walk",
-    type = list(basket.values()),
-    include_population = True
-)
-
-one_per_block = blocks.drop_duplicates("geoid").reset_index(drop = True)
+    center = {"lon": city["lon"], "lat": city["lat"]}, radius_m = 2500,
+    mode = "walk", type = basket["dest_type_id"].tolist(), 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 (`NaN` when it is more than 30 minutes
+away). One clean table, no leftover census-block metadata.
+
+```{code-cell} python
+one_per_block = blocks.drop_duplicates("geoid")[["geoid", "geometry"]].reset_index(drop = True)
+population = blocks.drop_duplicates("geoid").set_index("geoid")["population"]
+one_per_block["population"] = one_per_block["geoid"].map(population)
+for _, row in basket.iterrows():
+    times = blocks[blocks["dest_type_id"] == row["dest_type_id"]].set_index("geoid")["travel_time"]
+    one_per_block[f"{row['amenity']}_min"] = one_per_block["geoid"].map(times)
 total_pop = one_per_block["population"].sum()
+
+# The six walk-time columns as a boolean "within 15 minutes" table, for the tallies.
+min_cols = [f"{a}_min" for a in basket["amenity"]]
+within_15 = (one_per_block[min_cols] <= 15).fillna(False)
 ```
 
 ## Coverage, one amenity at a time
 
 For each amenity, a block counts as covered when it is within a 15-minute walk.
 
 ```{code-cell} python
-for name, type_id in basket.items():
-    covered = set(blocks.loc[(blocks.dest_type_id == type_id) &
-                             (blocks.travel_time <= 15), "geoid"])
-    pop = one_per_block.loc[one_per_block.geoid.isin(covered), "population"].sum()
-    print(f"{name:11} {100 * pop / total_pop:3.0f}%")
+for a in basket["amenity"]:
+    covered = within_15[f"{a}_min"]
+    pop = one_per_block.loc[covered, "population"].sum()
+    print(f"{a:11} {100 * pop / total_pop:3.0f}%")
 ```
 
 Parks and restaurants tend to be everywhere; supermarkets and libraries are usually
-the hardest to reach. Map one amenity — every block shown, the covered ones
-highlighted, the city boundary behind.
+the hardest to reach. Map the library coverage — every block shown, the covered ones
+highlighted, the library locations as points, the city boundary behind.
 
 ```{code-cell} python
-near_library = set(blocks.loc[(blocks.dest_type_id == basket["library"]) &
-                              (blocks.travel_time <= 15), "geoid"])
-one_per_block["has_library"] = one_per_block.geoid.isin(near_library)
-close_map(
-    one_per_block,
-    highlight = "has_library",
-    color = "#058040",
-    boundary = city_boundary
-)
+one_per_block["has_library"] = within_15["library_min"]
+library_type = int(basket.loc[basket["amenity"] == "library", "dest_type_id"].iloc[0])
+libraries = close.pois_search(lat = city["lat"], lon = city["lon"], radius_m = 2500,
+                              type = library_type)
+close_map(one_per_block, highlight = "has_library", color = "#058040",
+          points = libraries, boundary = city_boundary)
 ```
 
 ## The 15-minute-city score
@@ -113,34 +118,21 @@ That score, from 0 to 6, is the map planners reach for; blue marks the best-serv
 blocks. It reuses the data you already pulled, so it costs nothing more.
 
 ```{code-cell} python
-covered = blocks[blocks.travel_time <= 15]
-score = covered.groupby("geoid")["dest_type_id"].nunique()
-one_per_block["score"] = one_per_block.geoid.map(score).fillna(0).astype(int)
+one_per_block["score"] = within_15.sum(axis = 1)
 close_map(one_per_block, fill = "score", reverse = True, boundary = city_boundary)
 ```
 
 ## Who can reach all six
 
-A block is fully covered only if all six amenities are within 15 minutes.
+A block is fully covered only when all six amenities are within 15 minutes.
 
 ```{code-cell} python
-covered_all = set(one_per_block.geoid)
-for type_id in basket.values():
-    covered = set(blocks.loc[(blocks.dest_type_id == type_id) &
-                             (blocks.travel_time <= 15), "geoid"])
-    covered_all &= covered
-
-basket_pop = one_per_block.loc[one_per_block.geoid.isin(covered_all),
-                               "population"].sum()
+one_per_block["full_basket"] = one_per_block["score"] == 6
+basket_pop = one_per_block.loc[one_per_block["full_basket"], "population"].sum()
 print(f"All six amenities: {100 * basket_pop / total_pop:.0f}% of residents")
 
-one_per_block["full_basket"] = one_per_block.geoid.isin(covered_all)
-close_map(
-    one_per_block,
-    highlight = "full_basket",
-    color = "#f36e21",
-    boundary = city_boundary
-)
+close_map(one_per_block, highlight = "full_basket", color = "#f36e21",
+          boundary = city_boundary)
 ```
 
 ## Which amenity to add first
@@ -149,43 +141,29 @@ Count how many not-yet-covered residents are missing each amenity. The amenity t
 most people lack is the one to add first.
 
 ```{code-cell} python
-uncovered = set(one_per_block.geoid) - covered_all
-for name, type_id in basket.items():
-    covered = set(blocks.loc[(blocks.dest_type_id == type_id) &
-                             (blocks.travel_time <= 15), "geoid"])
-    lacking = uncovered - covered
-    pop = one_per_block.loc[one_per_block.geoid.isin(lacking), "population"].sum()
-    print(f"{name:11} {pop:6.0f} residents would gain access")
+for a in basket["amenity"]:
+    lacking = ~one_per_block["full_basket"] & ~within_15[f"{a}_min"]
+    print(f"{a:11} {one_per_block.loc[lacking, 'population'].sum():6.0f} residents would gain access")
 ```
 
 ## 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. Build a block-by-amenity
-coverage grid, 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.
+The residents worth targeting first are the ones *almost* there — a block with five
+of the six is a much easier win than one with 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.
 
 ```{code-cell} python
-import numpy as np
-
-names = list(basket)
-covered_by = {name: set(blocks.loc[(blocks.dest_type_id == basket[name]) &
-                                   (blocks.travel_time <= 15), "geoid"])
-              for name in names}
-geoids = one_per_block["geoid"].to_numpy()
+have = within_15.to_numpy()
+amenities = list(basket["amenity"])
+n_missing = 6 - one_per_block["score"].to_numpy()
+gap = np.array([" + ".join(a for a, h in zip(amenities, row) if not h) for row in have])
 pops = one_per_block["population"].to_numpy()
 
-# One row per block: which amenities are missing (not within 15 minutes)?
-missing = np.array([[g not in covered_by[name] for name in names] for g in geoids])
-n_missing = missing.sum(axis = 1)
-gap = np.array([" + ".join(n for n, m in zip(names, row) if m) for row in missing])
-
 almost = np.isin(n_missing, [1, 2])
 almost_pop = pops[almost].sum()
-print(
-    f"{100 * almost_pop / total_pop:.0f}% of residents are one or two amenities "
-    f"short of the full basket.\n"
-)
+print(f"{100 * almost_pop / total_pop:.0f}% of residents are one or two amenities "
+      f"short of the full basket.\n")
 
 by_gap = (pd.Series(pops[almost], index = gap[almost])
           .groupby(level = 0).sum().sort_values(ascending = False))
@@ -195,38 +173,31 @@ for g, p in by_gap.items():
 
 ## Site a new supermarket
 
-The counts above say *which* amenity to add; the next question is *where*. Take a
-candidate site on land north of the James River and ask how many residents would newly
-gain a supermarket within a 15-minute walk if one opened there. A `direction="to"`
-isochrone gives exactly the blocks that could reach the site on foot in 15 minutes.
+The counts above say *which* amenity to add; the next question is *where*. Pick the
+candidate site automatically: a populated block with 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.
 
 ```{code-cell} python
-site_lon, site_lat = -77.437, 37.548
-reachable = close.isochrone(
-    lon = site_lon,
-    lat = site_lat,
-    mode = "walk",
-    direction = "to",
-    minutes = 15,
-    format = "blocks"
-)
-
-near_supermarket = set(blocks.loc[(blocks.dest_type_id == basket["supermarket"]) &
-                                  (blocks.travel_time <= 15), "geoid"])
-newly_served = set(reachable.geoid) - near_supermarket
-gain_pop = one_per_block.loc[one_per_block.geoid.isin(newly_served),
-                             "population"].sum()
+uncovered = one_per_block[~within_15["supermarket_min"] & (one_per_block["population"] > 0)]
+points = uncovered.geometry.representative_point()
+middle = one_per_block.geometry.union_all().centroid
+site = points.iloc[(((points.x - middle.x) ** 2 + (points.y - middle.y) ** 2)).to_numpy().argmin()]
+
+reachable = close.isochrone(lon = site.x, lat = site.y, mode = "walk",
+                            direction = "to", minutes = 15, format = "blocks")
+near_supermarket = set(one_per_block.loc[within_15["supermarket_min"], "geoid"])
+newly_served = (set(reachable["geoid"]) & set(one_per_block["geoid"])) - near_supermarket
+gain_pop = one_per_block.loc[one_per_block["geoid"].isin(newly_served), "population"].sum()
 print(f"A supermarket here would newly serve {gain_pop:.0f} residents")
 ```
 
-Map the whole city and highlight the blocks that would newly gain access.
+Map the whole city and highlight the blocks that would newly gain access, with the
+candidate site marked by an X.
 
 ```{code-cell} python
-one_per_block["newly_served"] = one_per_block.geoid.isin(newly_served)
-close_map(
-    one_per_block,
-    highlight = "newly_served",
-    color = "#e8590c",
-    boundary = city_boundary
-)
+one_per_block["newly_served"] = one_per_block["geoid"].isin(newly_served)
+close_map(one_per_block, highlight = "newly_served", color = "#e8590c",
+          mark = (site.x, site.y), boundary = city_boundary)
 ```