MA_ZHAOXIN_MINI1_Attempt2

tidycensus::census_api_key(Sys.getenv("CENSUS_API_KEY"))

## To install your API key for use in future sessions, run this function with `install = TRUE`.

install = TRUE

Getting Census BG Boundary

# Get block groups
bg <- suppressMessages(
  tidycensus::get_acs(geography = "block group", 
                      state = "NY",
                      county = c("Broome"),
                      variables = c(hhincome = 'B19013_001'),
                      year = 2023,
                      survey = "acs5", # we typically use American Community Survey 5-year estimates
                      geometry = TRUE, # we need an sf object
                      output = "wide") # wide vs. long
)

# City of Binghamton boundary
binghamton <- tigris::places('NY') %>% filter(NAME == 'Binghamton')

## Retrieving data for the year 2022

# Get BGs intersecting with the City of Binghamton boundary
bg_binghamton <- bg[binghamton,]

# View the data
bg_binghamton %>% head() %>% knitr::kable() # Ignore kable(). This function is for neatly displaying tables on HTML document.

## Warning in attr(x, "align"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

## Warning in attr(x, "format"): 'xfun::attr()' is deprecated.
## Use 'xfun::attr2()' instead.
## See help("Deprecated")

	GEOID	NAME	hhincomeE	hhincomeM	geometry
