population mapping

2019 tract population

# pull 2019 ACS data
durham2019 <- get_acs(
  geography = "tract",
  variables = "B01001_001E",
  year = 2019,
  state = "NC",
  county = "Durham County",
  geometry = TRUE
)

nrow( durham2019 )
## [1] 60
durham2019 %>%
  summarize( mean( estimate ), sd( estimate ), sum( estimate )) %>%
  st_drop_geometry()
# create "tract_code" variable to identify tract numbers 
durham2019 <- durham2019 %>%
  select( -variable, -NAME ) %>%
  mutate( tract_code = substr( GEOID, 6, 11 )) %>%
  relocate( tract_code, .after = GEOID )

# create color palette for population values 
mypal2019 <- colorBin( palette = "YlGnBu",
                       domain = durham2019$estimate,
                       bins = 5,
                       pretty = TRUE )

# create popup content for the map
popup <- paste0( "tract: ", durham2019$tract_code, "<br>",
                 "population: ", durham2019$estimate )

# create bounding box for the map
bbox <- st_bbox( durham2019 ) %>%
  as.vector()

# create the leaflet map
mymap2019 <- leaflet( options = leafletOptions( zoomSnap = 0 )) %>%
  addProviderTiles( "CartoDB.Positron" ) %>%
  fitBounds( bbox[1], bbox[2], bbox[3], bbox[4] ) %>%
  addPolygons( data = durham2019,
               fillColor = ~mypal2019( estimate ), 
               color = "#b2aeae",
               fillOpacity = 0.7,
               weight = 1, 
               smoothFactor = 0.2,
               popup = popup ) %>%
  addLegend( pal = mypal2019, 
             values = durham2019$estimate,
             position = "bottomright",
             title = "population" )

# show map
mymap2019

In 2019, Durham County, North Carolina encompassed 60 census tracts with a total estimated population of 311,848 residents. The population distribution across tracts exhibited a mean of 5,197 inhabitants with a standard deviation of 2,455, reflecting notable demographic variation within the county.

2022 tract population

# pull 2022 ACS data
durham2022 <- get_acs(
  geography = "tract",
  variables = "B01001_001E",
  year = 2022,
  state = "NC",
  county = "Durham County",
  geometry = TRUE
)

nrow( durham2022 )
## [1] 68
durham2022 %>%
  summarize( mean( estimate ), sd( estimate ), sum( estimate )) %>%
  st_drop_geometry()
# create "tract_code" variable to identify tract numbers 
durham2022 <- durham2022 %>%
  select( -variable, -NAME ) %>%
  mutate( tract_code = substr( GEOID, 6, 11 )) %>%
  relocate( tract_code, .after = GEOID )

# create color palette for population values 
mypal2022 <- colorBin( palette = "YlGnBu",
                       domain = durham2022$estimate,
                       bins = 5,
                       pretty = TRUE )

# create popup content for the map
popup <- paste0( "tract: ", durham2022$tract_code, "<br>",
                 "population: ", durham2022$estimate )

# create bounding box for the map
bbox <- st_bbox( durham2022 ) %>%
  as.vector()

# create the leaflet map
mymap2022 <- leaflet( options = leafletOptions( zoomSnap = 0 )) %>%
  addProviderTiles( "CartoDB.Positron" ) %>%
  fitBounds( bbox[1], bbox[2], bbox[3], bbox[4] ) %>%
  addPolygons( data = durham2022,
               fillColor = ~mypal2022( estimate ), 
               color = "#b2aeae",
               fillOpacity = 0.7,
               weight = 1, 
               smoothFactor = 0.2,
               popup = popup ) %>%
  addLegend( pal = mypal2022, 
             values = durham2022$estimate,
             position = "bottomright",
             title = "population" )

# show map
mymap2022

Following significant population growth, Durham County expanded to 68 census tracts in 2022, with a total population of 325,101 residents. The mean population per tract decreased to 4,781, while maintaining comparable demographic dispersion with a standard deviation of 2,147.

population change

# create population change object 
durham_popchange <- durham2022 %>%
  inner_join(
    durham2019 %>%
      st_drop_geometry() %>%
      select( tract_code, estimate ) %>%
      rename( estimate2019 = estimate ), 
    by = "tract_code" ) %>%
  rename( estimate2022 = estimate ) %>%
  select( -moe ) %>%
  mutate( change = estimate2022 - estimate2019,
          percentchange = (change / estimate2019 ) * 100 )

