Introduction
The heat wave in Colorado a few weeks ago motivated me to think about access to air conditioning for people who cannot afford it at their own home or apartment. During heat waves Denver Denver opens up its recreation centers as cooling centers. I decided to use geospatial analysis to try to figure out how many people in Denver have access to a cooling center.
Goal
- The goal of this analysis is to identify locations/areas/populations in Denver that have lower or no access to cooling centers, in order to inform placement of additional cooling resources for these populations.
Methods Overview:
- Obtain and map Denver cooling center locations.
- Compute isochrones (the area that can be reached by walking/driving in a certain time) from each cooling center.
- Identity areas that do not have access to cooling centers, defined here as not being inside any of the isochrones.
Data
Denver Cooling Centers
I found a data set containing the location of all Denver recreation centers on the Denver Open Data Catalog.
The Denver recreation centers shapefile was read into R with the {sf} package (Pebesma and Bivand 2023).
Load cooling centers
# rec shape file crs = 2877 for Colorado; use for calculations later?
shp_file <- "./data/ODC_PARK_RECCENTER_P_-1008034318397091855/PARK_RECCENTER_P.shp"
rec_shp <- sf::read_sf(shp_file) |>
sf::st_transform(4326) |>
janitor::clean_names()
head(rec_shp)Simple feature collection with 6 features and 26 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: -105.0512 ymin: 39.7582 xmax: -104.9832 ymax: 39.78888
Geodetic CRS: WGS 84
# A tibble: 6 × 27
rec_name loc_code address_id address_li address_1 city state zip phone
<chr> <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <chr>
1 Scheitler Re… RB2 0 5031 W 46… <Null> Denv… CO <NA> 720-…
2 Aztlan Recre… RA3 0 4435 Nava… <NA> Denv… CO <NA> 720-…
3 Highland Rec… SH1 0 2880 Osce… <NA> Denv… CO <NA> 720-…
4 Ashland Recr… RA1 0 2475 W Du… <NA> Denv… CO <NA> 720-…
5 5090 Broadwa… RS2 0 5090 Broa… <NA> Denv… CO <NA> 720-…
6 Globeville R… RG2 0 4496 Gran… <NA> Denv… CO <NA> 303-…
# ℹ 18 more variables: rec_type <chr>, marketing <chr>, marketed_f <chr>,
# marketed_p <chr>, programs_l <chr>, pool_hours <chr>, news_link <chr>,
# hours <chr>, photo <chr>, year_built <chr>, year_remod <chr>,
# bldg_sqft <dbl>, label <chr>, latitude <dbl>, longitude <dbl>,
# facilities <chr>, globalid <chr>, geometry <POINT [°]>County boundary and census data
To calculate the population within/outside range of cooling centers, I decided to use census tracts from the US census. I obtained shapefiles for census tracts and their population using the {tidycensus} R package (Walker and Herman 2024). The shapefile for the entire county (without population) was downloaded using the {tigris} package (Walker 2024b).
Get census tracts
# Data from 2020 decennial US Census
den_tracts <- get_decennial(
geography = "tract",
variables = "P1_001N",
year = 2020,
state = "CO",
county = "Denver",
geometry = TRUE
) |> st_transform(4326)There are 178 tracts, with a total population of 715,522.
The data sets are visualized on an interactive map (figure Figure 1) using the {leaflet} package (Cheng et al. 2024).
Plot map of data
leaflet() |>
addTiles() |>
addMarkers(data = rec_shp,
label = ~rec_name) |>
addPolygons(data = den_tracts,
fillOpacity = 0.2,
color = "black",
weight = 1,
label = ~NAME,
highlight = highlightOptions(
weight = 3,
fillOpacity = 0,
color = "black",
opacity = 1.0,
bringToFront = TRUE,
sendToBack = TRUE)
)Analysis
Compute isochrones
I used the {mapboxapi} package (Walker 2024a) to calculate walking isochrones from each recreation center location. I chose a 20 minute walk as reasonable access.
Calculate isochrones
time_minutes <- 20
isos <- mapboxapi::mb_isochrone(
location = rec_shp,
profile = "walking",
time = time_minutes,
id_column = "rec_name"
)Figure 2 is an interactive map that shows all of the recreation centers and isochrones, as well as census tracts.
Map isochrones
leaflet() |>
addProviderTiles(providers$CartoDB.Voyager) |>
addPolygons(data = den_tracts, fillOpacity = 0.2, color = "black", weight = 1, label = ~NAME) |>
addPolygons(data = isos, label = ~id) |>
addMarkers(data = rec_shp,
label = ~rec_name) Calculate Access to Cooling Centers
The next step was to calculate what area and population are within/outside these isochrones. I tried two methods: one is simpler, and one is more complex but probably more accurate.
Simpler method
The simpler method is to keep census tracts whole and find which ones intersect any of the isochrones:
- Calculate union of all the isochrones (using st_union() from the {sf} package)
- Find which census tracts do not intersect isochrones at all (using st_join() from {sf})
Method 1
# Union of all isochrones
isos_union <- sf::st_union(isos, by_feature = FALSE) |> sf::st_as_sf()
isos_union$index <- 99 # need to add dummy column to identify where join is NA?
# join tracts and isocrhoes
tracts_iso_join <- sf::st_join(den_tracts, isos_union)
# find which tracts were not in isochrones (join is NA)
not_in <- tracts_iso_join |> filter(is.na(index))
#glue::glue('There are {nrow(not_in)} census tracts in Denver (out of total {nrow(den_tracts)}) that are not in range of a cooling center. The population of census tracts without access to cooling centers is {sum(not_in$value)}.')There are 36 census tracts in Denver (out of a total 178) that are not in range of a cooling center. The population of census tracts without access to cooling centers is 135,635.
Figure 3 shows the results of this method, with census tracts that do not intersect isochrones shaded in red.
Map method 1
leaflet() |>
addProviderTiles(providers$CartoDB.Voyager) |>
addPolygons(data = den_tracts, fillOpacity = 0.2, color = "black", weight = 1, label = ~NAME) |>
addPolygons(data = not_in, fillOpacity = 0.2, color = "red", weight = 1, label = ~NAME) |>
addPolygons(data = isos_union, color = "blue", weight = 2, label = "Isochrones")# |>Weighting by Area
The first method is much simpler but probably not very accurate, since some of the census tracts just barely intersect an isochrone. A more complex but more accurate method is to use areal interpolation that considers what part of each census tract intersects the isochrones. Note that this method assumes that population is distributed evenly across each census tract.
- Calculate union of all the isochrones
- For each tract, compute the percent of area that intersects with isochrones (using st_intersection() from {sf})
- Use that ratio to multiply the total population of the tract, to find population in/out.
- Sum the population estimates from all the tracts.
Method 2
# function to compute perc overlap and weighted population for 1 tract
calc_pop_out <- function(tract, isos) {
tract_intersect <- sf::st_intersection(tract, isos)
total_pop <- tract$value
if ( nrow(tract_intersect) == 0) {return(total_pop)}
area_tot <- sf::st_area(tract)
area_intersect <- sf::st_area(tract_intersect)
perc_intersect <- as.numeric(area_intersect/area_tot*100)
perc_area_not_in <- 100 - perc_intersect
pop_out <- perc_area_not_in/100*total_pop
}
# apply to each tract (could use *map* function instead of loop? )
pops_out <- rep(NA, times = nrow(den_tracts))
for (i in 1:nrow(den_tracts)) {
pops_out[i] <- calc_pop_out(den_tracts[i,], isos_union)
}
# Union of all Denver census tracts >> 1 multi-polygon
tracts_union <- sf::st_union(den_tracts)
tracts_isos_intersect <- sf::st_intersection(tracts_union, isos_union)
tracts_isos_diff <- sf::st_difference(tracts_union, isos_union)
area_tot <- sf::st_area(tracts_union)
area_in <- sf::st_area(tracts_isos_intersect)
area_out <- sf::st_area(tracts_isos_diff)
percent_out <- round(area_out/area_tot*100)
pop_out <- round(sum(pops_out))
#glue('Approximately {round(area_out/area_tot*100)} percent of the area of Denver is not within a 20 minute walk of a cooling center.')
#glue('An estimated {round(sum(pops_out))} people in Denver are not within a 20 minute walk of a cooling center.')Figure Figure 4 shows the results of this method. The union of all the isochrones is shown in blue, and the difference (the area not within any isochrones) is shown in blue.
Map method 2
leaflet() |>
addProviderTiles(providers$CartoDB.Voyager) |>
addPolygons(data = tracts_isos_intersect,
color = "blue",
weight = 1,
label = "Intersection") |>
addPolygons(data = tracts_isos_diff,
color = "red",
weight = 1,
label = "Difference")Conclusion
- 20-minute walking isochrones were calculated for all Denver cooling centers (recreation centers).
- The area and population not within any of these isochrones was computed.
- Approximately 69 percent of the area of Denver (County) is not within a 20 minute walk of a cooling center.
- An estimated 389,095 people in Denver (County) are not within a 20 minute walk of a cooling center.
Caveats/future work
This analysis provides a starting framework and estimate; however there are many assumptions and areas for further research and improvement. Below are some that come to mind:
- Areal interpolation assumes that population is evenly distributed across each tract area.
- Future work could try to incorporate estimates of how many people have air conditioning in their homes or apartments.
- Future work could attempt to also include other cooling locations besides recreation centers (such as public libraries).
- It is difficult to specify exactly what “access” to a cooling center means (here I defined it as a 20 minute walk). Future work could also try to take into account public transit (though currently the Mapbox isochrone API does not appear to have a public transit mode option).
- It would be interesting to compare the distribution of cooling centers with data on tree coverage and temperatures.
SessionInfo
R version 4.4.1 (2024-06-14)
Platform: x86_64-apple-darwin20
Running under: macOS Sonoma 14.6.1
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-x86_64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: America/Denver
tzcode source: internal
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] tigris_2.1 tidycensus_1.6.5 glue_1.7.0 janitor_2.2.0
[5] dplyr_1.1.4 sf_1.0-17 leaflet_2.2.2 mapboxapi_0.6
loaded via a namespace (and not attached):
[1] xfun_0.47 raster_3.6-26 htmlwidgets_1.6.4
[4] lattice_0.22-6 tzdb_0.4.0 leaflet.providers_2.0.0
[7] vctrs_0.6.5 tools_4.4.1 crosstalk_1.2.1
[10] generics_0.1.3 curl_5.2.1 tibble_3.2.1
[13] proxy_0.4-27 fansi_1.0.6 pkgconfig_2.0.3
[16] KernSmooth_2.23-24 protolite_2.3.0 uuid_1.2-0
[19] lifecycle_1.0.4 compiler_4.4.1 stringr_1.5.1
[22] munsell_0.5.1 terra_1.7-71 codetools_0.2-20
[25] snakecase_0.11.1 htmltools_0.5.8.1 class_7.3-22
[28] yaml_2.3.10 jquerylib_0.1.4 crayon_1.5.3
[31] pillar_1.9.0 tidyr_1.3.1 aws.s3_0.3.21
[34] classInt_0.4-10 slippymath_0.3.1 wk_0.9.1
[37] magick_2.8.4 tidyselect_1.2.1 rvest_1.0.4
[40] digest_0.6.36 stringi_1.8.4 purrr_1.0.2
[43] fastmap_1.2.0 grid_4.4.1 colorspace_2.1-0
[46] cli_3.6.3 magrittr_2.0.3 base64enc_0.1-3
[49] utf8_1.2.4 e1071_1.7-14 aws.signature_0.6.0
[52] withr_3.0.1 readr_2.1.5 scales_1.3.0
[55] rappdirs_0.3.3 sp_2.1-4 lubridate_1.9.3
[58] timechange_0.3.0 rmarkdown_2.28 httr_1.4.7
[61] jpeg_0.1-10 hms_1.1.3 png_0.1-8
[64] evaluate_0.24.0 knitr_1.48 s2_1.1.6
[67] rlang_1.1.4 Rcpp_1.0.13 DBI_1.2.2
[70] geojsonsf_2.0.3 xml2_1.3.6 renv_1.0.4
[73] rstudioapi_0.16.0 jsonlite_1.8.8 R6_2.5.1
[76] units_0.8-5 Updates
- 2024-08-19 : Update to use 2020 census data instead of 2010. Total population increases from 600,158 to 715,522.
- 2024-08-19 : Fix error in calculating population not within isochrones. For tracts that do not intersect at all, population was set to 0 instead of tract population.