# Get block groups for Jefferson County (Birmingham lies in this county)
bg <- suppressMessages(
tidycensus::get_acs(
geography = "block group",
state = "AL",
county = "Jefferson",
variables = c(hhincome = 'B19013_001'),
year = 2023,
survey = "acs5",
geometry = TRUE,
output = "wide"
)
)
# Get Birmingham boundary
bham <- tigris::places("AL") %>% filter(NAME == "Birmingham")
## Retrieving data for the year 2024
# Filter block groups that intersect Birmingham
bg_bham <- bg[bham,]
# Keep only relevant fields
bg_bham <- bg_bham %>% select(GEOID, hhincome = hhincomeE)
# Map to check boundary vs. BGs
tmap_mode("view")
## ℹ tmap mode set to "view".
tm_shape(bg_bham) + tm_borders(lwd = 2) +
tm_shape(bham) + 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`.
## This message is displayed once every 8 hours.
# Function to compute centroid + radius for Nearby Search queries
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_corner, bb_center)
bb_center <- bb_center %>% st_transform(gcs_id)
xy <- bb_center %>% st_coordinates() %>% as.vector()
lon_lat_radius <- data.frame(x = xy[1], y = xy[2], r = r)
return(lon_lat_radius)
}
# Define coordinate reference systems
gcs_id <- 4326 # geographic (lat/lon)
pcs_id <- 32616 # projected UTM zone for Alabama
# Pre-allocate results
bg_bham_xyr <- data.frame(x = numeric(nrow(bg_bham)),
y = NA,
r = NA)
# Apply function to each BG
for (i in 1:nrow(bg_bham)){
bg_bham_xyr[i,] <- bg_bham[i, ] %>%
getXYRadius(gcs_id = gcs_id, pcs_id = pcs_id)
}
# Visualize buffers to verify full coverage
tmap_mode("view")
## ℹ tmap mode set to "view".
bg_bham_xyr %>%
st_as_sf(coords = c("x", "y"), crs = st_crs(bg_bham)) %>%
st_buffer(dist = .$r) %>%
tm_shape(.) + tm_polygons(alpha = 0.3, col = 'red') +
tm_shape(bg_bham) + tm_borders(col= 'blue')
##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: use `fill_alpha` instead of `alpha`.
# Function to query Google Places Nearby Search
nearbySearch <- function(lat, lon, radius, types_vec, fieldmask_vec, google_api_key){
endpoint <- "https://places.googleapis.com/v1/places:searchNearby"
# Request body
body <- list(
includedTypes = as.list(types_vec),
locationRestriction = list(
circle = list(
center = list(latitude = lat, longitude = lon),
radius = radius
)
),
rankPreference = "DISTANCE"
)
# Send API request
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"
)
# Parse response into data frame
data <- content(resp, as="text") %>%
jsonlite::fromJSON(flatten = T) %>%
as.data.frame()
# Warning if 20-row limit reached
if (nrow(data) == 20){
print("WARNING: Response has 20 rows. Consider smaller radius.")
}
return(data)
}
# Load Google API key from environment
google_api_key <- Sys.getenv("google_API")
# Pre-allocate list for results
data_list <- vector("list", nrow(bg_bham_xyr))
# Query for each BG centroid
for (i in seq_len(nrow(bg_bham_xyr))) {
data_list[[i]] <- nearbySearch(
lat = bg_bham_xyr$y[i],
lon = bg_bham_xyr$x[i],
radius = bg_bham_xyr$r[i],
types_vec = c("school", "hospital"), # chosen POI types
fieldmask_vec = c("places.id",
"places.displayName",
"places.formattedAddress",
"places.location",
"places.types",
"places.priceLevel",
"places.rating",
"places.userRatingCount"),
google_api_key = google_api_key
)
Sys.sleep(1) # pause to avoid API rate limits
}
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
## [1] "WARNING: Response has 20 rows. Consider smaller radius."
# Combine all results
data_all <- dplyr::bind_rows(data_list)
# 1) Flatten the list-column into a single string per row
data_all$types_chr <- sapply(data_all$places.types, function(x) {
if (is.null(x)) return(NA_character_)
if (is.list(x)) return(paste(unlist(x), collapse = ";"))
return(as.character(x))
}, USE.NAMES = FALSE)
# 2) Create a category column based on types
data_all <- data_all %>%
mutate(category = case_when(
str_detect(types_chr, regex("\\bschool\\b", ignore_case = TRUE)) ~ "School",
str_detect(types_chr, regex("\\bhospital\\b", ignore_case = TRUE)) ~ "Hospital",
TRUE ~ "Other"
))
## Warning: There were 2 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `category = case_when(...)`.
## Caused by warning in `stri_detect_regex()`:
## ! argument is not an atomic vector; coercing
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 1 remaining warning.
# 3) Convert to sf (only rows with valid lat/lon)
pois_sf <- data_all %>%
filter(!is.na(places.location.longitude) & !is.na(places.location.latitude)) %>%
st_as_sf(coords = c("places.location.longitude", "places.location.latitude"), crs = 4326, remove = FALSE)
# 4) Counts for legend
school_count <- sum(pois_sf$category == "School", na.rm = TRUE)
hospital_count <- sum(pois_sf$category == "Hospital", na.rm = TRUE)
legend_labels <- c(paste0("Schools (", school_count, ")"),
paste0("Hospitals (", hospital_count, ")"))
# 5) Map POIs with different colors and legend box
tmap_mode("view")
## ℹ tmap mode set to "view".
tm_shape(bham) +
tm_borders(lwd = 2) +
tm_shape(pois_sf) +
tm_dots(col = "category",
palette = c("School" = "blue", "Hospital" = "red"),
size = 0.3,
title = "POI Type") +
tm_add_legend(type = "symbol",
labels = legend_labels,
col = c("blue", "red"),
title = "Identified POIs",
position = c("right", "bottom"))
##
## ── 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>).[tm_dots()] Argument `title` unknown.[v3->v4] `tm_add_legend()`: use `type = "symbols"` instead of `type =
## "symbol"`.[v3->v4] `tm_add_legend()`: use `fill` instead of `col` for the fill color of
## symbols.Multiple palettes called "blue" found: "kovesi.blue", "tableau.blue". The first one, "kovesi.blue", is returned.
## Registered S3 method overwritten by 'jsonify':
## method from
## print.json jsonlite