# inspect data 
colnames( durham_popchange %>% st_drop_geometry())
## [1] "GEOID"         "tract_code"    "estimate2022"  "estimate2019" 
## [5] "change"        "percentchange"
head( durham_popchange %>% st_drop_geometry() )
durham_popchange %>%
  summarize( mean( change ), 
             quantile( change, 0.5 ), 
             quantile( change, 0.8 )) %>% 
  st_drop_geometry()

Over the three-year period, Durham County census tracts experienced a median population increase of 162 residents, with a mean growth of 164 inhabitants. The most rapidly expanding tracts, representing the top quintile, gained 577 or more residents during this interval.

# calculate mean and standard deviation
meanchange <- mean( durham_popchange$change )
sdchange <- sd( durham_popchange$change )

meanchange
## [1] 163.9615
sdchange
## [1] 672.0677
# create a histogram of the population changes 
hist( durham_popchange$change,
      probability = TRUE,
      main = "distribution of population change across census tracts",
      xlab = "population change",
      ylab = "density" )

# overlay normal curve 
curve( dnorm( x, mean = meanchange, sd = sdchange ), 
       add = TRUE, col = "aquamarine2", lwd = 2 )

# create color palette for change values 
mypalchange <- colorBin( palette = "YlGnBu",
                         domain = durham_popchange$change,
                         bins = 5, 
                         pretty = TRUE )

# create popup content for the map
popup <- paste0( "tract: ", durham_popchange$tract_code, "<br>",
                 "population change: ", durham_popchange$change, "<br>",
                 "percentage change: ", durham_popchange$percentchange)

# create bounding box for the map
bbox <- st_bbox( durham_popchange ) %>%
  as.vector()

# create the leaflet map
mymapchange <- leaflet( options = leafletOptions( zoomSnap = 0 )) %>%
  addProviderTiles( "CartoDB.Positron" ) %>%
  fitBounds( bbox[1], bbox[2], bbox[3], bbox[4] ) %>%
  addPolygons( data = durham_popchange,
               fillColor = ~mypalchange( change ), 
               color = "#b2aeae",
               fillOpacity = 0.7,
               weight = 1, 
               smoothFactor = 0.2,
               popup = popup ) %>%
  addLegend( pal = mypalchange, 
             values = durham_popchange$change,
             position = "bottomright",
             title = "population" )

# show map
mymapchange

Population change exhibited substantial variation across census tracts, with a standard deviation of 672 inhabitants. Several tracts lack comparative data due to their creation between the 2019 and 2022 American Community Survey periods, reflecting the county’s dynamic demographic shifts and necessitating boundary adjustments.

tract changes between 2019 and 2022

# find tracts new to the 2022 ACS
newtracts2022 <- durham2022 %>%
  filter( !GEOID %in% durham2019$GEOID )

newtracts2022
# find tracts removed after the 2019 ACS
removedtracts2019 <- durham2019 %>%
  filter( !GEOID %in% durham2022$GEOID )

removedtracts2019

tract changes maps, 2019 and 2022 side by side

# create 2019 map with removed tracts filled in
removedmap2019 <- leaflet( durham2019 ) %>%
  addProviderTiles( "CartoDB.Positron" ) %>%
  addPolygons( fillColor = "#FFFFFF",
               color = "#808080",
               weight = 1) %>%
  addPolygons( data = removedtracts2019,
               fillColor = "cadetblue3",
               popup = ~tract_code,
               fillOpacity = 0.7 )

# create 2022 map with new tracts filled in
newmap2022 <- leaflet( durham2022 ) %>%
  addProviderTiles( "CartoDB.Positron" ) %>%
  addPolygons( fillColor = "#FFFFFF",
               color = "#808080",      # border color
               weight = 1) %>%  
  addPolygons( data = newtracts2022,
               fillColor = "cadetblue3",
               popup = ~tract_code,
               fillOpacity = 0.7 )

# create patchwork map sidebyside
tagList(
  tags$div( style = "display: flex;",
    tags$div( style = "width: 50%;", removedmap2019 ),
    tags$div( style = "width: 50%;", newmap2022 )
  )
)

tract matching & combination

