# 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.
# 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.
# 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.
# 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
# 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 )
)
)
# 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
# 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 <- 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