Data

• Census ACS “hhincome” data for 2022 in Fulton and Dekalb Counties
• Hospital POI data from Yelp within the two counties

tract <- suppressMessages(
  get_acs(geography = 'tract',
          state = 'GA',
          county = c('Fulton','Dekalb'),
          variables = c(hhincome = 'B19019_001'),
          year = 2022,
          survey = "acs5", 
          geometry = TRUE,
          output = "wide")) %>% 
  select(GEOID, hhincome = hhincomeE)

## Warning: • You have not set a Census API key. Users without a key are limited to 500
## queries per day and may experience performance limitations.
## ℹ For best results, get a Census API key at
## http://api.census.gov/data/key_signup.html and then supply the key to the
## `census_api_key()` function to use it throughout your tidycensus session.
## This warning is displayed once per session.

tract <- tract %>% 
  mutate(inc_category = case_when(
    hhincome >= 200000 ~ "200,000-300,000",
    hhincome >= 100000 ~ "100,000-200,000",
    TRUE ~ "0-100,000"))

print(table(tract$inc_category))

## 
##       0-100,000 100,000-200,000 200,000-300,000 
##             345             156              29

## tmap mode set to interactive viewing
tmap_mode("view")

tm_shape(tract) + tm_polygons(col = "hhincome")

tm_shape(tract) + tm_polygons(col = "inc_category")

# Reading the yelp data
yelp_hosi <- st_read("https://raw.githubusercontent.com/ujhwang/urban-analytics-2024/main/Assignment/mini_3/yelp_hospital.geojson")

## Reading layer `yelp_hospital' from data source 
##   `https://raw.githubusercontent.com/ujhwang/urban-analytics-2024/main/Assignment/mini_3/yelp_hospital.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 129 features and 23 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -84.56242 ymin: 33.60009 xmax: -84.08677 ymax: 34.0701
## Geodetic CRS:  WGS 84

skim(yelp_hosi)

tm_shape(yelp_hosi) + tm_dots(col="red") + tm_shape(tract) + tm_borders()

Appending Census data

1. make sure the CRS of the two sf objects are same
2. Calculate how many hospitals are within each census tract
3. Join the two files with the count of hospitals

## Check to see that the CRS for both the sf files are the same
head(tract$geometry)

## Geometry set for 6 features 
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -84.47194 ymin: 33.65339 xmax: -84.16538 ymax: 33.77463
## Geodetic CRS:  NAD83
## First 5 geometries:

## MULTIPOLYGON (((-84.18651 33.72737, -84.18381 3...

## MULTIPOLYGON (((-84.34921 33.76496, -84.34915 3...

## MULTIPOLYGON (((-84.34919 33.7481, -84.34919 33...

## MULTIPOLYGON (((-84.2965 33.76991, -84.29264 33...

## MULTIPOLYGON (((-84.4686 33.66347, -84.46678 33...

head(yelp_hosi$geometry)

## Geometry set for 6 features 
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -84.39409 ymin: 33.81016 xmax: -84.26324 ymax: 34.0701
## Geodetic CRS:  WGS 84
## First 5 geometries:

## POINT (-84.26324 34.0701)

## POINT (-84.27857 34.05106)

## POINT (-84.26324 34.0701)

## POINT (-84.28281 34.06029)

## POINT (-84.37021 33.84957)

tract <- tract %>% st_transform(crs=4326) 
yelp_hosi <- yelp_hosi %>% st_transform(crs=4326)

# Calculate the count of hospitals at the tract level
buffer_distance <- 402.336
tract_buffered <- tract %>%
  st_buffer(dist = buffer_distance)

intersections <- st_intersects(tract_buffered, yelp_hosi)

tract$buffered_hospital_count <- sapply(intersections, length)

# Join tract geometry with hospitals
tract_hosi <- st_join(tract, yelp_hosi %>% mutate(hos_count = 1))

tract_hosi_count <- tract_hosi %>%
  group_by(GEOID) %>%
  summarise(hos_count = sum(hos_count, na.rm = TRUE))

tract <- tract %>%
  left_join(tract_hosi_count %>% st_drop_geometry(), 
            by = "GEOID")

tract$hos_count[is.na(tract$hos_count)] <- 0

tract <- tract %>%
  mutate(total_hospital_count = hos_count + buffered_hospital_count)

skim(tract)

## Warning: Couldn't find skimmers for class: sfc_MULTIPOLYGON, sfc; No
## user-defined `sfl` provided. Falling back to `character`.

Variable type: character

Variable type: numeric

# Let's check to see if the points of yoga studios and their counts in the polygon data match
tm_shape(tract) + tm_polygons(col="hos_count") +
  tm_shape(yelp_hosi) + tm_dots()

Some questions

ggplot(tract, aes(x=hhincome, y=total_hospital_count)) +
  geom_point() +
  ylab("Number of Hospitals") +
  xlab("Household Income in Tract")

## Warning: Removed 11 rows containing missing values or values outside the scale range
## (`geom_point()`).

tm_shape(tract) + tm_polygons(col="inc_category") +
  tm_shape(yelp_hosi) +tm_dots()

boxplot(total_hospital_count ~ inc_category, data=tract, 
        main="Boxplot of Hospitals by Household Income Type", 
        xlab="Household Income types", 
        ylab="Hospitals count")

tract <- tract %>% 
  mutate(hosi_binary = total_hospital_count > 0)

hosi_binary_summary <- tract %>% 
  st_drop_geometry() %>% 
  group_by(hosi_binary, inc_category) %>% 
  summarize(n = n()) %>% 
  ungroup() %>% 
  mutate(proportion = n / sum(n))

## `summarise()` has grouped output by 'hosi_binary'. You can override using the
## `.groups` argument.

hosi_binary_summary %>% 
  ggplot(aes(x = inc_category, y = proportion, fill = hosi_binary)) +
  geom_bar(stat = 'identity', position = 'fill')