# install
install.packages("tidyverse") # cleaning data
install.packages("janitor") # cleaning data
install.packages("tigris") # downloading Census boundaries
install.packages("sf") # working with spatial data
install.packages("spdep") # defining spatial neighbors
install.packages("bisdep") # bivariate analysis
install.packages("spam") # bivariate analysis
install.packages("osrm") # creating isochrones
install.packages("ggplot2") # plotting
install.packages("patchwork") # combining plotsMeasuring Park Access in King County, WA
Introduction
The following vignette provides an example of the spatial isochrone analysis used to measure park access in King County, Washington. The data used are publicly available and drawn from the King County Geographic Information System (GIS) Open Data website and the United States Census Bureau.
To build isochrones, the R package osrm (“Open Source Routing Machine”) is used. The osrm package draws on publicly available OpenStreetMap data to calculate routes, distances, and travel times based on road networks. Importantly, travel time calculations account for street network configuration and topographic constraints (e.g., highways and lakes).
R Packages
You will need the following R packages. If you have already installed them, skip the install.packages lines.
# load
library(tidyverse)
library(tidycensus)
library(janitor)
library(tigris)
library(sf)
library(spdep)
library(bispdep)
library(spam)
library(osrm)
library(ggplot2)
library(patchwork)Data Preparation
King County Parks
To begin, download the Parks in King County spatial layer (select Download and then Shapefile Download). Once downloaded, unzip the file and store it in your project folder.
Note
The unzipped file has been renamed to kc_parks for simplicity.
Before importing a spatial layer, it is good practice to first check the name of the layer with st_layers.
st_layers("data/kc_parks")$name[1] "Parks_in_King_County"This layer name is then entered into read_sf.
Note
Layer names can change depending on changes made by the data host/provider.
Now read in the spatial data using the read_sf function.
## load
parks <- read_sf("data/kc_parks", # shapefile location
"Parks_in_King_County") # layer nameAnd clean the data using tidyverse and janitor functions.
## clean
parks <- parks |>
# set CRS
st_transform("EPSG:2926") |>
# clean variable names
clean_names() |>
# filter to parks (removing pools, trails, and shopping sites)
filter(sitetype == "Park Site") |>
# filter out parks managed above county level
filter(!manager %in% c("US Forest Service",
"Washington State DNR",
"Other Federal",
"Washington State Parks",
"Army Corps of Engineers",
"US National Park Service")) |>
# select variables (and rename a few)
select(park_id = kc_fac_fid,
kcparkfid,
name = sitename)
Note
Setting the Coordinate Reference System (CRS) is an important but location-dependent step. EPSG:2926 is used here because it best reflects King County, WA. If you are working in a different region, search for an appropriate CRS at epsg.io or spatialreference.org by entering your location or projection name.
Census Block Groups
Census block group boundaries are easily accessed using the tigris package and block_groups function.
kc_bg <- block_groups(state = "WA",
county = "King",
year = 2023,
progress_bar = FALSE) |>
st_transform("EPSG:2926") |>
clean_names() |>
# erase water for cleaner mapping later
erase_water() |>
select(geoid, geometry)
Note
Be sure to adjust this code to reflect the Census administrative unit and geographic region you are working with. In addition to the block_groups function, tigris offers functions for states, counties, tracts, and blocks.
A Census API key may be required to access Census data. If you do not have one, request here (free). If you already have a key, install with tidycensus::census_api_key("YOUR_KEY_HERE", install = TRUE)
optional: click here to plot parks and block groups
When working with spatial data, producing plots as you go is a useful way to make sure your data are loading properly.
kc_bg |>
ggplot() +
geom_sf(fill = NA, # hollow block groups
color = "black",
linewidth = 0.01) +
# add parks
geom_sf(data = parks,
aes(fill = "Parks"),
color = NA, # no border
linewidth = 0.2) +
labs(title = "King County, WA",
subtitle = "Parks & Census Block Groups") +
# legend
scale_fill_manual(name = "", # leave blank
values = c("Parks" = "darkgreen")) +
theme_void() +
theme(plot.title = element_text(face = "bold",
hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
legend.text = element_text(size = 15),
legend.box.margin = margin(0, 0, 0, -50), # force legend to left for cleaner look
legend.key.size = unit(0.5, 'inches')) # increase legend size
Isochrones
With the data prepared, isochrones can be computed to determine park access. For visual simplicity, the following example uses the subset of King County block groups that comprise Vashon Island (vashon_bg).
# extract Vashon Island block groups
vashon_bg <- kc_bg |>
filter(grepl(5303302770, geoid))Vashon Island provides a compact, manageable area to demonstrate the workflow while avoiding long computation times. The full sample of King County block groups (n = 1,544) is run analogously but takes several hours to complete.
optional: click here to plot Vashon Island and King County
vashon_bg |>
ggplot() +
geom_sf(aes(fill = "Vashon Island"),
color = NA,
linewidth = 0.01) +
# add full block groups
geom_sf(data = kc_bg,
fill = NA, # hollow block groups
color = "black",
linewidth = 0.01) +
# add parks
geom_sf(data = parks,
aes(fill = "Parks"),
color = NA, # no border
linewidth = 0.2) +
labs(title = "King County, WA",
subtitle = "Parks & Census Block Groups") +
# legend
scale_fill_manual(name = "", # leave blank
values = c("Parks" = "darkgreen",
"Vashon Island" = "pink")) +
theme_void() +
theme(plot.title = element_text(face = "bold",
hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
legend.text = element_text(size = 12),
legend.box.margin = margin(0, 0, 0, -50), # force legend to left for cleaner look
legend.key.size = unit(0.3, 'inches')) # increase legend size
Generating isochrones takes four steps:
- Identify block group centroids (isochrones require a point location)
- Run
osrmIsochronefunction - Overlay parks layer
- Identify overlap
Step 1: Identify Centroids
vashon_centroids <- vashon_bg |>
mutate(geometry = st_centroid(geometry))Step 2: Run Isochrone Function
Isochrones are generated using the osrmIsochrone function. Given a starting point location, osrmIsochrone computes the area reachable within a certain travel time by drawing on actual road networks. The result is an irregularly shaped polygon, reflecting the fact that available road networks and topography will restrict travel distances in certain directions.
Here, you can specify a travel time of 30 minutes (breaks = 30) and the mode of travel as walking (osrm.profile = "foot"). The loc argument is the starting point location, to which you iteratively feed in the centroid of each block group. The n argument specifies the resolution of the resultant isochrone polygons; higher values for n, such as 500, produce higher quality isochrones but take longer to run.
As noted, you will iteratively loop through each block group. The osrmIsochrone function evaluates isochrones one at a time. The code below takes the list of block groups (temp), grabs a single block group (by referencing its geoid), and enters it into the osrmIsochrone function. Once an isochrone is built (i.e., a polygon spatial feature), it is stored in empty_list and the loop returns to the start of the function to run the next block group. Once all block groups are run, empty_list is bound together (bind_rows(empty_list)) providing a complete set of isochrones for each of the block groups. As a reminder, this process is identical to the full King County analysis, with the only difference being the time required to complete.
# transform CRS to WGS84 (necessary for `osrmIsochrone` function)
temp <- st_transform(vashon_centroids, "EPSG:4326")
empty_list <- list()
for (i in temp$geoid) {
tryCatch({
iso <- osrmIsochrone(loc = temp |> filter(geoid == i),
breaks = 30,
n = 500,
osrm.profile = "foot") |>
mutate(geoid = i) |>
select(geoid, geometry)
# add to list
empty_list[[as.character(i)]] <- iso
}, error = function(e) {
message(sprintf("Error at geoid %s: %s", i, e$message))
# add NULL to keep index consistent
empty_list[[as.character(i)]] <- NULL
})
}
# bind together
vashon_iso <- bind_rows(empty_list)
# return spatial feature to desired CRS
vashon_iso <- vashon_iso |>
st_as_sf() |>
st_transform("EPSG:2926")
optional: click here for notes on tryCatch function
As an added measure of precaution, the osrmIsochrone loop above is wrapped in a tryCatch function. tryCatch allows each loop to run without interruption, regardless of whether an error occurs. The osrmIsochrone function depends on external routing services (OSRM API), which may intermittently fail due to network timeouts, server load, or invalid routing queries.
If an error occurs, tryCatch prints a message indicating the geoid of the block group that did not generate an isochrone (message(sprintf("Error at geoid %s: %s", i, e$message))) and assigns a NULL value to that observation (empty_list[[as.character(i)]] <- NULL). This ensures that computation continues without interruption and that all successful isochrones are retained for downstream analysis.
This process allows failed osrmIsochrone requests to be explicitly identified and inspected after the initial run. These cases can then be reprocessed individually in subsequent runs, at which point they typically complete successfully and are appended to the final dataset. Importantly, these failures are attributable to intermittent limitations of the OSRM service rather than issues with the input data.
In brief, tryCatch allows you to identify errors without losing the isochrones that were already built. This approach is particularly valuable given the computational cost (approximately 4-6 hours) of generating isochrones for the full sample (1,544 block groups), where a single failure without error handling would otherwise terminate the entire process.
optional: click here for osrmIsochrone troubleshooting
osrmIsochrone requires a connection to an OSRM routing server, which performs network-based travel calculations. On some local machines (especially Windows PCs), you may encounter SSL/TLS errors that prevent R from connecting to the public OSRM servers. To run the code successfully, you will need either a machine with internet access that can reach a public OSRM server or a locally running OSRM server.
If you are having trouble running osrmIsochrone, the simplest solution is trying on a different computer or through a different internet connection. However, if you continue to encounter trouble, here are a few helpful sources of information:
Step 3: Overlay Parks
optional: click here to view code for plot above
First, create a new spatial layer with st_intersection that includes only parks on Vashon Island. st_intersection returns the overlap (i.e., intersection) of two spatial layers. Here, that is the full parks layer (parks) and Vashon Island block groups (vashon_bg).
vashon_parks <- st_intersection(vashon_bg, parks) |>
distinct(name, .keep_all = TRUE) |>
select(name, geometry)Second, plot the Vashon Island block groups (vashon_bg) and overlay the Vashon Island parks (vashon_parks).
vashon_bg |>
ggplot() +
geom_sf(fill = NA,
color = "black",
linewidth = 0.01) +
# add Vashon Island parks
geom_sf(data = vashon_parks,
aes(fill = "Parks"),
color = NA, # no border
linewidth = 0.2) +
labs(title = "Vashon Island\nParks & Block Groups") +
# legend
scale_fill_manual(name = "", # leave blank
values = c("Parks" = "darkgreen")) +
theme_void() +
theme(plot.title = element_text(face = "bold",
hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
legend.text = element_text(size = 12),
legend.box.margin = margin(-50, 0, 0, -35), # force legend to left and up for cleaner look
legend.key.size = unit(0.3, 'inches')) # increase legend sizeStep 4: Identify Park & Isochrone Overlap
Finally, generate a value for park access. “Park access” is operationalized as the number of park polygons intersecting a block group’s isochrone. To generate this value, use the st_intersection function, which returns the overlap (i.e., intersection) of two spatial layers.
iso_parks <- st_intersection(vashon_iso, parks) |>
select(name, geometry) |>
# return to an sf object with correct CRS
st_as_sf() |>
st_transform("EPSG:2926")Here, the parks (parks) and Vashon Island isochrone (vashon_iso) layers are supplied to st_intersection, returning the total number of parks that intersect each isochrone.
A visualization of these four steps is provided below, using a single isochrone to demonstrate the process. Panel D shows the parks (highlighted in yellow) that intersect by the isochrone.
optional: click here to view code for side-by-side plot above
iso_parks <- st_intersection(
# filter to a single isochrone (for demonstration)
vashon_iso |> filter(geoid == 530330277014),
vashon_parks
) |>
select(name, geometry) |>
# return to an sf object with correct CRS
st_as_sf() |>
st_transform("EPSG:2926")
# vashon block groups
v1 <- vashon_bg |>
ggplot() +
geom_sf(fill = NA) +
geom_sf(data = vashon_centroids,
size = 0.5) +
theme_void() +
labs(title = "Block Group\nCentroids",
caption = "Panel A") +
theme(plot.title = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0.5,
size = 12))
# vashon isochrone
v2 <- vashon_bg |>
ggplot() +
geom_sf(fill = NA) +
geom_sf(data = vashon_centroids,
size = 0.5) +
geom_sf(data = vashon_iso |> filter(geoid == 530330277014),
color = '#e7298a',
fill = NA,
linewidth = 0.5) +
theme_void() +
labs(title = "30 Minute\nIsochrone",
caption = "Panel B") +
theme(plot.title = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0.5,
size = 12))
# vashon parks
v3 <- vashon_bg |>
ggplot() +
geom_sf(fill = NA) +
geom_sf(data = vashon_parks,
fill = 'darkgreen') +
geom_sf(data = vashon_iso |> filter(geoid == 530330277014),
color = '#e7298a',
fill = NA,
linewidth = 0.5) +
theme_void() +
labs(title = "Public\nParks",
subtitle = paste0("n = ", nrow(vashon_parks)),
caption = "Panel C") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0.5,
size = 12))
# vashon parks in isochrone
v4 <- vashon_bg |>
ggplot() +
geom_sf(fill = NA) +
geom_sf(data = vashon_parks,
fill = 'darkgreen') +
geom_sf(data = vashon_iso |> filter(geoid == 530330277014),
color = '#e7298a',
fill = NA,
linewidth = 0.5) +
geom_sf(data = iso_parks,
fill = 'yellow') +
theme_void() +
labs(title = "Parks in\nIsochrone",
subtitle = paste0("n = ", nrow(iso_parks)),
caption = "Panel D") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0.5,
size = 12))
# side-by-side plot
(v1 | v2 | v3 | v4)Park Access
Once isochrones are built, you can assign a park access value to each block group (i.e., perform Step 4 across all block groups) by using the st_buffer and st_intersects functions. (Above, only st_intersection was needed because you were only visualizing the overlap.)
st_buffer builds a buffer around a spatial object. The spatial object supplied will be the isochrones. Importantly, st_buffer requires a distance value (dist) that reflects the distance the buffer should extend from the spatial object. Because the size of each isochrone polygon already reflects that desired distance, set dist to zero.
buffer <- st_buffer(vashon_iso, dist = 0)st_intersects is a Boolean version of st_intersection, returning TRUE if an intersection is present and FALSE otherwise. Because each isochrone will serve as a buffer, st_intersects asks do isochrone (i.e., buffer) polygons intersect park polygons?
st_intersects is followed by length, which tells R to sum the total number of TRUE (which mathematically are equal to one). This provides the sum total of intersected parks for each isochrone, which is then applied to the block group that the isochrone was built around.
vashon_iso <- vashon_iso |>
# `st_intersects` is Boolean: do buffer polygons intersect park polygons?
# for each row (`sapply`), asks whether a buffer intersects a park polygon
# `length` returns total number of elements for which this is TRUE
mutate(parks_in_buffer = sapply(st_intersects(buffer, parks),
length)) |>
# geometry no longer needed
st_drop_geometry()Finally, join these park access values (vashon_iso) with your block groups (vashon_bg).
vashon_bg <- left_join(vashon_bg, vashon_iso, by = "geoid")optional: click here to plot Vashon Island park access values
vashon_bg |>
ggplot() +
geom_sf(aes(fill = parks_in_buffer),
color = NA) +
scale_fill_gradientn(colors = c("#d95f02", "#e08214", "#a6d854", "#1b9e77"),
name = "Total\nParks",
breaks = c(1, 3, 5, 7)) +
theme_void() +
labs(title = "Vashon Island Block Group\nPark Access via Isochrones") +
theme(
plot.title = element_text(hjust = 0.5,
vjust = 0.2,
size = 15),
legend.title = element_text(hjust = 0.5,
vjust = 4,
size = 12),
legend.title.position = 'top',
legend.text = element_text(size = 10),
legend.text.position = 'right',
legend.position = 'right', # position relative to figure
legend.box.margin = margin(-50, 0, 0, -35),
legend.key.size = unit(0.1, 'inches')
)
Spatial Weights
In order to run a cluster analysis, a spatial weights matrix must be generated for the data. This matrix reveals which block groups are “neighbors” and how much influence (i.e., weight) each neighbor has on a block group.
Identifying Neighbors
The first step is to determine each block group’s neighbors. This will be done using the spdep package and poly2nb function. Neighbors are what allow you to determine whether clusters of similar park access values exist. In other words, cluster analysis answers the question, is a block group’s park access value similar to that of its neighbors?
Here, neighbors are defined using a Queen’s contiguity weights matrix, where block groups sharing any portion of a boundary are considered neighbors. For each block group, poly2nb returns a list of all neighboring block groups (i.e., all block groups that share a border).
neighbors <- poly2nb(vashon_bg$geometry)optional: click here to plot these neighbor connections
# polygon centroids
centroids <- st_centroid(vashon_bg)
# centroid coordinates (for cleaner connections on plot)
coords <- st_coordinates(centroids)
# df of edges for neighbors
edges_list <- lapply(seq_along(neighbors), function(i) {
nbrs <- neighbors[[i]]
if(length(nbrs) == 0) return(NULL)
# for each neighbor, create a pair of points (from i to neighbor)
data.frame(
from_x = coords[i,1],
from_y = coords[i,2],
to_x = coords[nbrs,1],
to_y = coords[nbrs,2]
)
})
edges <- do.call(rbind, edges_list)
# neighbor links map
vashon_bg |>
ggplot() +
# block groups
geom_sf(fill = NA,
color = "black",
linewidth = 0.5) +
# neighbor links
geom_segment(data = edges,
aes(x = from_x, y = from_y, xend = to_x, yend = to_y),
color = "red",
linewidth = 1) +
# centroids
geom_point(data = as.data.frame(coords),
aes(X, Y),
color = "black",
size = 2) +
theme_void() +
labs(title = "Neighbor Links for\nVashon Island Block Groups") +
theme(plot.title = element_text(face = "bold",
hjust = 0.5))
The red lines connecting each block group reflect the number of neighbors that block group has. For instance, the northern most block group shares a border with two other block groups, and thus has two red lines representing these connections.
Assigning Weights
Next, a numerical “weight” must be ascribed to every neighbor that a block group has. Weights reflect the influence each neighbor has over the park access value of the focal block group. In the context of cluster analysis, weights are necessary to compute the spatially weighted average of neighboring park access values for each block group.
To calculate these weights, supply the neighbors list (generated above) to the nb2listw function.
weights_matrix <- nb2listw(neighbors, style = "W", zero.policy = TRUE)Here, weights are row-standardized (style = "W"), meaning that the weights assigned to each block group’s neighbors will sum to one. For example, a block group with two neighbors assigns a weight of 0.5 to each, whereas a block group with four neighbors assigns 0.25 to each (see below).
In the above plots, the focal block group is in green and its neighbors are shades of red, with hues reflecting their weights. The two important takeaways from this figure are that all neighbors of a single block group carry the same influence, and block groups that are not neighbors carry no weight (shown in white).
optional: click here to view code for weights plot above
n_bg <- length(neighbors) # number of block groups
# convert weights to long form
weights_long <- lapply(1:n_bg, function(i) {
w <- rep(0, n_bg)
if(length(neighbors[[i]]) > 0) {
w[neighbors[[i]]] <- weights_matrix$weights[[i]] # assign neighbor weights
}
data.frame(focal_bg = i, target_bg = 1:n_bg, weight = w)
}) |>
bind_rows()
# join with spatial data
weights_long_sf <- weights_long |>
left_join(vashon_bg |>
mutate(target_bg = 1:nrow(vashon_bg)), by = "target_bg") |>
st_as_sf() |>
mutate(fill_weight = case_when(
focal_bg == target_bg ~ NA_real_, # focal bg (will be colored separately)
weight == 0 ~ NA_real_, # non-neighbors (white)
TRUE ~ weight # neighbors (gradient)
),
is_focal = focal_bg == target_bg)
# small multiples plot
weights_long_sf |>
ggplot() +
geom_sf(aes(fill = fill_weight),
color = "black") + # neighbors colored by weight
geom_sf(data = subset(weights_long_sf,
is_focal),
fill = "darkgreen",
color = "black") +
scale_fill_gradient(low = "lavenderblush",
high = "red",
na.value = "white",
name = "Weight") +
facet_wrap(~focal_bg, ncol = 4) +
theme_void() +
labs(title = "Neighor Weights for Vashon Island Block Groups",
subtitle = "focal block groups in green") +
theme(plot.title = element_text(face = "bold",
hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
strip.text = element_blank(),
strip.background = element_blank(),
legend.title = element_text(vjust = 4),
legend.box.margin = margin(0, 0, -130, -65),
legend.text = element_text(size = 10),
legend.key.size = unit(0.2, 'inches'))Analysis
Now return to the full King County sample to run the spatial cluster analysis. While Vashon Island provides a convenient and computationally simple way to demonstrate the workflow, it contains too few observations (n = 7) to run a cluster analysis.
Preparing the full sample for analysis follows the preceding steps and code identically. The only change is that generating isochrones for 1,544 block groups requires 4-6 hours.
optional: click here to view code for full sample preparation
# generate centroids
kc_centroids <- kc_bg |>
mutate(geometry = st_centroid(geometry))
# transform CRS to WGS84 (necessary for `osrmIsochrone` function)
temp <- st_transform(kc_centroids, "EPSG:4326")
empty_list <- list()
for (i in temp$geoid) {
tryCatch({
iso <- osrmIsochrone(loc = temp |> filter(geoid == i),
breaks = 30,
n = 500,
osrm.profile = "foot") |>
mutate(geoid = i) |>
select(geoid, geometry)
# add to list
empty_list[[as.character(i)]] <- iso
}, error = function(e) {
message(sprintf("Error at geoid %s: %s", i, e$message))
# add NULL to keep index consistent
empty_list[[as.character(i)]] <- NULL
})
}
# bind together
kc_iso <- bind_rows(empty_list)
# return spatial feature to desired CRS
kc_iso <- kc_iso |>
st_as_sf() |>
st_transform("EPSG:2926")
# park access counts
buffer <- st_buffer(kc_iso, dist = 0)
kc_iso <- kc_iso |>
# `st_intersects` is Boolean: do buffer polygons intersect park polygons?
# for each row (`sapply`), asks whether a buffer intersects a park polygon
# `length` returns total number of elements for which this is TRUE
mutate(parks_in_buffer = sapply(st_intersects(buffer, parks),
length)) |>
# drop geometry (no longer needed)
st_drop_geometry()
# join to King County block groups
kc_bg <- left_join(kc_bg, kc_iso, by = "geoid")
# identify neighbors
neighbors <- poly2nb(kc_bg$geometry)
# build spatial weights matrix
weights_matrix <- nb2listw(neighbors, style = "W", zero.policy = TRUE)Cluster Analysis
The aim of the cluster analysis is to identify areas in King County where park access is consistently high or low relative to neighboring block groups. In other words, which neighborhoods are surrounded by similarly high- or low-access areas, and where do access disparities appear concentrated?
To begin, Global Moran’s I provides a full sample assessment of spatial dependence in park access values. This statistic is generated using the moran.test function and indicates whether block groups with similar access levels tend to cluster together, spread apart, or are randomly distributed.
moran.test takes as arguments your park access values (parks_in_buffer column of kc_bg) and the spatial weights matrix (weights_matrix). Resultant values range from -1 to 1: values near -1 indicate strong negative spatial dependence (high-access areas near low-access areas), values near 1 indicate strong positive dependence (similar-access areas clustered together), and values near 0 suggest little or no spatial pattern.
# Global Moran's I
g_moran <- moran.test(kc_bg$parks_in_buffer,
listw = weights_matrix,
zero.policy = TRUE)Next, calculate Local Moran’s I using the localmoran function. Local Moran’s I evaluates the degree to which each block group’s park access value correlates with the weighted average park access of its neighbors. Positive values indicate that a block group has access levels similar to surrounding block groups (e.g., high access surrounded by high access or low access surrounded by low access), while negative values indicate spatial dissimilarity (e.g., high access isolated among low access or low access isolated among high access).
Like moran.test, localmoran takes as arguments your park access values (parks_in_buffer column of kc_bg) and the spatial weights matrix (weights_matrix).
# Local Moran's I
local_moran <- localmoran(kc_bg$parks_in_buffer,
listw = weights_matrix,
zero.policy = TRUE)
# add Local Moran's I values to kc_bg
kc_bg$local_moran <- local_moran[, "Ii"]Extract the p-value of each Local Moran’s I to distinguish statistically significant clusters (p < 0.05) when visualizing.
# add p-value to kc_bg
kc_bg$local_moran_p <- local_moran[, "Pr(z != E(Ii))"]Before classifying Local Moran’s I results into cluster types, the park access values are standardized into z-scores. Standardization centers the data around a mean of zero and scales it by the standard deviation, allowing each block group’s park access to be interpreted relative to the county-wide average rather than in terms of raw park counts. This makes it straightforward to distinguish observations with above-average (positive z-scores) and below-average (negative z-scores) park access.
# standardized park access
kc_bg$std_parks <- (
(kc_bg$parks_in_buffer - mean(kc_bg$parks_in_buffer, na.rm = TRUE)) /
sd(kc_bg$parks_in_buffer, na.rm = TRUE)
)Using these standardized values, calculate a spatial lag for each block group. The term “spatial lag” is adapted from time-series statistics, where a “lag” refers to a variable observed at an earlier point in time. In spatial analysis, the concept is extended to represent the value of a variable in neighboring locations rather than at earlier times. Thus, in the present context, a spatial lag represents the average park access value of a block group’s neighbors.
The lag.listw function takes the spatial weights matrix (weights_matrix) and the (standardized) park access values (std_parks column of kc_bg) as inputs. For each observation (i.e., block group), lag.listw computes the weighted average of the standardized park access values for that block group’s neighboring block groups, where the weights are defined by the spatial weights matrix.
# spatial lag of standardized park count
kc_bg$lag_parks <- lag.listw(weights_matrix, kc_bg$std_parks, zero.policy = TRUE)Cluster types are then assigned by comparing each block group’s standardized park access with the average standardized park access of its neighboring block groups. Statistically significant observations (p < 0.05) are classified as High-High, Low-Low, High-Low, or Low-High, while all others are labeled Not Significant. High-High clusters indicate areas of consistently high park access, whereas Low-Low clusters indicate areas of consistently low park access.
# identify cluster types
kc_bg <- kc_bg |>
mutate(
cluster_type = case_when(
(std_parks >= 0) & (lag_parks >= 0) & (local_moran_p < 0.05) ~ "High-High",
(std_parks <= 0) & (lag_parks <= 0) & (local_moran_p < 0.05) ~ "Low-Low",
(std_parks >= 0) & (lag_parks <= 0) & (local_moran_p < 0.05) ~ "High-Low",
(std_parks <= 0) & (lag_parks >= 0) & (local_moran_p < 0.05) ~ "Low-High",
TRUE ~ "Not Significant"
),
cluster_type = factor(cluster_type, levels = c("High-High",
"Low-Low",
"High-Low",
"Low-High",
"Not Significant"))
)
The figure above plots park access clusters based on a 30-minute isochrone interval. Global Moran’s I is nearly equal to 1 (0.94) and highly statistically significant (p < 0.001), indicating that block groups with high access tend to be found clustered together, while block groups with low access are also found near to one another.
This conclusion is corroborated by the Local Moran’s I values, which are displayed in the maps themselves. Specifically, high-access areas (colored in teal), represent regions where the park access of a block group is relatively high and the park access of neighboring block groups is relatively high. High-access clustering appears primarily within Seattle block groups, in Northwestern King County. Conversely, orange areas reflect the opposite – low-access regions where a block group has relatively low park access and is surrounded by block groups with similarly low park access. Low-access clustering is found in the eastern and southern regions of King County.
optional: click here to view code for cluster analysis plot
# define a p-value function (for map text)
p_value_text <- function(x) { # x = p-value
ifelse(x < 0.001, 'p < 0.001',
ifelse(x < 0.01, 'p < 0.01',
ifelse(x < 0.05, 'p < 0.05', 'insignificant')
)
)
}
kc_bg |>
ggplot() +
geom_sf(aes(fill = cluster_type), color = NA, size = 0.001) +
scale_fill_manual(
# plot colors
values = c("High-High" = "#1b9e77",
"Low-Low" = "#d95f02",
"High-Low" = "#7570b3",
"Low-High" = "#e7298a70",
"Not Significant" = "lightgrey"),
# legend labels
labels = c("High-High" = "high access",
"Low-Low" = "low access",
"High-Low" = "high-low outlier",
"Low-High" = "low-high outlier",
"Not Significant" = "insignificant")
) +
labs(title = "King County Park Access",
subtitle = "30 Minute Isochrone Interval",
caption = paste0("Global Moran's I: ",
round(g_moran$estimate[1], 2),
"\n",
p_value_text(g_moran$p.value)),
fill = "") +
theme_void() +
theme(plot.title = element_text(face = "bold",
hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0.5,
vjust = 25),
legend.box.margin = margin(0, 0, -10, 0),
legend.position = 'bottom',
legend.direction = 'horizontal',
legend.key.size = unit(0.3, 'inches'))Bivariate Association
You can expand this analysis and incorporate an equity lens by considering how park access patterns relate spatially to residential sociodemographic characteristics. For instance, are areas with comparatively low park access located near areas with higher concentrations of children? To examine the spatial relationship between park access and youth population proportions, conduct a bivariate spatial association analysis.
Here, youth proportions are calculated as the proportion of residents under 18 years of age in each block group. These data are obtained from the 2023 U.S. Census Bureau’s American Community Survey population estimates.
Note
A Census API key may be required to access Census data. If you do not have one, request here (free). If you already have a key, install with: tidycensus::census_api_key("YOUR_KEY_HERE", install = TRUE)
# youth population proportions
kc_pop <- get_acs(geography = "block group",
state = "WA",
county = "King",
year = 2023,
survey = "acs5",
geometry = FALSE,
variables = c(
# total population
"B01001_001",
# under 18 male
"B01001_003", "B01001_004", "B01001_005", "B01001_006",
# under 18 female
"B01001_027", "B01001_028", "B01001_029", "B01001_030"
)) |>
clean_names() |>
group_by(geoid) |>
summarize(total_pop = sum(estimate[variable == "B01001_001"]),
under_18 = sum(estimate[variable != "B01001_001"]),
.groups = "drop") |>
# youth rate (accounts for population dense areas)
mutate(youth_rate = case_when(total_pop == 0 ~ 0,
TRUE ~ under_18 / total_pop))Next, merge these youth proportions into the park access data and standardize.
kc_bg <- kc_bg |>
left_join(kc_pop, by = "geoid")
# standardize youth population
kc_bg$std_youth <- (
(kc_bg$youth_rate - mean(kc_bg$youth_rate, na.rm = TRUE)) /
sd(kc_bg$youth_rate, na.rm = TRUE)
)Using these standardized values and the lag.listw function, calculate a spatial lag for each block group (same as above).
# spatial lag of standardized youth proportion
kc_bg$lag_youth <- lag.listw(weights_matrix, kc_bg$std_youth, zero.policy = TRUE)Using the moran.bi function, perform a bivariate analysis between park access (std_parks) and youth proportions (std_youth). moran.bi evaluates the relationship between park access in each block group and youth concentrations in surrounding block groups.
# bivariate analysis
bv <- moran.bi(varX = kc_bg$std_parks,
varY = kc_bg$std_youth,
listw = weights_matrix)Finally, classify the into four quadrants based on standardized park access and the spatial lag of youth proportion. The first term refers to the focal block group’s park access, while the second term refers to the youth proportion of neighboring block groups: High-High (above-average park access and above-average neighboring youth proportion), Low-Low (below-average park access and below-average neighboring youth proportion), High-Low (above-average park access and below-average neighboring youth proportion), and Low-High (below-average park access and above-average neighboring youth proportion).
# cluster types
kc_bg <- kc_bg |>
mutate(cluster_type = case_when(std_parks >= 0 & lag_youth >= 0 ~ "High Parks - High Youth",
std_parks <= 0 & lag_youth <= 0 ~ "Low Parks - Low Youth",
std_parks >= 0 & lag_youth <= 0 ~ "High Parks - Low Youth",
std_parks <= 0 & lag_youth >= 0 ~ "Low Parks - High Youth"),
cluster_type = factor(cluster_type,
levels = c("High Parks - High Youth",
"Low Parks - Low Youth",
"High Parks - Low Youth",
"Low Parks - High Youth")))
optional: click here to view code for bivariate association plot
kc_bg |>
ggplot() +
geom_sf(aes(fill = cluster_type), color = NA, size = 0.001) +
scale_fill_manual(
# plot colors
values = c("High Parks - High Youth" = "#1b9e77",
"Low Parks - Low Youth" = "#d95f02",
"High Parks - Low Youth" = "#7570b3",
"Low Parks - High Youth" = "#e7298a",
"Not Significant" = "lightgrey"),
# legend labels
labels = c("High Parks - High Youth" = "High Parks &\nHigh Youth",
"Low Parks - Low Youth" = "Low Parks &\nLow Youth",
"High Parks - Low Youth" = "High Parks & Low Youth\n(potential oversupply)",
"Low Parks - High Youth" = "Low Parks & High Youth\n(potential access gap)")
) +
labs(title = "Bivariate Association: Park Access & Youth Share",
subtitle = "30 Minute Isochrone Interval",
caption = paste0("Global Bivariate Moran's I: ",
round(bv$I, 2)),
fill = element_blank()) +
theme_void() +
theme(plot.title = element_text(face = "bold",
hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5),
plot.caption = element_text(hjust = 0.5,
vjust = 25),
legend.position = 'bottom',
legend.direction = 'horizontal',
legend.key.size = unit(0.3, 'inches'))The analysis reveals evidence of a significant bivariate spatial association across King County. The dominant clustering pattern is regions with low park access and higher-than-average neighboring youth proportions, suggesting areas where limited park accessibility overlaps with nearby concentrations of youth populations. Conversely, regions of high park access with lower-than-average neighboring youth proportions are also present, primarily in the Seattle area. These clusters indicate areas where park accessibility may exceed the relative concentration of youth populations.
This analysis represents one of many possible extensions of the spatial framework presented in this study. By comparing the spatial distribution of park access with a sociodemographic variable such as youth proportion, the analysis helps identify areas where spatial mismatches between park accessibility and community characteristics may be most pronounced and where park provision could be most beneficial. Similar bivariate spatial approaches could be extended to other vulnerable or policy-relevant populations, including older adults, low-income households, communities of color, and households without access to private transportation. In addition, future work could incorporate health outcomes, such as physical activity levels or chronic disease prevalence, or examine access to other types of amenities (e.g., schools, transit, or healthcare facilities) to further contextualize patterns of spatial equity across King County.