Grid Search: Finding the Most Convenient Block to Live on in New York City

library(tidyverse)
library(patchwork)
library(sf)
library(osmdata)
library(leaflet)
library(leaflet.extras)
library(htmltools)
library(knitr)
library(kableExtra)
library(stringi)

options(dplyr.summarise.inform = FALSE)
expand.grid.df <- function(...) Reduce(function(...) merge(..., by=NULL), list(...))

coord_sys <- 4326

# shape files
boros <- read_sf("Borough Boundaries")
streets <- read_sf("SimplifiedStreets.shp")
zips <- read_sf("ZIP_CODE_040114") 
zips <- zips %>% st_transform(coord_sys) %>% select(geometry, zip = ZIPCODE) 

# location grading system
points_sys <- read_csv("points_sys.csv")
points_sys <- points_sys %>% mutate(category = str_replace_all(category, "_", " "))

boro_names <- c("New York", "Bronx", "Queens", "Kings")

# take out SI (too suburban), and only take largest land mass of remaining four
boros <- boros %>%
  filter(boro_name != "Staten Island") %>%
  st_cast("POLYGON") %>%
  mutate(area = st_area(.)) %>%
  arrange(desc(area)) %>%
  group_by(boro_name) %>%
  slice(1) %>%
  select(geometry, borough = boro_name)

# we only want NYC streets
streets_nyc <- streets %>% 
  filter(LeftCounty %in% boro_names | RightCount %in% boro_names) %>%
  select(Label, geometry) %>%
  arrange(str_detect(Label, "Ave"))

# share a coordinate system so we can compare/intersect them
boros <- st_transform(boros, st_crs(streets))

# grids are what remains after we remove the (buffed) streets
grids_nyc <- st_difference(
  boros, 
  streets_nyc %>% st_buffer(10) %>% st_union()
) %>%
  st_cast("POLYGON") %>%
  mutate(area = round(st_area(.))) %>%
  # area restrictions should limit us to true blocks
  mutate(grid_idx = row_number()) %>%
  # add the street intersections for identification purposes
  {.$streets <- (.) %>%
    st_intersects(streets_nyc %>% 
                    st_buffer(10 + 1)) %>%
    # paste together the indexed labels
    map(~streets_nyc$Label[.x])
  .
  } %>%
  st_transform(crs = st_crs(coord_sys))



# bounding box for NYC
nyc_bb <- c(-74.2588430, 40.4765780, -73.7002330, 40.9176300)

# places recognized by OSM
all_points <- opq(bbox = nyc_bb, nodes_only = TRUE, timeout = 60*60)

# get every location that falls into any of amenity, shop, leisure
amenity_points <- all_points %>%
  add_osm_feature("amenity", available_tags("amenity")) %>%
  osmdata_sf()

shop_points <- all_points %>%
  add_osm_feature("shop", available_tags("shop")) %>%
  osmdata_sf()

leisure_points <- all_points %>%
  add_osm_feature("leisure", available_tags("leisure")) %>%
  osmdata_sf()

combo_points <- bind_rows(
  leisure_points$osm_points,
  shop_points$osm_points,
  amenity_points$osm_points
) %>%
  select(name, amenity, leisure, shop, geometry) %>%
  filter(!is.na(name)) %>%
  distinct() %>%
  transmute(name, category = coalesce(amenity, leisure, shop)) %>%
  # combine some very similar categories
  mutate(category = case_when(
    category == "pub" ~ "bar",
    category == "chemist" ~ "pharmacy",
    category == "pastry" ~ "bakery",
    TRUE ~ category
  ))

# intersections between grid squares and points
grid_point_ix <- grids_nyc %>% 
  # get coordinates of overlap
  st_intersects(combo_points) %>%
  # cross join each grid to its overlapping points
  imap_dfr(~{
    
    grids_nyc[.y, ] %>% 
      rename(geometry_grid = geometry) %>%
      as.data.frame() %>%
      expand.grid.df(combo_points %>%
                       slice(.x) %>% 
                       rename(geometry_pt = geometry) %>%
                       as.data.frame()
      )
  }
  )

save(grids_nyc, grid_point_ix, file = "grid_items.Rdata")
load("grid_items.Rdata")

# OSM string encoding came in wonky (Latin1?)
fix_string <- function(strings) {
  map_chr(strings, ~
            # no clue why this works
            rawToChar(charToRaw(stri_encode(charToRaw(.x), from = "UTF-8", to = "latin1")))
          )
}