# match new 2022 census tracts to removed 2019 census tracts
newcensustracts2022 <- data.frame(
  removedtract_code = c( "002027", "002018", "001502", "001707", 
                    "002017", "001807", "002028", "002016" ), 
  newGEOID1 = c( "002036", "002033", "001504", "001712",
                 "002029", "001810", "002037", "002031" ), 
  newGEOID2 = c( "002035", "002034", "001505", "001713",
                 "002030", "001811", "002038", "002032" )
)

# combine new census tracts to get a total comparable to the old census tract
result <- newcensustracts2022 %>%
  left_join( durham2022 %>%
               select( tract_code, estimate ) %>%
               rename( estimate1 = estimate ),
             by = c( "newGEOID1" = "tract_code" )) %>%
  left_join( durham2022 %>%
               select( tract_code, estimate ) %>%
               rename( estimate2 = estimate ), 
             by = c( "newGEOID2" = "tract_code" )) %>%
  mutate( newtotal = estimate1 + estimate2 ) %>%
  select( removedtract_code, newtotal )

# view results
result 
# update durham_popchange to include all 2019 tracts  
durham_popchange_new <- durham2019 %>%
  left_join(
    durham2022 %>%
      st_drop_geometry() %>%
      select( tract_code, estimate ) %>%
      rename( estimate2022 = estimate ), 
    by = "tract_code" ) %>%
  rename( estimate2019 = estimate ) %>%
  select( -moe ) %>%
  mutate( change = estimate2022 - estimate2019,
          percentchange = (change / estimate2019 ) * 100 )

# update durham_popchange_new with the new combined estimates
durham_popchange_new <- durham_popchange_new %>%
  left_join( result, by = c( "tract_code" = "removedtract_code" )) %>%
  mutate( estimate2022 = coalesce( estimate2022, newtotal )) %>%
  select( -newtotal ) %>% 
  mutate( change = estimate2022 - estimate2019,
          percentchange = (change / estimate2019 ) * 100 )

colnames( durham_popchange_new )
## [1] "GEOID"         "tract_code"    "estimate2019"  "estimate2022" 
## [5] "change"        "percentchange" "geometry"
durham_popchange_new

updated population change

# calculate mean and standard deviation
meanchange_new <- mean( durham_popchange_new$change )
sdchange_new <- sd( durham_popchange_new$change )

meanchange_new
## [1] 220.8833
sdchange_new
## [1] 735.5108
# create a histogram of the population changes 
hist( durham_popchange_new$change,
      probability = TRUE,
      main = "distribution of population change across census tracts",
      xlab = "population change",
      ylab = "density" )

# overlay normal curve 
curve( dnorm( x, mean = meanchange_new, sd = sdchange_new ), 
       add = TRUE, col = "aquamarine2", lwd = 2 )

# recreate color palette for change values 
mypalchange <- colorBin( palette = "YlGnBu",
                         domain = durham_popchange_new$change,
                         bins = 5, 
                         pretty = TRUE )

# recreate popup content for the map
popup <- paste0( "tract: ", durham_popchange_new$tract_code, "<br>",
                 "population change: ", durham_popchange_new$change, "<br>",
                 "percentage change: ", durham_popchange_new$percentchange)

# recreate bounding box for the map
bbox <- st_bbox( durham_popchange_new ) %>%
  as.vector()

# recreate the leaflet map
mymapchange_new <- leaflet( options = leafletOptions( zoomSnap = 0 )) %>%
  addProviderTiles( "CartoDB.Positron" ) %>%
  fitBounds( bbox[1], bbox[2], bbox[3], bbox[4] ) %>%
  addPolygons( data = durham_popchange_new,
               fillColor = ~mypalchange( change ), 
               color = "#b2aeae",
               fillOpacity = 0.7,
               weight = 1, 
               smoothFactor = 0.2,
               popup = popup ) %>%
  addLegend( pal = mypalchange, 
             values = durham_popchange_new$change,
             position = "bottomright",
             title = "population" )

# show map
mymapchange_new

After accounting for the redistricted census tracts, the mean population growth increased to 221 residents per tract, with a standard deviation of 736. This substantial variance suggests the redistricting primarily occurred in areas experiencing significant population expansion.

# calculate total population estimates for 2019 and 2022
totalest2019 <- sum( durham_popchange_new$estimate2019 )
totalest2022 <- sum( durham_popchange_new$estimate2022 )

totalest2019
## [1] 311848
totalest2022
## [1] 325101
# calculate percentage change from 2019 to 2022

totalpercentchange <- (( totalest2022 - totalest2019 ) / totalest2019 ) * 100
totalpercentchange
## [1] 4.249827
uspopulationest2019 <- 328239523
uspopulationest2022 <- 332403650