1	360070128003	Block Group 3; Census Tract 128; Broome County; New York	61771	10296	MULTIPOLYGON (((-75.90157 4…
2	360070016003	Block Group 3; Census Tract 16; Broome County; New York	120577	25839	MULTIPOLYGON (((-75.94818 4…
3	360070014021	Block Group 1; Census Tract 14.02; Broome County; New York	54716	21607	MULTIPOLYGON (((-75.93081 4…
6	360070017005	Block Group 5; Census Tract 17; Broome County; New York	33904	10568	MULTIPOLYGON (((-75.9159 42…
7	360070012001	Block Group 1; Census Tract 12; Broome County; New York	NA	NA	MULTIPOLYGON (((-75.91007 4…
10	360070003003	Block Group 3; Census Tract 3; Broome County; New York	45313	39113	MULTIPOLYGON (((-75.92243 4…

bg_binghamton <- bg_binghamton %>%
  select(GEOID,
         hhincome = hhincomeE) # New name = old name

tmap_mode("view")

## ℹ tmap mode set to "view".

tm_shape(bg_binghamton) + tm_borders(lwd = 2) +
  tm_shape(binghamton) + tm_polygons(col = 'red', alpha = 0.4)

## 
## ── 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`.

Preparing location parameters for the Nearby Search

Using Census BGs

# Function: Get XY coordinates and radius
getXYRadius <- function(polygon, gcs_id, pcs_id){

  # Transform the CRS to PCS. WHY?
  if (st_crs(polygon) != st_crs(pcs_id)){
    polygon <- polygon %>% st_transform(pcs_id)
  }
  
  # Get bounding box of a given polygon
  bb <- st_bbox(polygon)
  
  # Get XY coordinates of any one corner of the bounding box.
  bb_corner <- st_point(c(bb[1], bb[2])) %>% st_sfc(crs = pcs_id)
  
  # Get centroid of the bb
  bb_center <- bb %>% st_as_sfc() %>% st_centroid()

  # Get the distance between bb_center and bb_corner
  r <- st_distance(bb_corner, bb_center)
  
  # Convert the CRS of centroid to GCS. 
  bb_center <- bb_center %>% st_transform(gcs_id)
  
  # Get longitude and latitude
  xy <- bb_center %>% st_coordinates() %>% as.vector()
  
  lon_lat_radius <- data.frame(x = xy[1],
                               y = xy[2],
                               r = r)

  return(lon_lat_radius)
}

We can apply this function to each BG.

# Define EPSG codes for GCS (WGS 84) and PCS 
#UTM 18N for upstate NY
gcs_id <- 4326
pcs_id <- 32618 

# Pre-allocate a data frame. Results will fill this data frame
bg_binghamton_xyr <- data.frame(x = numeric(nrow(bg_binghamton)),
                             y = NA,
                             r = NA)

# Do a for-loop
for (i in 1:nrow(bg_binghamton)){
  bg_binghamton_xyr[i,] <- bg_binghamton[i, ] %>%
    getXYRadius(gcs_id = gcs_id,
                pcs_id = pcs_id)
}

Let’s visualize what we’ve just done.

tmap_mode('view')

## ℹ tmap mode set to "view".

bg_binghamton_xyr %>%
  # Convert the data frame into an sf object
  st_as_sf(coords = c("x", "y"), crs = st_crs(bg_binghamton)) %>% 
  # Draw a buffer centered at the centroid of BG polygons.
  # The buffer distance is the radius we just calculated
  st_buffer(dist = .$r) %>%
  # Display this buffer in red
  tm_shape(.) + tm_polygons(alpha = 0.3, col = 'red') +
  # Display the original polygon in blue
  tm_shape(bg_binghamton) + tm_borders(col= 'blue')

## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.

Another way to divide an area: creating a fishnet grid

# Define EPSG codes for GCS (WGS 84) and PCS (WGS 84 / UTM zone 16N)
gcs_id <- 4326
pcs_id <- 32618

# Get bbox for Binghamton
binghamton_bb <- binghamton %>% st_transform(pcs_id) %>% st_bbox()

# Find coordinates for the four sides of the bbox
west <- binghamton_bb[1]
east <- binghamton_bb[3]
south <- binghamton_bb[2]
north <- binghamton_bb[4]

# Split the bbox into a grid
fishnet_n <- 4
dist <- abs(east - west)/fishnet_n

# Fishnet points
fish_x <- seq(from = west, to = east, by = dist)
fish_y <- seq(from = south, to = north, by = dist)

fishnet <- expand.grid(fish_x, fish_y) %>%
  rename(x = Var1, y = Var2) %>%
  st_as_sf(coords = c('x', 'y'), crs = pcs_id)

fishnet$x <- fishnet %>% st_transform(gcs_id) %>% st_coordinates() %>% .[,1]
fishnet$y <- fishnet %>% st_transform(gcs_id) %>% st_coordinates() %>% .[,2]
fishnet$r <- dist*(2^(1/2)/2) 

# Visualize it
fishnet %>% 
  st_buffer(dist = .$r) %>% 
  tm_shape() +
  tm_polygons(alpha = 0.1, col = 'red') +
  tm_shape(binghamton) + tm_borders(lwd = 3)

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.

Requesting the Google Places API

API Request

library(httr)

google_api_key <- Sys.getenv("GOOGLE_API_KEY")
lat <- 42.0987   # Binghamton city center
lon <- -75.9180 # Bing city center
radius <- 1609.34 # 1 mile in meters

endpoint <- "https://places.googleapis.com/v1/places:searchNearby" # Nearby Search endpoint

body <- list(
  includedTypes = list("mexican_restaurant"),
  locationRestriction = list(
    circle = list(
      center = list(latitude = lat, longitude = lon),
      radius = radius
    )
  )
)

resp <- POST(
  endpoint,
  add_headers(
    "Content-Type" = "application/json",
    "X-Goog-Api-Key" = google_api_key,
    "X-Goog-FieldMask" = "places.displayName,places.formattedAddress,places.types"),
  body = body,
  encode = "json"
)

API Response

HTTP response object

print(resp) # metadata + body

## Response [https://places.googleapis.com/v1/places:searchNearby]
##   Date: 2025-10-03 00:51
##   Status: 200
##   Content-Type: application/json; charset=UTF-8
##   Size: 368 B
## {
##   "places": [
##     {
##       "types": [
##         "mexican_restaurant",
##         "restaurant",
##         "food",
##         "point_of_interest",
##         "establishment"
##       ],
## ...

Content from the response

# Extract content from a request
data <- content(resp, as="text")

# Parse JSON into a list and turn it into a data frame
data <- jsonlite::fromJSON(data, flatten = T) %>% as.data.frame()

names(data)

## [1] "places.types"                    "places.formattedAddress"        
## [3] "places.displayName.text"         "places.displayName.languageCode"

print(data)

##                                                             places.types
## 1 mexican_restaurant, restaurant, food, point_of_interest, establishment
##                              places.formattedAddress
## 1 1171 Vestal Ave suite b, Binghamton, NY 13903, USA
##       places.displayName.text places.displayName.languageCode
## 1 Hacienda Mexican Restaurant                              en

Custom API function

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"
  )
  
  data <- content(resp, as="text") %>% 
    jsonlite::fromJSON(flatten = T) %>% 
    as.data.frame()
  
  if (nrow(data) == 20){
    print("WARNING: The response has 20 rows! Consider using a smaller spatial unit.")
  }
  
  return(data)
}

nearbySearch(lat = 42.0987,
             lon = -75.9180,
             radius = 1609.34,
             types_vec = c("mexican_restaurant"),
             fieldmask_vec = c("places.displayName",
                               "places.formattedAddress",
                               "places.types"),
             google_api_key = Sys.getenv("GOOGLE_API_KEY"))

Downloading and saving the POI data

Applying the function to all Census BGs

# pre-allocate list
data_list <- vector("list", nrow(bg_binghamton_xyr))

for (i in seq_len(nrow(bg_binghamton_xyr))) {
  
  data_list[[i]] <- nearbySearch(
    lat = bg_binghamton_xyr$y[i],
    lon = bg_binghamton_xyr$x[i],
    radius = bg_binghamton_xyr$r[i],
    types_vec = c("mexican_restaurant"), # hacienda mexican place
    fieldmask_vec = c("places.id",
                      "places.displayName",
                      "places.formattedAddress",
                      "places.location",
                      "places.types",
                      "places.primaryType",
                      "places.businessStatus",
                      "places.priceLevel",
                      "places.priceRange",
                      "places.rating",
                      "places.userRatingCount",
                      "places.reviews",
                      "places.reviewSummary",
                      "places.delivery",
                      "places.dineIn",
                      "places.takeout",
                      "places.menuForChildren",
                      "places.outdoorSeating",
                      "places.allowsDogs"),
    google_api_key = Sys.getenv("GOOGLE_API_KEY")
  )
  
  Sys.sleep(1) 
}

# Combine all data frames
data_all <- dplyr::bind_rows(data_list)

Saving the data

saveRDS(data_all, here('google_poi_data.rds'))

Glimpse

Convert to sf using the columns produced by js$places

data_all_sf <- data_all %>%
  rename(x = places.location.longitude, y = places.location.latitude) %>% 
  filter(!is.na(x) & !is.na(y)) %>%
  st_as_sf(coords = c("x", "y"), crs = 4326)

# Map
tm_shape(data_all_sf) + 
  tm_dots(col = "places.rating", 
          size = "places.userRatingCount",
          palette = "magma",
          popup.vars = c("Name" = "places.displayName.text",
                         "Rating" = "places.rating",
                         "Rating Count" = "places.userRatingCount")) +
  tm_shape(binghamton) + 
  tm_borders()

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_dots()`: migrate the argument(s) related to the scale of the
## visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).

#final questions for mini assignment answer
#What city did you choose? ANSWER: Binghamton, NY :)
#Which two place types? ANSWER: I chose parks and museums.
cat("How many rows?             ", nrow(data_all), "\n", sep = "")

## How many rows?             21

#Briefly describe the clusters/gaps I see (spatial patterns)
#ANSWER: I see a blue boundary that represents Binghamton, NY city with orange circles which represent the buffers from my block group (BG) centroids and finally I see black dots that represent centroids!!! This was a cool lab to do. Thanks!