# giant pipe chain that formats blocks for mapping
g <- grid_point_ix %>%
  group_by(grid_idx, category) %>%
  mutate(name = fix_string(name)) %>%
  # total points and list by grid and category
  summarize(n = n(), locs = list(name)) %>%
  mutate(category = str_replace_all(category, "_", " ")) %>%
  inner_join(points_sys %>% select(-n), by = "category") %>%
  mutate(points = pmin(max, n*per)) %>%
  group_by(grid_idx) %>%
  # total points and distinct categories by grid
  mutate(total_points = sum(points), distinct_cat = n_distinct(category)) %>%
  mutate(category = factor(category, levels = points_sys$category)) %>%
  arrange(grid_idx, category) %>%
  group_by(grid_idx, category) %>%
  # gluing together html labels
  mutate(cell1 = paste0(category, " (", n, "): "),
         cell2 = paste0(unlist(locs), collapse = ", "),
         both_cells = paste0("<td><nobr>", cell1, "</nobr></td>", "<td>", cell2, "</td>"),
         full_row = paste0("<tr>", both_cells, "</tr>")) %>%
  group_by(grid_idx, total_points, distinct_cat) %>%
  # collapsing into html table for each grid
  summarize(full_table = paste0('<div class="scrollable">', '<table>', paste0(full_row, collapse = ""), '</table>', '</div>')) %>%
  ungroup() %>%
  left_join(grids_nyc, ., by = c("grid_idx")) %>%
  mutate(total_points = coalesce(total_points, 0),
         distinct_cat = coalesce(distinct_cat, 0),
         full_table = coalesce(full_table, "")) %>%
  # various ranks
  mutate(percentile = ntile(total_points, 100),
         rank = rank(-total_points, ties.method = "min"),
         max_rank = max(rank)) %>%
  group_by(grid_idx) %>%
  # convert area to sq ft for text
  mutate(area_text = formattable::comma(as.numeric(area)*10.764, digits = 0)) %>%
  # final html blob
  mutate(summary = paste0(
    "<b>Block #", grid_idx, " in ", borough, " (", area_text, " sq ft)</b><br>",
    "Enclosed by ", paste0(head(unlist(streets), 6), collapse = ", "), "<br><br>",
    distinct_cat, " location types, ", total_points, " total points (", toOrdinal::toOrdinal(rank), " out of ", max_rank, ")<br><br>")) %>%
  mutate(popup = paste0(summary, full_table)) %>%
  arrange(desc(total_points)) %>%
  # best binning for colors
  group_by(five_plus = total_points >= 5) %>%
  mutate(points_bin = ifelse(!five_plus, 
                             as.numeric(total_points > 0), 
                             ntile(total_points, 13) + 1)) %>%
  ungroup()

colfunc <- colorRampPalette(c("gray",
                              "darkolivegreen4",
                              "lawngreen",
                              "yellow",
                              "orange",
                              "firebrick1",
                              "darkred"))

factpal <- colorFactor(colfunc(15), sort(unique(g$points_bin)))

simple_levels <- factor(c('0', '5', '10', '15', '20+'), levels =  c('0', '5', '10', '15', '20+'))
factpal2 <- colorFactor(colfunc(length(simple_levels)), simple_levels)

block_map <- leaflet() %>% 
  addProviderTiles(providers$Stadia.StamenTonerLite) %>% 
  addCircles(data = grid_point_ix %>% 
               filter(grid_idx == 15883) %>%
               slice(-22) %>% # chase dupe
               st_set_geometry("geometry_pt"),
             label = ~name, radius = 1,
             labelOptions = labelOptions(noHide = TRUE, textOnly = TRUE, offset = c(-6, 0),
                                         direction = "left",
                                         textsize = "15px"))

big_map <- leaflet(height = "700px") %>% 
  addProviderTiles(providers$Stadia.StamenTonerLite) %>% 
  # zip data for search function
  addPolygons(data = zips,
              label = ~zip, group = "hidden",
              stroke = 0, fillOpacity = 0,
              labelOptions = labelOptions(opacity = 0)) %>%
  # primary block data
  addPolygons(data = g[, ],
              popup = ~popup,
              color = ~factpal(points_bin),
              stroke = FALSE,
              fillOpacity = 0.7,
              popupOptions = popupOptions(minWidth = 500)) %>%
  # center at NYC
  setView(lng = -73.95, lat = 40.75, zoom = 13) %>%
  # settings for the zip search functionality
  addSearchFeatures(targetGroups  = 'hidden',
                    options = searchFeaturesOptions(
                      textPlaceholder = 'Search by ZIP',
                      zoom = 15,
                      collapsed = FALSE,
                      position = "topright")
  ) %>%
  addResetMapButton() %>%
  addLegend("bottomright", pal = factpal2, values = simple_levels,
            title = "total points",
            opacity = 1
  )

table <- g %>%
  # only want real blocks no huge lots
  filter(as.numeric(area)*10.764 < 2e5) %>% 
  head(10) %>% 
  st_set_geometry(NULL) %>%
  # can't remember why double quotes don't work
  mutate(full_table = str_replace_all(full_table, '"', "'")) %>%
  select(borough, streets, points = total_points, full_table) %>%
  rowwise() %>%
  mutate(streets = paste0(unlist(streets), collapse = ", ")) %>%
  ungroup() %>%
  # details element allows for accordion reveal
  mutate(points = paste0(points, "<details><summary>locations</summary>", full_table, "</details>")) %>%
  select(-full_table) %>%
  kbl(escape = FALSE) %>%
  kable_styling()

#save(table, big_map, file = "C:/Users/Walker Harrison/Documents/GitHub/website/content/posts/2024-01-13-block-draft/page_elements.Rdata")
save(table, big_map, file = "page_elements.Rdata")