uspercentchange <- (( uspopulationest2022 - uspopulationest2019 ) / 
                      uspopulationest2019 ) * 100
uspercentchange
## [1] 1.268624

During the three-year period from 2019 to 2022, Durham’s population grew at more than three times the national rate, with a 4.25% increase compared to the U.S. growth rate of 1.27%.

This significant population growth rate warrants a comprehensive review of Durham EMS resource allocation. The data suggests the need to evaluate current unit deployment and station locations, with consideration for expanding capacity through additional units and facilities to maintain effective emergency response coverage.

EMS data and mapping

ems_data <- read_excel( "~/Documents/Github/capstone-project/ems_data.xlsx" )

colnames( ems_data )
##  [1] "date"          "time"          "latitude"      "longitude"    
##  [5] "call_type"     "response_unit" "home_station"  "response_time"
##  [9] "action_taken"  "patient_age"   "patient_race"  "patient_sex"  
## [13] "hospital"
head( ems_data )
# set durham city center in latitude and longitude coordinates
durham_center <- data.frame(
  lon = -78.8986, 
  lat = 35.9940
) %>% 
  st_as_sf(coords = c("lon", "lat"), crs = 4326)

# load the 8 existing EMS station locations
ems_stations <- data.frame(
  address = c(
    "2725 Holloway St, Durham, NC 27703",
    "615 Old Fayetteville St, Durham, NC 27701",
    "2715 S Alston Ave, Durham, NC 27713",
    "2910 Latta Rd, Durham, NC 27712",
    "2212 Chapel Hill Rd, Durham, NC 27707",
    "226 Milton Rd, Durham, NC 27712",
    "402 Stadium Dr, Durham, NC 27704",
    "225 Infinity Rd, Durham, NC 27712"
  )
)

# convert addresses to coordinates 
ems_locations <- ems_stations %>%
  geocode(address = address, method = "osm") %>%
  st_as_sf(coords = c("long", "lat"), crs = 4326)

ems_locations
# create a map of population change with the ems stations added 
emslocationsmap <- mymapchange_new %>%
  addCircleMarkers( 
   data = ems_locations,
   radius = 5,
   color = "red", 
   fillColor = "red",
   fillOpacity = 0.7,
   weight = 2,
   popup = ~address
)

emslocationsmap
# create Voronoi polygons
voronoi <- st_voronoi( do.call( c, st_geometry( ems_locations ))) %>%
  st_collection_extract() %>%
  st_sf( crs = 4326 ) %>%
  st_transform( st_crs( durham2022 )) %>%
  st_intersection( st_union( durham2022 ))

# add to map
emslocationsmap <- emslocationsmap %>%
  addPolygons( data = voronoi,
               color = "purple",
               fillColor = "purple",
               fillOpacity = 0.2,
               weight = 2,
               opacity = 0.8
 )

emslocationsmap
# convert ems_data to sf object with correct CRS
ems_data_sf <- ems_data %>%
  st_as_sf( coords = c( "longitude", "latitude" ), crs = 4326 ) %>%
  st_transform( st_crs( durham2022 ))

# filter points within Durham boundaries
ems_data_filtered <- ems_data_sf %>%
 st_intersection( st_union( durham2022 ))

# create new color palette 
mypalcalls <- colorFactor( palette = "YlGnBu",
                           domain = ems_data$response_unit )

# create bounding box for the map
bbox <- st_bbox( durham2022 ) %>%
  as.vector()

# create the leaflet map
emscallmap <- leaflet( options = leafletOptions( zoomSnap = 0 )) %>%
  addProviderTiles( "CartoDB.Positron" ) %>%
  fitBounds( bbox[1], bbox[2], bbox[3], bbox[4] ) %>%
  addCircleMarkers( data = ems_data_filtered,
                    radius = 2,
                    fillColor = ~mypalcalls( response_unit ),
                    color = ~mypalcalls( response_unit ),
                    fillOpacity = 0.7,
                    weight = 1 ) %>%
  addLegend(
    pal = mypalcalls,
    values = ems_data_filtered$response_unit,
    title = "response unit"
  ) %>%
  addPolygons( data = voronoi,
               color = "purple",
               fillColor = "purple",
               fillOpacity = 0.2,
               weight = 2,
               opacity = 0.8 )

# show map
emscallmap