Isochrone Mapping
Jonathan Simkin
2022-10-27
Intro
- Drive times are an important indicator/exposure in health care analyses / research as a proxy for access to care (specifically spatial access).
- Isochrone maps are a very cool way of visualizing drive times and assigning exposure to data points.
- In this doc, I’ll show what they are and how to produce them in R
Objective
- To describe access to care through an isochrone map
Analysis
Libraries
Here are some R packages we need
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.0 ✔ stringr 1.4.1
## ✔ readr 2.1.2 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(sf)## Linking to GEOS 3.10.2, GDAL 3.4.2, PROJ 8.2.1; sf_use_s2() is TRUElibrary(osrm)## Data: (c) OpenStreetMap contributors, ODbL 1.0 - http://www.openstreetmap.org/copyright
## Routing: OSRM - http://project-osrm.org/library(tmap)Our cohort
Here’s a random set of point locations in Berling - we are using the apotheke spatial file that comes with this package
Let’s load the data
apotheke.sf <- st_read(system.file("gpkg/apotheke.gpkg", package = "osrm"),
quiet = TRUE)Let’s map these locations
tmap_mode("view")## tmap mode set to interactive viewingtm_shape(apotheke.sf) +
tm_dots(size = 0.1,
col = "black")How close are they to a health care facility?
We’re going to pick one health facility and in this example, we’ll pick VGH.
You can get the location from google maps
health_location <-
tibble(
case_id = c(1),
long = c(13.406638967397514
),
lat = c(
52.526597905944406
)
) %>%
st_as_sf(., coords = c("long", "lat"), crs = 4326, agr = "constant") %>%
st_transform(., crs = 4326)
tm_shape(apotheke.sf) +
tm_dots(size = 0.1,
col = "black") +
tm_shape(health_location) +
tm_dots(size = 0.2,
col = "blue") Exposure assessment
Here we want to measure an exposure… and we will choose travel time on the road (driving a car)
And we’ll do that with the osrm R package
iso_map <-
osrmIsochrone(loc = health_location,
breaks = seq(from = 0, to = 60, 5),
res = 30
)
tm_shape(iso_map %>% rename(drive_time = isomax)) +
tm_polygons("drive_time",
palette = "viridis",
alpha = 0.8)And now let’s add the health events to the map
tm_shape(iso_map %>% rename(drive_time_max = isomax)) +
tm_polygons("drive_time_max",
n = 10,
palette = "viridis",
alpha = 0.6) +
tm_shape(apotheke.sf) +
tm_dots(size = 0.1,
col = "black") +
tm_shape(health_location) +
tm_dots(size = 0.2,
col = "blue") Linking spatial data
Now we can link data from one map to the other. Let’s link the drive time categories to the point locations
int <-
st_intersects(apotheke.sf %>% st_transform(crs = 4326), iso_map, prepared = TRUE) #Wherever there are intersection between the two layers (the CHSA Map and the study data), we can link them! and store them in "int"
linked_data <-
apotheke.sf %>%
mutate(
drive_time_max = as.character(iso_map$isomax[as.numeric(as.character(int))]),
drive_time_min = as.character(iso_map$isomin[as.numeric(as.character(int))])
) %>%
st_drop_geometry()
linked_data## id drive_time_max drive_time_min
## 1 440338666 20 15
## 2 538057637 30 25
## 3 977657079 25 20
## 4 3770254015 20 15
## 5 364363337 15 10
## 6 3666540067 20 15
## 7 4234162689 15 10
## 8 1492336686 10 0
## 9 1298511709 15 10
## 10 268726812 15 10
## 11 371157962 10 0
## 12 332555440 15 10
## 13 1632227368 35 30
## 14 573452398 20 15
## 15 2424327600 30 25
## 16 715203882 20 15
## 17 1869023719 20 15
## 18 4936408725 35 30
## 19 565034853 25 20
## 20 2476268050 20 15
## 21 4010310541 25 20
## 22 392352485 25 20
## 23 1400224784 25 20
## 24 308719063 30 25
## 25 440339696 20 15
## 26 576608994 15 10
## 27 76596388 10 0
## 28 571800623 30 25
## 29 3343659161 15 10
## 30 473696292 35 30
## 31 674086086 35 30
## 32 1124612201 35 30
## 33 703042593 30 25
## 34 336117213 20 15
## 35 3073150198 10 0
## 36 1539485716 20 15
## 37 2713665583 20 15
## 38 287852225 25 20
## 39 1914403852 20 15
## 40 609146610 30 25
## 41 3012143573 25 20
## 42 1382999356 25 20
## 43 2516120286 15 10
## 44 833390882 25 20
## 45 769365479 15 10
## 46 2574350272 20 15
## 47 81362602 25 20
## 48 671210857 25 20
## 49 1963129662 20 15
## 50 1641113914 25 20
## 51 3337776414 35 30
## 52 646579859 15 10
## 53 428114308 20 15
## 54 270906159 45 40
## 55 282744493 15 10
## 56 459166572 15 10
## 57 765341424 15 10
## 58 602531459 15 10
## 59 3775184617 15 10
## 60 448818334 15 10
## 61 664244049 25 20
## 62 465560644 15 10
## 63 1397070332 20 15
## 64 440330654 15 10
## 65 671390314 20 15
## 66 2361562308 20 15
## 67 272385919 35 30
## 68 3353236089 35 30
## 69 472372248 25 20
## 70 3157917186 25 20
## 71 499046359 20 15
## 72 469201022 15 10
## 73 670658040 15 10
## 74 3632641432 20 15
## 75 3338651114 35 30
## 76 582854286 30 25
## 77 4396193589 40 35
## 78 645677466 20 15
## 79 1833509167 25 20
## 80 588552733 15 10
## 81 459785177 20 15
## 82 2288123586 15 10
## 83 366186442 20 15
## 84 1828707071 15 10
## 85 304367397 40 35
## 86 428114310 20 15
## 87 314752277 15 10
## 88 419873047 30 25
## 89 267860831 35 30
## 90 1365702036 30 25
## 91 2484338009 15 10
## 92 2722616195 15 10
## 93 661126972 20 15
## 94 609142032 30 25
## 95 2864491450 30 25
## 96 1169920988 25 20
## 97 1419134521 25 20
## 98 506338850 20 15
## 99 441514950 20 15
## 100 1336607982 30 25Let’s do some simple stats
linked_data %>%
count(drive_time_max)## drive_time_max n
## 1 10 4
## 2 15 26
## 3 20 27
## 4 25 18
## 5 30 12
## 6 35 10
## 7 40 2
## 8 45 1ggplot(linked_data, aes(x = drive_time_max)) +
geom_histogram(stat = "count", colour = "black", alpha =0.8) +
theme_bw() +
labs(
x = "Drive time from health care location",
y = "Count",
title = "Distribution of drive time (minutes) to health care location"
)## Warning: Ignoring unknown parameters: binwidth, bins, pad