#I choose Cambridge, MA, as my location, and museum and peruvain_restaurant as my place types
#Cambridge boundary
cambridge <- tigris::places(state = "MA") %>%
  filter(NAME == "Cambridge")
## Retrieving data for the year 2024
#Block Groups for Middlesex County (ACS 2024 5-year)
bg_middlesex <- suppressMessages(
  tidycensus::get_acs(
    geography = "block group",
    state = "MA",
    county = "Middlesex",
    variables = c(hhincome = "B19013_001"),
    year = 2024,
    survey = "acs5",
    geometry = TRUE,
    output = "wide"
  )
)

#Keep only BGs that intersect Cambridge
cambridge_trans <- cambridge %>% st_transform(st_crs(bg_middlesex))
bg_cambridge <- bg_middlesex[cambridge_trans, ]
# Coordinate reference systems and radius calculation function
gcs_id <- 4326
pcs_id <- 26919

getXYRadius <- function(polygon, gcs_id, pcs_id){
  if (st_crs(polygon) != st_crs(pcs_id)){
    polygon <- polygon %>% st_transform(pcs_id)
  }
  
  bb <- st_bbox(polygon)
  bb_corner <- st_point(c(bb[1], bb[2])) %>% st_sfc(crs = pcs_id)
  bb_center <- bb %>% st_as_sfc() %>% st_centroid()
  r <- st_distance(bb_center, bb_corner)
  
  bb_center_gcs <- bb_center %>% st_transform(gcs_id)
  xy <- bb_center_gcs %>% st_coordinates() %>% as.vector()
  
  data.frame(x = xy[1], y = xy[2], r = as.numeric(r))
}

# Pre-allocate and calculate for all Cambridge BGs
bg_cambridge_xyr <- data.frame(
  x = numeric(nrow(bg_cambridge)),
  y = numeric(nrow(bg_cambridge)),
  r = numeric(nrow(bg_cambridge))
)

for (i in 1:nrow(bg_cambridge)){
  bg_cambridge_xyr[i, ] <- getXYRadius(bg_cambridge[i, ], gcs_id, pcs_id)
}
#Create a point sf object in PCS, and buffer the points
search_points_pcs <- bg_cambridge_xyr %>%
  st_as_sf(coords = c("x", "y"), crs = gcs_id) %>%
  st_transform(pcs_id)

search_buffers_pcs <- st_buffer(search_points_pcs, dist = bg_cambridge_xyr$r)

#visualize it
tm_shape(cambridge) +
  tm_borders(col = "black", lwd = 2) +
  tm_shape(search_buffers_pcs) +
  tm_polygons(col = "blue", alpha = 0.2, border.col = "blue")
## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
## This message is displayed once every 8 hours.
#ensure all 8 fields are included, with help from R assistant (i could not get the search circle to work)
nearbySearch <- function(lat, lon, radius, types_vec, fieldmask_vec, google_api_key){
  endpoint <- "https://places.googleapis.com/v1/places:searchNearby"
  
  body <- list(
    includedTypes = as.list(types_vec),
    locationRestriction = list(
      circle = list(
        center = list(latitude = lat, longitude = lon),
        radius = radius
      )
    ),
    rankPreference = "DISTANCE"
  )
  
  resp <- POST(
    endpoint,
    add_headers(
      "Content-Type" = "application/json",
      "X-Goog-Api-Key" = google_api_key,
      "X-Goog-FieldMask" = paste(fieldmask_vec, collapse = ",")
    ),
    body = body,
    encode = "json"
  )
  
  parsed <- content(resp, as = "text") %>% jsonlite::fromJSON(flatten = TRUE)
  
  if (is.null(parsed$places) || length(parsed$places) == 0) {
    return(data.frame())
  }
  
  as.data.frame(parsed$places)
}

# Required field mask elements
fields_required <- c(
  "places.id",
  "places.displayName",
  "places.formattedAddress",
  "places.location",
  "places.types",
  "places.priceLevel",
  "places.rating",
  "places.userRatingCount"
)

# Call API across all search circles
cambridge_poi_list <- vector("list", nrow(bg_cambridge_xyr))

for (i in seq_len(nrow(bg_cambridge_xyr))) {
  cambridge_poi_list[[i]] <- nearbySearch(
    lat = bg_cambridge_xyr$y[i],
    lon = bg_cambridge_xyr$x[i],
    radius = bg_cambridge_xyr$r[i],
    types_vec = c("peruvian_restaurant", "museum"),
    fieldmask_vec = fields_required,
    google_api_key = Sys.getenv("GOOGLE_API_KEY")
  )
  Sys.sleep(0.5) # Courtesy pause between requests
}

# Combine and remove duplicate places returned across overlapping buffers
cambridge_pois_raw <- dplyr::bind_rows(cambridge_poi_list)

# De-duplicate by unique place ID
cambridge_pois <- cambridge_pois_raw %>%
  distinct(id, .keep_all = TRUE)

# Export to RDS
readr::write_rds(cambridge_pois, "cambridge_pois.rds")
#latitude and longitude
cambridge_pois_sf <- cambridge_pois %>%
  filter(!is.na(location.longitude) & !is.na(location.latitude)) %>%
  st_as_sf(coords = c("location.longitude", "location.latitude"), crs = 4326)

#Clip to city boundary
cambridge_pois_sf <- cambridge_pois_sf[st_transform(cambridge, 4326), ]

tm_shape(cambridge) +
  tm_borders(lwd = 1.5) +
  tm_shape(cambridge_pois_sf) +
  tm_dots(fill = "rating",
          fill.scale = tm_scale_continuous(values = "viridis"),
          size = "userRatingCount",
          size.scale = tm_scale_continuous(values.scale = 2),
          col = "black",
          shape = 21,
          popup.vars = c("Name" = "displayName.text",
                         "Address" = "formattedAddress",
                         "Rating" = "rating",
                         "Reviews" = "userRatingCount"))
#I choose Cambridge, MA, as my location, and museum and peruvian_restaurant as my place types
nrow(cambridge_pois_sf)
## [1] 23
#there are 23 rows in my dataset
#museums cluster around Harvard Square/Cambridge Common and MIT/Kendall Square, which is not surprising as this is something you would expect within a university
#Peruvian restaurants tend to be located along major commercial corridors like Massachusetts Avenue, which is generally where restaurants tend to be. I don't think Cambridge is large enough (only around 6.5 sq miles) to notice any diaspora trends. I would need to observe the entire Boston metropolitan area in order to likely get noticeable clusters.