The Center for Spatial Data Science at UChicago created several tutorials for their Opioid Environment Toolkit focusing on how methadone clinics serve the chicago area. This code is adapted from their toolkit for use in the HEAL platform.
The Opioid Environment Policy Scan (OEPS) is a database providing access to data at multiple spatial scales to help characterize the multi-dimensional risk environment impacting opioid use in justice populations across the United States. The OEPS and the Opioid Environment Toolkit, from which this script was adapted, were developed for the JCOIN network by Marynia Kolak, Qinyun Lin, Susan Paykin, Moksha Menghaney, and Angela Li at the Center for Spatial Data Science at the University of Chicago as part of the Methodology and Advanced Analytics Resource Center (MAARC).
Citation:
Marynia Kolak, Qinyun Lin, Susan Paykin, Moksha Menghaney, & Angela Li. (2021, May 11). GeoDaCenter/opioid-policy-scan: Opioid Environment Policy Scan Data Warehouse (Version v0.1-beta). Zenodo. http://doi.org/10.5281/zenodo.4747876
While much of the code is directly from the Toolkit tutorial, I also expanded to make it usable within the HEAL platform and within a R Kernel Jupyter Notebook environment and to learn a few additional functions to quantify results.
This notebook shows how to use both locally uploaded files in combination with HEAL platform specific files.
Specifically, COVID data has been uploaded locally while geometric and clinic files are taken from the OEPS database on the heal platform.
#suppress warnings (to include warnings,set warn to 0)
options(warn=-1)
library(sf)
library(tidygeocoder)
library(tmap)
library(tidyverse)
library(readxl)
library(units)
# Pull file objects using Gen3 SDK
system("gen3 drs-pull object dg.6VTS/eccae5e6-d9be-403a-b30f-4df1434a2921")
system("gen3 drs-pull object dg.6VTS/f30a3855-9093-473f-8ce6-784fc9b53066")
system("gen3 drs-pull object dg.6VTS/1eb8e61b-328f-4cd0-a7e2-df2f7deec846")
system("gen3 drs-pull object dg.6VTS/e4c13f4f-7466-4631-bf9e-c83710abcd32")
system("gen3 drs-pull object dg.6VTS/33b0fa36-f58d-4a49-a1f7-376376c4ffb1")
system("gen3 drs-pull object dg.6VTS/12d90fed-c2d7-402d-9fa6-abd78142191e")
#local data on HEAL platform
local_data_dir = '/home/jovyan/local_data'
A common goal in opioid environment research is to calculate and compare access metrics to different providers of Medications for Opioid Overuse Disorder (MOUDs). Before we can run any analytics on the resource location data, we need to convert resource addresses to spatial data points, which can be then used to calculate access metrics.
Geocoding is the process of converting addresses (like a street address) into geographic coordinates using a known coordinate reference system. We can then use these coordinates (latitude, longitude) to spatially enable our data. This means we convert to a spatial data frame (sf) within R for spatial analysis within our R session, and then save as a shapefile (a spatial data format) for future use. In this tutorial we demonstrate how to geocode resource location addresses and convert to spatial data points that can be used for future mapping and geospatial analysis.
#get zip - county crosswalk to filter zip codes for Cook County
#Cook County (county Chicago is in) FIPS code = 17-031
cook_county_zips <- read_excel("./ZIP_COUNTY.xlsx") %>%
filter(COUNTY==17031)
allClinics <- read.csv('./us-wide-moudsCleaned.csv')
names(allClinics) <-c(
"name1",
"name2" ,
"street1" ,
"street2" ,
"city" ,
"state" ,
"zip" ,
"zip4" ,
"category",
"countyGEOID" ,
"countyName" ,
"source",
"geom1" ,
"geom2",
"Longitude",
"Latitude"
)
methadoneClinics <- allClinics %>%
filter(countyGEOID==17031 & category=='methadone')
#project to coordinate system
#longitude and latitude contrary to common lat/long
methadoneSf <- st_as_sf(methadoneClinics,
coords = c("Longitude", "Latitude"),
crs = 4326)
#geocode projection in geometry field
head(tibble(methadoneSf),3)
name1 | name2 | street1 | street2 | city | state | zip | zip4 | category | countyGEOID | countyName | source | geom1 | geom2 | geometry |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
<chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <chr> | <int> | <chr> | <int> | <chr> | <chr> | <chr> | <chr> | <POINT [°]> |
Elite Houses of Sober Living | Elite Treatment Center | 395 West Lincoln Highway | NA | Chicago Heights | IL | 60411 | NA | methadone | 17031 | Cook | SAMHSA | c(-87.6576384 | 41.5062933) | POINT (-87.65764 41.50629) |
Radius Foundation Inc | NA | 11952 South Harlem Avenue | Suite 100 | Palos Heights | IL | 60463 | NA | methadone | 17031 | Cook | SAMHSA | c(-87.7972744 | 41.6737572) | POINT (-87.79727 41.67376) |
Symetria Recovery | Palos Heights | 11925 South Harlem Avenue | NA | Palos Heights | IL | 60463 | 1138 | methadone | 17031 | Cook | SAMHSA | c(-87.7963882 | 41.6745493) | POINT (-87.79639 41.67455) |
Show projection of points on interactive Open Street map:
tmap_mode("view")
## tmap mode set to interactive viewing
#notebook edit: not rendering interactive so need to save rendered hmtl as variable with print fxn
clinic_map <- tm_shape(methadoneSf) +
tm_dots() +
tm_basemap("OpenStreetMap") +
tm_view(
set.bounds=c(-87.94011,41.64454,-87.52414,42.02304),
alpha = 1
)
#tmap_save(clinic_map,'clinic_map')
clinic_rendered <- print(clinic_map)
tmap_save(tm=clinic_map,file='clinic_render.html')
#clinic_rendered
tmap mode set to interactive viewing Interactive map saved to /home/jovyan/clinic_render.html
Once we have spatially referenced resource locations, it's helpful to plot the data in the community of interest for some preliminary analysis. In this tutorial we will plot Methadone Providers in Chicago and community areas to provide some context. We will also generate a simple 1-mile buffer service area around each provider to highlight neighborhoods with better, and worse, access to resources. In order to accomplish this task, we will need to standardize our spatial data (clinic points, and community areas) with an appropriate coordinate reference system. Finally, we'll make some maps!
Our objectives are thus to:
shapefile <- read_sf('./zctas2018.shp')
#filter shapefile with only cook county zips
areas <- shapefile %>% filter(ZCTA5CE10 %in% cook_county_zips$ZIP)
#local file
##can get directly from github or uploaded locally
cityBoundary <- st_read(paste0(local_data_dir,"/boundaries_chicago.geojson"))
Reading layer `boundaries_chicago' from data source `/home/jovyan/local_data/boundaries_chicago.geojson' using driver `GeoJSON' Simple feature collection with 1 feature and 4 fields Geometry type: MULTIPOLYGON Dimension: XY Bounding box: xmin: -87.94011 ymin: 41.64454 xmax: -87.52414 ymax: 42.02304 Geodetic CRS: WGS 84
#plot to show geocoding performed as desired
tmap_mode("plot")
## 1st layer (gets plotted first)
tm_shape(areas) +
tm_borders(alpha = 0.4) +
## 2nd layer (overlay)
tm_shape(methadoneSf) +
tm_dots(size = 0.4, col="red")
tmap mode set to plotting
#coordinate reference system that preserves distance
CRS.new <- st_crs("EPSG:3435")
#transform datasets to new CRS
metClinics.3435 <- st_transform(methadoneSf, CRS.new)
areas.3435 <- st_transform(areas, CRS.new)
#generate buffer with new CRS but use 5280 ft to equal 1 mile
metClinic_buffers <- st_buffer(metClinics.3435, 5280)
#plot clinic points and 1 mile buffers around clinic
tmap_mode("plot")
tm_shape(areas.3435) +
tm_borders() +
tm_shape(metClinics.3435) +
tm_dots(col = "red") +
tm_shape(metClinic_buffers) +
tm_borders(col = "blue")
tmap mode set to plotting
#fill in buffers for easier viz
tm_shape(areas) +
tm_borders(alpha = 0.6) +
tm_shape(metClinic_buffers) +
tm_fill(col = "blue", alpha = .4) +
tm_borders(col = "blue") +
tm_shape(metClinics.3435) +
tm_dots(col = "red", size = 0.2)
#does an area have coverage (ie falls within any one area?)
unionBuffers <- st_union(metClinic_buffers)
tm_shape(areas) +
tm_borders()+
tm_shape(unionBuffers) +
tm_fill(col = "blue", alpha = .2) +
tm_borders(col = "blue") +
tm_shape(metClinics.3435) +
tm_dots(col = "red", size = 0.1)+
tm_shape(cityBoundary)+
tm_borders(lwd =3)
metClinic_2mbuffers <- st_buffer(metClinics.3435, 10560)
tmap_mode("plot")
tm_shape(areas) +
tm_borders(alpha = 0.4) +
tm_shape(cityBoundary) +
tm_borders(lwd = 3) +
tm_shape(metClinic_2mbuffers) +
tm_fill(col = "gray10", alpha = .4) +
tm_borders(col = "dimgray", alpha = .4) +
tm_shape(metClinic_buffers) +
tm_fill(col = "gray90", alpha = .4) +
tm_borders(col = "darkslategray") +
tm_shape(metClinics.3435) +
tm_dots(col = "red", size = 0.2) +
tm_layout(main.title = "Methadone Clinic Service Areas in Cook County and Chicago",
main.title.position = "center",
main.title.size = 1,
frame = FALSE)
tmap mode set to plotting
#does an area have coverage (ie falls within any one area?)
union2mBuffers <- st_union(metClinic_2mbuffers)
tmap_mode("plot")
tm_shape(areas) +
tm_borders(alpha = 0.4) +
tm_shape(cityBoundary) +
tm_borders(lwd = 3) +
tm_shape(union2mBuffers) +
tm_fill(col = "gray10", alpha = .2) +
tm_borders(col = "dimgray", alpha = .2) +
tm_shape(unionBuffers) +
tm_fill(col = "gray90", alpha = .4) +
tm_borders(col = "darkslategray") +
tm_shape(metClinics.3435) +
tm_dots(col = "red", size = 0.2) +
tm_layout(main.title = "Methadone Clinic Service Areas in Cook County/Chicago",
main.title.position = "center",
main.title.size = 1,
frame = FALSE)
tmap mode set to plotting
tm_shape(areas,bbox = cityBoundary) + #zoom in on just chicago
tm_borders(alpha = 0.5, col="gray")+
tm_text("GEOID10", size = 0.6) +
tm_shape(cityBoundary) +
tm_borders() +
tm_shape(unionBuffers) +
tm_fill(col = "blue", alpha = .2) +
tm_borders(col = "blue") +
tm_shape(metClinics.3435) +
tm_dots(col = "red")
Can do conditional buffers:
You can then quantify what you are seeing above by creating a simple feature object from the intersection of the city and the buffers with the overall city boundaries using st_intersection
and st_area
#apply the CRS of chicago Centroids to methadone to ensure they are both in same CRS
newCRS <- st_crs(cityBoundary)
unionBuffers_newCRS <- st_transform(unionBuffers, newCRS)
unionBuffers2mi_newCRS <- st_transform(union2mBuffers, newCRS)
#calculate area of chicago with chicago boundary
city_area <- st_area(cityBoundary)
#calculate area of 1 mi buffer and chicago boundary intersection
buffers_and_city_intersection <- st_intersection(unionBuffers_newCRS,cityBoundary)
buffers_2mi_and_city_intersection <- st_intersection(unionBuffers2mi_newCRS,cityBoundary)
#calculate area of 2 mi buffer and chicago boundary intersection
buffer_area <- st_area(buffers_and_city_intersection)
buffer2mi_area <- st_area(buffers_2mi_and_city_intersection)
paste0(
as.character(
as.integer(
(buffer_area/city_area)*100)),
"% within 1 mile of a resource across Chicago"
)
paste0(
as.character(
as.integer(
(buffer2mi_area/city_area)*100)),
"% within 2 miles of a resource across Chicago"
)
link additional pieces of information with the buffer analysis to target key vulnerable populations. The OEPS tutorial uses the COVID pandemic as a case study.
connect Chicago COVID-19 Case data by ZIP Code, available as a flat file on the city’s data portal, to our environment.
overlap the 1-mile buffers representing walkable access to the Methadone providers in the city. COVID impacts travel so conservative, "walkable" threshold used.
covid_dir <- paste0(local_data_dir,"/COVID-19_Cases__Tests__and_Deaths_by_ZIP_Code.csv")
#read in data for COVID info by zip code
#local file
COVID <- read.csv(covid_dir) %>%
select(c("ZIP.Code", "Case.Rate...Cumulative")) %>%
mutate(GEOID10=as.character(ZIP.Code))
#merge COVID with zips
zipsMerged <- merge(areas, COVID, by = "GEOID10")
tmap_mode("plot")
tm_shape(zipsMerged,bbox=cityBoundary) +
tm_fill(
"Case.Rate...Cumulative", style="quantile", n=6, pal="-RdYlGn",
title = "COVID Case Rate",alpha = 0.8
) +
tm_borders(lwd = 0) +
tm_shape(metClinic_buffers,bbox=cityBoundary) +
tm_borders(col = "gray") +
tm_fill(alpha=0.5) +
tm_shape(methadoneSf,bbox=cityBoundary) +
tm_dots(col = "black", size = 0.1) +
tm_layout(main.title = "Walkable Methadone Service Areas\nin and around Chicago",
main.title.position = "center",
main.title.size = 1.5,
frame = FALSE)
tmap mode set to plotting
tmap_mode("view")
community_interactive <-
tm_shape(zipsMerged) +
tm_fill(
"Case.Rate...Cumulative", style="quantile", n=6, pal="-RdYlGn",
title = "COVID Case Rate",alpha = 0.8
) +
tm_borders(lwd = 0) +
tm_shape(metClinic_buffers) +
tm_borders(col = "gray") +
tm_fill(alpha=0.5) +
tm_shape(methadoneSf) +
tm_dots(col = "black", size = 0.01) +
tm_layout(main.title = "Walkable Methadone Service Areas",
main.title.position = "center",
main.title.size = 1,
frame = FALSE)
community_interactive_rendered <- print(community_interactive)
tmap_save(tm=community_interactive,file='covid_render.html')
tmap mode set to interactive viewing Interactive map saved to /home/jovyan/covid_render.html
Of course, this visualization needs a more complete picture of opioid outcomes by zip. But, the red zip codes far from walking distance may be at risk for increased opioid overdoses without intervention (ie from pop up shops).
This visualization illustrates that many of the high risk areas do not have walkable access to a clinic for the majority of that area.
Obviously, other pieces of community info need to be linked (such as baseline use of opioids for these areas to determine risk) -- see intro for recommendations
Perhaps an action item could be to set up "pop up shop" methadone clinics.
#change projection of data to a new CRS
chicagoZips <- st_transform(areas, 3435)
#calculate centroids for each zip code region
chicagoCentroids <- st_centroid(chicagoZips)
#show the centroids
plot(st_geometry(chicagoZips))
plot(st_geometry(chicagoCentroids), add = TRUE, col = "red")
#apply the CRS of chicago Centroids to methadone to ensure they are both in same CRS
newCRS <- st_crs(chicagoCentroids)
methadoneSf <- st_transform(methadoneSf, newCRS)
#now plot both the methadone clinics and the zip centroids
plot(st_geometry(chicagoZips))
plot(st_geometry(chicagoCentroids), col = "red", add = TRUE)
plot(st_geometry(methadoneSf), col = "blue", add = TRUE)
#identify the nearest clinic to each zip centroid
nearestClinic_indexes <- st_nearest_feature(chicagoCentroids, methadoneSf)
nearestClinic <- methadoneSf[nearestClinic_indexes,]
#calculate distance to nearest clinic
minDist <- st_distance(chicagoCentroids, nearestClinic, by_element = TRUE)
#set units to miles
minDist_mi <- set_units(minDist, "mi")
#join zip codes with the minimum distance data
minDistSf <- cbind(chicagoZips, minDist_mi)
#visualize minimum distance to each centroid
tmap_mode("plot")
tm_shape(minDistSf) +
tm_polygons("minDist_mi", style = 'quantile', n=5,
title = "Minimum Distance (mi)") +
tm_shape(methadoneSf) +
tm_dots(size = 0.2) +
tm_layout(main.title = "Minimum Distance from Zip Centroid\n to Methadone Clinic",
main.title.position ="center",
frame = FALSE,
main.title.size = 1)
tmap mode set to plotting
To do for future mapping other data