Senior Accessibility in Hamilton, Ontario

Introduction

This analysis aims to provide a simple example of how OpenTripPlanner can be used to calculate accessibility scores for a target population given a target municipal service. In this case, we will calculate accessibility scores based on the number of free public municipal facilities (libraries, community centres, recreation centres, and service centres) located within a 15-minute commute of the current senior population (age 65+) in Hamilton, Ontario. This analysis is based on the principles of aging-in-place and the 15-minute city.

This RPub is organized into five parts as follows:

• PART 1 - Setup Environment
• PART 2 - Load and Prepare Data
• PART 3a - Setup OpenTripPlanner
• PART 3b - Connect and Test OpenTripPlanner
• PART 3c - Query OpenTripPlanner for Isochrones
• PART 4 - Calculate Accessibility Scores
• PART 5 - Plot Results

To complete this analysis, we will be using the following data:

• Census Division and Dissemination Area shapefiles from the 2021 Canadian Census downloaded from Statistics Canada here.
• Census of Population data at the Dissemination Area level downloaded from Statistics Canada here.
• Library, Community Centre, Recreation Centre, and Service Centre shapefiles downloaded from Open Hamilton here.

NOTE: For a great guide on accessibility analysis in R using OpenTripPlanner see Hussein Mahfouz’s excellent RPub here.

PART 1 - Setup Environment

1.1 - Load required library packages.

# for rpubs
library(rmarkdown)
library(knitr)

# for input / outputs
library(here)

# for data manipulation
library(tidyverse)
library(tidyselect)
library(magrittr)
library(janitor)

# for spatial data manipulation
library(sf)

# for OpenTripPlanner
library(osmextract)
library(opentripplanner)
library(RcppSimdJson)

# for plotting
library(ggsn)
library(ggspatial)
library(tmap)
library(tmaptools)
library(cowplot)

1.2 - Set preferences.

# set working directory for our inputs/outputs (optional)
#setwd()

# OTP requires coordinates to be in CRS 4326 (WGS 84) so we will set a common crs
common_crs = 4326

# set tmap view mode to static or interactive ("plot" or "view")
tmap_mode("plot")

tmap mode set to plotting

1.3 - Set functions.

Thank you again to Hussein Mahfouz’s RPub for this handy code.

# OTP requires latitude and longitude coordinates as numeric values in two separate columns
# we will use the function below split the geometry column of spatial data needed for OTP

sfc_as_cols <- function(x, names = c("lon","lat")) {
  stopifnot(inherits(x,"sf") && inherits(sf::st_geometry(x),"sfc_POINT"))
  ret <- sf::st_coordinates(x)
  ret <- tibble::as_tibble(ret)
  stopifnot(length(names) == ncol(ret))
  x <- x[ , !names(x) %in% names]
  ret <- setNames(ret,names)
  dplyr::bind_cols(x,ret)
}

PART 2 - Load and Prepare Data

To make this example as easy as possible, data has already been downloaded and pre-cleaned from the links above as follows:

• cd21_ham = Census Division shapefile reduced to City of Hamilton boundaries and columns not needed for this analysis removed.
• da21_ham = Dissemination Area shapefile reduced to City of Hamilton boundaries and columns not needed for this analysis removed.
• da21_data = Census of Population data by Dissemination Area reduced to City of Hamilton boundaries and the senior population aged 65+ in 2021, 2026, and 2031.
• cspaces = Library, Community Centre, Recreation Centre, and Service Centre shapefiles merged into a single Community Spaces shapefile.

2.1 - Load cleaned data.

# Census Division shapefile
cd21_ham <- sf::st_read(here::here("data/clean/cd21_ham.gpkg"),
                        stringsAsFactors = FALSE)
# Dissemination Area shapefile
da21_ham <- sf::st_read(here::here("data/clean/da21_ham.gpkg"),
                        stringsAsFactors = FALSE)
# Community Spaces shapefile
cspaces <- sf::st_read(here::here("data/clean/cspaces.gpkg"),
                       stringsAsFactors = FALSE)
# Dissemination Area data
da21_data <- read_csv(here::here("data/clean/da21_ham_srpop.csv"),
                      locale = locale(encoding = "latin1"),
                      na = "n/a")

Let’s have a look at our data.

# Census Division shapefile
head(cd21_ham)

Simple feature collection with 1 feature and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -80.24849 ymin: 43.05052 xmax: -79.62221 ymax: 43.47106
Geodetic CRS:  WGS 84
  CDUID         DGUID   CDNAME CDTYPE LANDAREA PRUID
1  3525 2021A00033525 Hamilton    CDR  1118.31    35
                            geom
1 MULTIPOLYGON (((-79.79927 4...

# Dissemination Area shapefile
head(da21_ham)

Simple feature collection with 6 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -79.94766 ymin: 43.26654 xmax: -79.88898 ymax: 43.30781
Geodetic CRS:  WGS 84
     DAUID             DGUID                           geom
1 35250030 2021S051235250030 MULTIPOLYGON (((-79.92591 4...
2 35250031 2021S051235250031 MULTIPOLYGON (((-79.89038 4...
3 35250032 2021S051235250032 MULTIPOLYGON (((-79.93015 4...
4 35250033 2021S051235250033 MULTIPOLYGON (((-79.93278 4...
5 35250034 2021S051235250034 MULTIPOLYGON (((-79.94082 4...
6 35250035 2021S051235250035 MULTIPOLYGON (((-79.94382 4...

# Community Spaces shapefile
head(cspaces)

Simple feature collection with 6 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -79.98109 ymin: 43.1217 xmax: -79.80343 ymax: 43.39605
Geodetic CRS:  WGS 84
                NAME                       geom
1   Ancaster Library POINT (-79.97666 43.22531)
2     Barton Library POINT (-79.84122 43.25803)
3   Binbrook Library  POINT (-79.80343 43.1217)
4   Carlisle Library POINT (-79.98109 43.39605)
5 Concession Library POINT (-79.85137 43.24104)
6     Dundas Library POINT (-79.95495 43.26548)

# Dissemination Area data
head(da21_data)

2.2 - Dissemination Area Data

We will start preparing data by merging the Dissemination Area shapefile with the senior population data.

# merge DA data with DA shapefile on DGUID
da21_ham <- merge(x = da21_ham,
                  y = da21_ham_srpop,
                  by.x = "DGUID",
                  by.y = "DGUID")

# final data check to confirm merge
str(da21_ham)

Classes ‘sf’ and 'data.frame':  891 obs. of  6 variables:
 $ DGUID   : chr  "2021S051235250030" "2021S051235250031" "2021S051235250032" "2021S051235250033" ...
 $ DAUID   : chr  "35250030" "35250031" "35250032" "35250033" ...
 $ srpop21 : num  90 120 115 160 105 75 95 140 320 215 ...
 $ srpop26 : num  120 145 180 205 155 105 145 185 370 240 ...
 $ srpop31 : num  165 190 230 250 210 120 185 225 435 280 ...
 $ geometry:sfc_MULTIPOLYGON of length 891; first list element: List of 1
  ..$ :List of 1
  .. ..$ : num [1:188, 1:2] -79.9 -79.9 -79.9 -79.9 -79.9 ...
  ..- attr(*, "class")= chr [1:3] "XY" "MULTIPOLYGON" "sfg"
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA
  ..- attr(*, "names")= chr [1:5] "DGUID" "DAUID" "srpop21" "srpop26" ...

# we will also make some statistical variables for later
srpop21total = sum(da21_ham_srpop$srpop21)
srpop26total = sum(da21_ham_srpop$srpop26)
srpop31total = sum(da21_ham_srpop$srpop31)

print(paste0("Hamilton's total senior population in 2021: ", srpop21total))

[1] "Hamilton's total senior population in 2021: 104260"

print(paste0("Hamilton's total senior population in 2026: ", srpop26total))

[1] "Hamilton's total senior population in 2026: 142190"

print(paste0("Hamilton's total senior population in 2031: ", srpop31total))

[1] "Hamilton's total senior population in 2031: 183095"

2.3 - Dissemination Area Centroids

Now we will create centroids from Dissemination Areas that we will use to query with OTP.

# create centroids of Hamilton's DA boundaries
# another name for centroid is 'geometric center' so we will use (gc) here
da21_ham_gc <- st_centroid(da21_ham)

# recall OTP requires latitude and longitude coordinates as numeric values in two separate columns
# we will use the function from above to split geometry and add results as new columns
da21_ham_gc <- sfc_as_cols(da21_ham_gc)

# confirm split with a final data check
#summary(da21_ham_gc)
str(da21_ham_gc)

Classes ‘sf’ and 'data.frame':  891 obs. of  8 variables:
 $ DGUID   : chr  "2021S051235250030" "2021S051235250031" "2021S051235250032" "2021S051235250033" ...
 $ DAUID   : chr  "35250030" "35250031" "35250032" "35250033" ...
 $ srpop21 : num  90 120 115 160 105 75 95 140 320 215 ...
 $ srpop26 : num  120 145 180 205 155 105 145 185 370 240 ...
 $ srpop31 : num  165 190 230 250 210 120 185 225 435 280 ...
 $ lon     : num  -79.9 -79.9 -79.9 -79.9 -79.9 ...
 $ lat     : num  43.3 43.3 43.3 43.3 43.3 ...
 $ geometry:sfc_POINT of length 891; first list element:  'XY' num  -79.9 43.3
 - attr(*, "sf_column")= chr "geometry"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA
  ..- attr(*, "names")= chr [1:7] "DGUID" "DAUID" "srpop21" "srpop26" ...

2.4 - Community Spaces

Let’s have another look at our Community Spaces data.

# final data check
str(cspaces)

Classes ‘sf’ and 'data.frame':  97 obs. of  2 variables:
 $ NAME: chr  "Ancaster Library" "Barton Library" "Binbrook Library" "Carlisle Library" ...
 $ geom:sfc_POINT of length 97; first list element:  'XY' num  -80 43.2
 - attr(*, "sf_column")= chr "geom"
 - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
  ..- attr(*, "names")= chr "NAME"

# there are 97 cspaces, let's make a variable to store this
cspacestotal <- 97

2.5 - Plot

Finally, let’s take a look at the study area together with all our data.

# plot together
tm_shape(da21_ham) +
  tm_borders(col = "black", lwd = 0.1) +
tm_shape(da21_ham_gc) +
  tm_bubbles(col = "blue", scale = .3) +
tm_shape(cspaces) +
  tm_bubbles(col = "red", scale = .5)

Perfect!

PART 3a - Setup OpenTripPlanner (OTP)

Sections 3a, 3b, and 3c are based on the following OTP package vignettes:

• Prerequisites
• Getting Started
• Advanced Features

NOTE: Java 8 is required to run OTP. Download Java 8 here.

3.1 Create the required OTP folder and data structure.

As outlined in the vignettes above, OTP requires a certain folder structure to run. An example folder structure is provided below for reference. We will recreate this structure step-by-step with the exception of Steps 6 and 8 which we will skip as they are not required for this analysis.

Step	Example	Description
	/wd	Working directory
1	/otp	Top folder for storing all OTP data
2	/graphs	Subfolder containing all graphs
3	/default	Subfolder that acts as OTP router
4	osm.pbf	Required OpenStreetMap road map
5	router-config.json	Optional config file
6	build-config.json	Optional config file
7	gtfs.zip	Optional GTFS data
8	dem.tif	Optional Elevation data

Step 1: Create the OTP folder.

# confirm correct version of Java (optional)
#otp_check_java()

# Step 1 - create otp folder
path_data <- file.path(here::here(),
                       "otp")
dir.create(path_data)

# download the required OTP jar file
path_otp <- otp_dl_jar(path_data, cache = FALSE)

Step 2: Create the graphs folder.

For some reason I could not sort out, OTP required the “graphs” folder to be created on my computer first through the demo function prior to being able to create my own graph or router. This function also downloads the “default” folder automatically

otp_dl_demo(path_data)

Step 3: Create the router folder.

dir.create(file.path(here::here("otp/graphs/hamilton")))

Step 4: Download OpenStreetMap (OSM) data.

This step uses the osmextract package and is based on the Introducing osmextract Vignette.

NOTE: Downloading OSM data can take a few minutes as the files can be quite large.

# geocode hamilton's location
hamilton = tmaptools::geocode_OSM("Hamilton, Canada")$coords

# Check osm data for hamilton's location
oe_match(hamilton, provider = "geofabrik")

The input place was matched with Ontario. 

$url
[1] "https://download.geofabrik.de/north-america/canada/ontario-latest.osm.pbf"

$file_size
[1] 762458320

oe_match(hamilton, provider = "openstreetmap_fr")

The input place was matched with Golden Horseshoe. 

$url
[1] "http://download.openstreetmap.fr/extracts/north-america/canada/ontario/golden_horseshoe-latest.osm.pbf"

$file_size
[1] 129940802

Here we can see that geofrabrik shows data is only available for the province, whereas openstreetmap_fr has data available at the regional level. We will download the data from the regional level as it has a significantly smaller file size.

# download osm data for Hamilton
oe_get("Hamilton", 
       provider = "openstreetmap_fr",
       download_directory = here::here("otp/graphs/hamilton"),
       boundary = cd21_ham,
       boundary_type = c("spat", "clipsrc"),
       download_only = TRUE, 
       skip_vectortranslate = TRUE)

Step 5: Configure the router.

The router-config file dictates how OTP will calculate route options. For example, the router-config file specifies the default walking speed from Point A to Point B. Since our analysis is on the senior population, and it is reasonable to assume that average walking speed decreases with age, we will configure our router to reflect this. But by how much? Luckily researchers at McMaster University in Hamilton have already completed the leg work for us here.

# make a config object
router_config <- otp_make_config("router")     

# check default walking speed
router_config$routingDefaults$walkSpeed

[1] 1.34

The default walking speed of OTP is 1.34 metres per second. Let’s see what the McMaster researchers found.

# let's create a dataframe of McMaster University's findings on senior walking speed
walk_speed <- data.frame("age" = c("60-69", "70-79", "80-89"),
                         "srmale" = c(1.34, 1.26, 0.97),
                         "srfemale" = c(1.24, 1.13, 0.94))
walk_speed

# now let's calculate the average senior walking speed
srmale_avg <- (sum(walk_speed$srmale)/3)
srfemale_avg <- (sum(walk_speed$srfemale)/3)
sravg <- ((srmale_avg + srfemale_avg)/2)
sravg

[1] 1.146667

The average senior walking speed is 1.15 metres per second. Let’s recalibrate OTP’s default walking speed.

# calibrate the otp router
router_config$routingDefaults$walkSpeed <- 1.15 

# confirm change
router_config$routingDefaults$walkSpeed

[1] 1.15

# Save the config file
otp_write_config(router_config,
                 dir = path_data,
                 router = "hamilton")

Step 7: Download the GTFS feed.

GTFS stands for General Transit Feed Specification and, according to gtfs.org, “is a data specification that allows public transit agencies to publish their transit data in a format that can be consumed by a wide variety of software applications”. If you use Google Maps to plan your route using transit, you have benefited from a GTFS feed.

Hamilton’s transit agency is called the Hamilton Street Railway (HSR), which ironically has not operated any rail related transit since 1951. Tom Luton maintains a neat website called Hamilton Transit History which is worth a look if interested. The hyperlink for HSR’s GTFS feed can be found here.

NOTE: The GTFS feed must be in a zipped folder with “gtfs” in the name.

url <- "https://googlehsrdocs.hamilton.ca/"  
destfile = here::here("otp",
                      "graphs",
                      "hamilton",
                      "ham_gtfs.zip")
download.file(url, destfile, mode="wb")

trying URL 'https://googlehsrdocs.hamilton.ca/'
Content type 'application/x-zip-compressed' length 8507622 bytes (8.1 MB)
downloaded 8.1 MB

Our OTP router is now ready!

PART 3b - Connect and Test OTP

Now that our data is in order and the OTP is configured, we can test our connection. A few public service announcements you should be aware of before we move any further:

1. Building and querying OTP can be quite computationally intensive and several OTP functions have options to increase the amount memory or the number of cores your computer uses in order to speed up the process.

2. It is entirely up to your discretion, given the size of your data and your computers abilty, whether to increase the memory or cores. For reference, this entire script uses the OSM defaults and runs in around five minutes from start to finish.

3. Where applicable below, I have indicate how to increase memory and cores should you wish to do so.

NOTE: Once connected, OTP will automatically load an interactive web interface in your browser. You can play with the interface as you wish or ignore it entirely, it is not required for the analysis.

WARNING: Never set ‘ncores’ to more cores than are available on your computer. To check the number of cores on your computer, hold ctrl + shift + esc to bring up task manager. Click the performance tab and the number of cores should be indicated in the bottom right.

3.2 - Build the OTP graph.

# create the Hamilton OTP graph file
# optional memory bump up - see explanation above
log1 <- otp_build_graph(otp = path_otp, 
                        dir = path_data, 
                        memory = 2048,        #default
                        router = "hamilton")

# launch otp and load the hamilton graph (the server)
log2 <- otp_setup(otp = path_otp, 
                  dir = path_data,
                  router = "hamilton")

3.3 - Make the connection.

# connect to OTP from R
otpcon <- otp_connect(router = "hamilton")

Router http://localhost:8080/otp/routers/hamilton exists

NOTE: the OTP web interface will auto-load when ready at http://localhost:8080.

3.4 - Test the connection.

# test by getting a route from McMaster University to Tim Hortons Field
ham_test1 <- otp_plan(otpcon, 
                      fromPlace = c(-79.917054, 43.258032), 
                      toPlace = c(-79.830855, 43.251014))

  |                                                  | 0 % ~calculating  
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s  

# view the route
osm_ham <- read_osm(ham_test1, ext=1.1)
tm_shape(osm_ham) + tm_rgb() +
  tm_shape(ham_test1) +
  tm_lines(col = 'blue', lwd = 2, alpha = 0.8)

3.5 - Test the isochrone function.

# test by getting an isochrone for McMaster University
# optional ncores bump up - see explanation above
ham_test2  <- otp_isochrone(otpcon = otpcon,
                            fromPlace = c(-79.917151, 43.258050),
                            mode = c("WALK", "TRANSIT"),
                            date_time = as.POSIXct(strptime("2022-08-03 09:00", "%Y-%m-%d %H:%M")),
                            maxWalkDistance = 1000,
                            cutoffSec = (c(5, 10, 15) * 60),
                            ncores = 1)    #default

  |                                                  | 0 % ~calculating  
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=01s  

# convert to minutes for plotting
ham_test2$minutes = (ham_test2$time / 60)

# view the isochrone
osm_ham <- read_osm(ham_test2, ext=1.1)
tm_shape(osm_ham) + tm_rgb() +
  tm_shape(ham_test2) +
  tm_fill("minutes",
          breaks = c(0, 5.01, 10.01, 15.01),
          style = "fixed",
          palette = "-RdYlBu") +
  tm_borders()

Brilliant, it works!

PART 3c - Query OpenTripPlanner for Isochrones

3.6 - Query OTP.

We will now query OTP to get 15-minute isochrones from all 891 Dissemination Area centroids.

NOTE: The function below uses the default ncores setting of ncores = 1. Therefore, expect the function to take a few minutes to compute. You can also increase ncores where indicated below if desired, however, you should read the warning at the beginning of this section first.

# get isochrone polygons for each DA centroid via loop

# variable with number of DAs
nrows <- nrow(da21_ham_gc)

# empty list to store output
da21_ham_gc_iso <- list()

# get isochrone polygon for each DA centroid and store in da21_ham_iso list
# cutoffsec is the max time in seconds for the isochrone (15 minutes for this analysis)
# 
for (i in 1:nrows){
  da21_ham_gc_iso[[i]] <- otp_isochrone(otpcon = otpcon,
                                        fromPlace = c(da21_ham_gc$lon[i], da21_ham_gc$lat[i]),
                                        fromID = c(da21_ham_gc$DAUID[i]),
                                        mode = c("WALK", "TRANSIT"),
                                        date_time = as.POSIXct(strptime("2022-08-03 09:00", "%Y-%m-%d %H:%M")),
                                        maxWalkDistance = 1000,
                                        cutoffSec = (15 * 60),
                                        ncores = 1)}  #default

# check first record to confirm complete
head(da21_ham_gc_iso, 1)

[[1]]
Simple feature collection with 1 feature and 3 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -79.9254 ymin: 43.294 xmax: -79.9117 ymax: 43.3058
Geodetic CRS:  WGS 84
  id time fromPlace                       geometry
1  1  900  35250030 MULTIPOLYGON (((-79.9236 43...

3.7 - Clean the results.

OTP is typically not able to query from all centroids. There are ways to address this, but to keep things simple we will not go into that level of detail here.

# these would show as NA values - lets check
sum(is.na(da21_ham_gc_iso))

[1] 32

32 centroids from our 891 dissemination areas have zero values. We will remove these to facilitate cleaning.

# remove NA values
da21_ham_gc_iso <- da21_ham_gc_iso[da21_ham_gc_iso != "NA"]

# convert result to a single simple feature
da21_ham_gc_iso = do.call(rbind, da21_ham_gc_iso)

# select desired columns 
da21_ham_gc_iso <- da21_ham_gc_iso %>%
  dplyr::select(c(3, 4))

# rename columns to facilitate merger to make data whole again
da21_ham_gc_iso <- da21_ham_gc_iso %>%
  dplyr::rename("DAUID"= 1)

# check
head(da21_ham_gc_iso, 3)

Simple feature collection with 3 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -79.941 ymin: 43.2761 xmax: -79.9056 ymax: 43.3058
Geodetic CRS:  WGS 84
     DAUID                       geometry
1 35250030 MULTIPOLYGON (((-79.9236 43...
2 35250031 MULTIPOLYGON (((-79.9198 43...
3 35250032 MULTIPOLYGON (((-79.9399 43...

3.8 - Practice plots.

Plot a single isochrone to see what it looks like.

# test plot a single isochrone
osm_ham <- read_osm(cd21_ham, ext=1.1)
tm_shape(osm_ham) + tm_rgb() +
  tm_shape(da21_ham_gc_iso[da21_ham_gc_iso$DAUID == "35250468", ]) +
  tm_fill(col = "#E5854B") +
  tm_borders()

Now let’s try plotting all the isochrones at once.

# test plot all isochrones
osm_ham <- read_osm(cd21_ham, ext=1.1)
tm_shape(osm_ham) + tm_rgb() +
  tm_shape(da21_ham_gc_iso) +
  tm_fill(col = "#E5854B") +
  tm_borders()

3.9 - Disconnect from OTP.

Make sure to disconnect from the OTP server.

NOTE: OTP will keep running in the background of your computer until you disconnect.

# otp no longer needed - stop connection to otp router
otp_stop(warn=FALSE)

[1] "SUCCESS: The process \"java.exe\" with PID 16532 has been terminated."

PART 4 - Calculate Accessibility Scores

We will calculate accessibility scores as follows:

Score	Equation	Indicator
Accessibility Score	Number of community spaces reached per dissemination area divided by the total number of community spaces in the city	Percentage of total community spaces in the city reachable per dissemination area
Average Accessibility Score	Sum of accessibility scores for all dissemination areas divided by the total number dissemination areas in the city	Average percentage of total community spaces in the city reachable per dissemination area
Weighted Accessibility Score	Number of community spaces reachable per dissemination area multiplied by its senior population and then divided by the total number of community spaces in the city	Percentage of total community spaces reachable per senior per dissemination area
Average Weighted Accessibility Score	Sum of weighted accessibility scores for all dissemination areas divided by the total number dissemination areas in the city	Average percentage of total community spaces reachable per senior per dissemination area

4.1 - Calculate Community Spaces Reached.

Determine the number of community spaces reached for each dissemination area isochrone and add it as a column to our data.

# perform intersection of 15 minute isochrones and cspaces, then add point count as numeric to each polygon
da21_ham_gc_iso_int <- da21_ham_gc_iso %>%
  mutate(cspaces_15 = (as.numeric(lengths(st_intersects(da21_ham_gc_iso, cspaces)))))

# confirm
da21_ham_gc_iso_int

Simple feature collection with 859 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -80.2156 ymin: 43.0702 xmax: -79.6186 ymax: 43.4248
Geodetic CRS:  WGS 84
First 10 features:
      DAUID                       geometry cspaces_15
1  35250030 MULTIPOLYGON (((-79.9236 43...          0
2  35250031 MULTIPOLYGON (((-79.9198 43...          0
3  35250032 MULTIPOLYGON (((-79.9399 43...          0
4  35250033 MULTIPOLYGON (((-79.9439 43...          0
5  35250034 MULTIPOLYGON (((-79.9495 43...          0
6  35250035 MULTIPOLYGON (((-79.9448 43...          2
7  35250036 MULTIPOLYGON (((-79.9598 43...          0
8  35250037 MULTIPOLYGON (((-79.9634 43...          2
9  35250038 MULTIPOLYGON (((-79.9629 43...          3
10 35250039 MULTIPOLYGON (((-79.9601 43...          2

# remove geometry first to facilitate merge
st_geometry(da21_ham_gc_iso_int) <- NULL

# confirm
da21_ham_gc_iso_int

# Merge back into the DA shape file on DGUID to plot results
da21_ham <- left_join(x = da21_ham,
                      y = da21_ham_gc_iso_int,
                      by.x = "DGUID",
                      by.y = "DGUID")

Joining, by = "DAUID"

# confirm
da21_ham

Simple feature collection with 891 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -80.24849 ymin: 43.05052 xmax: -79.62221 ymax: 43.47106
Geodetic CRS:  WGS 84
First 10 features:
               DGUID    DAUID srpop21 srpop26 srpop31 cspaces_15
1  2021S051235250030 35250030      90     120     165          0
2  2021S051235250031 35250031     120     145     190          0
3  2021S051235250032 35250032     115     180     230          0
4  2021S051235250033 35250033     160     205     250          0
5  2021S051235250034 35250034     105     155     210          0
6  2021S051235250035 35250035      75     105     120          2
7  2021S051235250036 35250036      95     145     185          0
8  2021S051235250037 35250037     140     185     225          2
9  2021S051235250038 35250038     320     370     435          3
10 2021S051235250039 35250039     215     240     280          2
                         geometry
1  MULTIPOLYGON (((-79.92591 4...
2  MULTIPOLYGON (((-79.89038 4...
3  MULTIPOLYGON (((-79.93015 4...
4  MULTIPOLYGON (((-79.93278 4...
5  MULTIPOLYGON (((-79.94082 4...
6  MULTIPOLYGON (((-79.94382 4...
7  MULTIPOLYGON (((-79.94652 4...
8  MULTIPOLYGON (((-79.94968 4...
9  MULTIPOLYGON (((-79.95334 4...
10 MULTIPOLYGON (((-79.94972 4...

4.2 - Calculate Accessibility Scores.

Calculate accessibility score of each dissemination area by dividing the results from STEP 1 by the total # of community spaces in Hamilton.

# calcualte score (% of the total community spaces in Hamilton that are reachable from each DA) as new column
da21_ham <- da21_ham %>%
  mutate(DA_SSOS = ((da21_ham$cspaces_15 / cspacestotal) * 100))

# recall that we dropped all cases where OTP could not route from a DA centroid
# let's have a look
sum(is.na(da21_ham))

[1] 64

64 NA values. Let’s check how many distinct values are in the DAUID column to make sure we have 891 - one for each dissemination area.

# check for 891 unique values
n_distinct(da21_ham$DAUID)

[1] 891

891 values. Perfect, our data is whole again.

NOTE: The key assumption here is that if OTP can not route from a location, then there are no community spaces reachable from said location.

# now let's perform a master replace of "NA" with "0" for all other columns
da21_ham[is.na(da21_ham)] <- 0

# confirm
sum(is.na(da21_ham))

[1] 0

No more NA values!

4.3 - Calculate Weighted Accessibility Scores.

Calculate weighted accessibility scores of each dissemination area by multiplying the results from STEP 2 by the senior population of each dissemination area and then dividing by the total senior population in Hamilton.

# use mutate to calculate for all three time periods
# by multiplying (weighing) each DA SSOS by its senior population
# then dividing the result by total senior population in Hamilton
da21_ham <- da21_ham %>%
  mutate(DA_SSOSw = (((DA_SSOS * srpop21)/srpop21total)*100),
         DA_SSOSw_5y = (((DA_SSOS * srpop26)/srpop26total)*100),
         DA_SSOSw_10y = (((DA_SSOS * srpop31)/srpop31total)*100))

# use mutate to calculate difference between Weighted SSOS now and 5y /10y
da21_ham <- da21_ham %>%
  mutate(DA_SSOSw_5y_chg = (DA_SSOSw_5y - DA_SSOSw),
         DA_SSOSw_10y_chg = (DA_SSOSw_10y - DA_SSOSw))

# check dataframe
head(da21_ham)

Simple feature collection with 6 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -79.94766 ymin: 43.26654 xmax: -79.88898 ymax: 43.30781
Geodetic CRS:  WGS 84
              DGUID    DAUID srpop21 srpop26 srpop31 cspaces_15
1 2021S051235250030 35250030      90     120     165          0
2 2021S051235250031 35250031     120     145     190          0
3 2021S051235250032 35250032     115     180     230          0
4 2021S051235250033 35250033     160     205     250          0
5 2021S051235250034 35250034     105     155     210          0
6 2021S051235250035 35250035      75     105     120          2
                        geometry  DA_SSOS  DA_SSOSw DA_SSOSw_5y
1 MULTIPOLYGON (((-79.92591 4... 0.000000 0.0000000   0.0000000
2 MULTIPOLYGON (((-79.89038 4... 0.000000 0.0000000   0.0000000
3 MULTIPOLYGON (((-79.93015 4... 0.000000 0.0000000   0.0000000
4 MULTIPOLYGON (((-79.93278 4... 0.000000 0.0000000   0.0000000
5 MULTIPOLYGON (((-79.94082 4... 0.000000 0.0000000   0.0000000
6 MULTIPOLYGON (((-79.94382 4... 2.061856 0.1483207   0.1522574
  DA_SSOSw_10y DA_SSOSw_5y_chg DA_SSOSw_10y_chg
1    0.0000000     0.000000000       0.00000000
2    0.0000000     0.000000000       0.00000000
3    0.0000000     0.000000000       0.00000000
4    0.0000000     0.000000000       0.00000000
5    0.0000000     0.000000000       0.00000000
6    0.1351335     0.003936721      -0.01318721

4.4 - Calculate Average Weighted Accessibility Score.

# use summarize
ham_scores_final <- da21_ham %>%
  st_drop_geometry() %>%
  summarise(DA_SSOS_avg = mean(DA_SSOS),
            DA_SSOSw_avg = mean(DA_SSOSw),
            DA_SSOSw_5y_avg = mean(DA_SSOSw_5y),
            DA_SSOSw_10y_avg = mean(DA_SSOSw_10y)) 

# Use mutate to calculate difference between Weighted SSOS now and 5y /10y
ham_scores_final <- ham_scores_final %>%
  mutate(DA_SSOSw_5y_chg = (DA_SSOSw_5y_avg - DA_SSOSw_avg),
         DA_SSOSw_10y_chg = (DA_SSOSw_10y_avg - DA_SSOSw_avg))

# check dataframe
ham_scores_final

PART 5 - Plot Results

5.1 - Number of Community Spaces Reached by Dissemination Area.

cspaces_reach <- tm_shape(da21_ham) +
  tm_fill("cspaces_15",
          style = "jenks",
          legend.hist = TRUE,
          title = "Community Spaces within 15 min.\n(Transit, max. 1000m walk, AM Peak)") +
  tm_layout(frame = FALSE,
            legend.outside = TRUE,
            legend.outside.position = 'bottom',
            legend.stack = 'horizontal',
            legend.title.fontface = 'bold',
            legend.hist.width = 1,
            legend.hist.height = 0.6) +
  tm_borders(col = "grey40", lwd = 0.1) +
  tm_legend(title.size=0.9,
            text.size = 0.6,
            position = c("left", "bottom")) +
  tm_shape(cd21_ham, color = 'black', fill = NA) +
  tm_borders(col = "black", lwd = 0.1) +
  tm_compass(north = 0, type = 'arrow', show.labels =0, position = c('right','bottom')) + 
  tm_scale_bar(color.dark = "black",
               position = c("left", "bottom"))
cspaces_reach

5.2 - Weighted Accessibility Score by Dissemination Area, 2021.

DA_SSOSw_plot <- tm_shape(da21_ham) +
  tm_polygons("DA_SSOSw", 
              style="jenks", 
              title="Weighted SSOS",
              border.alpha = 0) +
  tm_shape(cd21_ham, color = 'black', fill = NA) +
  tm_borders(col = "black", lwd = 0.1) +
  tm_layout(frame = FALSE) +
  tm_legend(position = c("right", "top")) +
  tm_compass(position = c('right','bottom')) + 
  tm_scale_bar(position = c("left", "bottom"), width = 0.15)
DA_SSOSw_plot

5.3 - Change in Weighted Accessibility Score by Dissemination Area, 2021 to 2026.

DA_SSOSw_5y_chg_plot <- tm_shape(da21_ham) +
  tm_polygons("DA_SSOSw_5y_chg", 
              style="jenks", 
              title="Five-Year\nChange (2026)",
              border.alpha = 0) +
  tm_shape(cd21_ham, color = 'black', fill = NA) +
  tm_borders(col = "black", lwd = 0.1) +
  tm_layout(frame = FALSE) +
  tm_legend(position = c("right", "top")) +
  tm_compass(position = c('right','bottom')) + 
  tm_scale_bar(position = c("left", "bottom"), width = 0.15)
DA_SSOSw_5y_chg_plot

5.4 - Change in Weighted Accessibility Score by Dissemination Area, 2021 to 2031.

DA_SSOSw_10y_chg_plot <- tm_shape(da21_ham) +
  tm_polygons("DA_SSOSw_10y_chg", 
              style="jenks", 
              title="Ten-Year\nChange (2031)",
              border.alpha = 0) +
  tm_shape(cd21_ham, color = 'black', fill = NA) +
  tm_borders(col = "black", lwd = 0.1) +
  tm_layout(frame = FALSE) +
  tm_legend(position = c("right", "top")) +
  tm_compass(position = c('right','bottom')) + 
  tm_scale_bar(position = c("left", "bottom"), width = 0.15)
DA_SSOSw_10y_chg_plot

Conclusion

Congratulations on making it to the end. Now run the same analysis, but this time, use a different facility data set and select a different population to study!

The figure below provides a nice visual summary of the analysis we just completed.

# 1 - DA Zones
DAplot <- ggplot() + 
  ggtitle("1) Load Dissemination Areas") +
  theme(plot.title = element_text(size = 10, hjust = 0, vjust = -1)) +
  geom_sf(data = da21_ham, color = 'grey', fill = NA, lwd = 0.3) +
  geom_sf(data = cd21_ham, color = 'black', fill = NA, lwd = 0.2) +
  blank()

# 2 - DA Centroids
DAGCplot <- ggplot() + 
  ggtitle("2) Calculate Centroids") +
  theme(plot.title = element_text(size = 10, hjust = 0, vjust = -1)) +
  geom_sf(data = da21_ham_gc, color = '#F8C471') +
  geom_sf(data = da21_ham, color = 'grey', fill = NA, lwd = 0.3) +
  geom_sf(data = cd21_ham, color = 'black', fill = NA, lwd = 0.2) +
  blank()

# 3 - 15min Centroid Isochrones
DAGCISOplot <- ggplot() + 
  ggtitle("3) Query OTP for 15min Isochrones") +
  theme(plot.title = element_text(size = 10, hjust = 0, vjust = -1)) +
  geom_sf(data = da21_ham_gc_iso, fill = '#F8C471', color = '#F8C471') +
  #geom_sf(data = da21_ham, color = 'grey', fill = NA, lwd = 0.3) +
  geom_sf(data = cd21_ham, color = 'black', fill = NA, lwd = 0.2) +
  blank()+
  annotation_scale(location = "bl",height = unit(0.1, "cm")) 

# 4 - Intersection Isochrones by Community Spaces
CSPACESplot <- ggplot() + 
  ggtitle("4) Intersect Community Spaces") +
  theme(plot.title = element_text(size = 10, hjust = 0, vjust = -1)) +
  geom_sf(data = da21_ham_gc_iso, fill = '#F8C471', color = '#F8C471') +
  #geom_sf(data = da21_ham, color = 'grey', fill = NA, lwd = 0.3) +
  geom_sf(data = cd21_ham, color = 'black', fill = NA, lwd = 0.2) +
  geom_sf(data = cspaces, color = '#A04000', shape=18, size=1.2) +
  blank()+
  annotation_north_arrow(location = "br", which_north = "true",
                         height = unit(0.8, "cm"),
                         width = unit(0.8, "cm"),
                         pad_y = unit(0.1, "in"),
                         style = north_arrow_fancy_orienteering)

# plot together
studyarea_method <- plot_grid(ncol = 2,
                              align = "hv",
                              DAplot, DAGCplot, DAGCISOplot, CSPACESplot,
                              label_size = 9) +
  theme(plot.title = element_text(hjust = 0.5, vjust = 1))
studyarea_method