# 7 Nearest Resource Analysis

## 7.1 Overview

Spatial Access to specific resource is often considered a multidimensional concept, though geographic distance is often central to the topic. Distance to the nearest resource is a common metric used to capture the availability of a resource, and in this tutorial we demonstrate how to calculate a minimum distance value from a ZCTA centroid to a set of resources, such as locations of methadone clinics. Each zip code will be assigned a “minimum distance access metric” as a value that indicates access to resources from that zip code. Our objectives are thus to:

• Generate centroids from areal data
• Calculate minimum distance from resources to area centroids
• Overlay resources and new minimum distance metric

## 7.2 Environment Setup

To replicate the codes & functions illustrated in this tutorial, you’ll need to have R and RStudio downloaded and installed on your system. This tutorial assumes some familiarity with the R programming language.

### 7.2.1 Packages used

We will use the following packages in this tutorial:

• `sf`: to manipulate spatial data
• `tmap`: to visualize and create maps
• `units`: to convert units within spatial data

### 7.2.2 Input/Output

Our inputs will be the clinic points we generated in previous tutorials, and the ZCTA boundary.

• Chicago Methadone Clinics, `methadone_clinics.shp`
• Chicago Zip Codes, `chicagoZips.shp`

We will calculate the minimum distance between the resources and the centroids of the zip codes, then save the results as a shapefile and as a CSV. Our final result will be a shapefile/CSV with the minimum distance value for each zip.

If you don’t have a shapefile of your data, but already have geographic coordinates as two columns in your CSV file, you can still use this tutorial. A reminder of how to transform your CSV with coordinates into a spatial data frame in R can be found here.

### 7.2.3 Load the packages

Load the libraries for use.

``````library(sf)
library(tmap)
library(units)``````

## 7.3 Load data

First, load in the MOUD resources shapefile. Let’s take a look at the first few rows of the dataset.

``methadoneSf <- st_read("data/methadone_clinics.shp")``
``````## Reading layer `methadone_clinics' from data source `/Users/yashbansal/Desktop/CSDS_RA/opioid-environment-toolkit/data/methadone_clinics.shp' using driver `ESRI Shapefile'
## Simple feature collection with 27 features and 8 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -87.7349 ymin: 41.68698 xmax: -87.57673 ymax: 41.96475
## geographic CRS: WGS 84``````
``head(methadoneSf)``
``````## Simple feature collection with 6 features and 8 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: -87.67818 ymin: 41.89268 xmax: -87.63409 ymax: 41.96475
## geographic CRS: WGS 84
##   X                                                         Name
## 1 1                Chicago Treatment and Counseling Center, Inc.
## 2 2                      Sundace Methadone Treatment Center, LLC
## 3 3 Soft Landing Interventions/DBA Symetria Recovery of Lakeview
## 4 4                                        PDSSC - Chicago, Inc.
## 5 5                          Center for Addictive Problems, Inc.
## 6 6                                Family Guidance Centers, Inc.
##                   Address    City State   Zip
## 1 4453 North Broadway st. Chicago    IL 60640
## 2 4545 North Broadway St. Chicago    IL 60640
## 3    3934 N. Lincoln Ave. Chicago    IL 60613
## 4     2260 N. Elston Ave. Chicago    IL 60614
## 5        609 N. Wells St. Chicago    IL 60654
## 6     310 W. Chicago Ave. Chicago    IL 60654
## 1 4453 North Broadway st. Chicago IL 60640        osm
## 2 4545 North Broadway St. Chicago IL 60640        osm
## 3    3934 N. Lincoln Ave. Chicago IL 60613     census
## 4     2260 N. Elston Ave. Chicago IL 60614     census
## 5        609 N. Wells St. Chicago IL 60654     census
## 6     310 W. Chicago Ave. Chicago IL 60654     census
##                     geometry
## 1 POINT (-87.65566 41.96321)
## 2 POINT (-87.65694 41.96475)
## 3 POINT (-87.67818 41.95331)
## 4 POINT (-87.67407 41.92269)
## 5 POINT (-87.63409 41.89268)
## 6 POINT (-87.63636 41.89657)``````

Next, we’ll load Chicago zip code boundaries.

``chicagoZips <- read_sf("data/chicago_zips.shp")``

We can quickly plot our data for to confirm they loaded correctly, here using an interactive map:

``````tm_shape(chicagoZips) +
tm_borders() +
tm_dots(col = "blue", size = 0.2)``````

## 7.4 Calculate centroids

Now, we will calculate the centroids of the zip code boundaries. We will first need to project our data, which means change it from latitude and longitude to meaningful units, like ft or meters, so we can calculate distance properly. We’ll use the Illinois State Plane projection, with an EPSG code of 3435.

Aside: To find the most appropriate projection for your data, do a Google Search for which projection works well - for state level data, each state has a State Plane projection with a specific code, known as the EPSG. I use epsg.io to check projections - here’s the New York State Plane page.

Use the `st_transform` function to change the projection of the data. Notice how the values in `geometry` go from being relatively small (unprojected, lat/long) to very large (projected, in US feet).

``````chicagoZips <- st_transform(chicagoZips, 3435)

chicagoZips``````
``````## Simple feature collection with 85 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 1058388 ymin: 1791133 xmax: 1205317 ymax: 1966816
## projected CRS:  NAD83 / Illinois East (ftUS)
## # A tibble: 85 x 10
##    ZCTA5CE10 GEOID10 CLASSFP10 MTFCC10 FUNCSTAT10 ALAND10 AWATER10
##  * <chr>     <chr>   <chr>     <chr>   <chr>      <chr>   <chr>
##  1 60501     60501   B5        G6350   S          125322… 974360
##  2 60007     60007   B5        G6350   S          364933… 917560
##  3 60651     60651   B5        G6350   S          9052862 0
##  4 60652     60652   B5        G6350   S          129878… 0
##  5 60653     60653   B5        G6350   S          6041418 1696670
##  6 60654     60654   B5        G6350   S          1464813 113471
##  7 60655     60655   B5        G6350   S          114080… 0
##  8 60656     60656   B5        G6350   S          8465226 0
##  9 60657     60657   B5        G6350   S          5888324 2025836
## 10 60659     60659   B5        G6350   S          5251086 2818
## # … with 75 more rows, and 3 more variables: INTPTLAT10 <chr>,
## #   INTPTLON10 <chr>, geometry <MULTIPOLYGON [US_survey_foot]>``````

Then, we will calculate the centroids:

``chicagoCentroids <- st_centroid(chicagoZips)``
``````## Warning in st_centroid.sf(chicagoZips): st_centroid assumes
## attributes are constant over geometries of x``````
``chicagoCentroids``
``````## Simple feature collection with 85 features and 9 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 1076716 ymin: 1802621 xmax: 1198093 ymax: 1956017
## projected CRS:  NAD83 / Illinois East (ftUS)
## # A tibble: 85 x 10
##    ZCTA5CE10 GEOID10 CLASSFP10 MTFCC10 FUNCSTAT10 ALAND10 AWATER10
##  * <chr>     <chr>   <chr>     <chr>   <chr>      <chr>   <chr>
##  1 60501     60501   B5        G6350   S          125322… 974360
##  2 60007     60007   B5        G6350   S          364933… 917560
##  3 60651     60651   B5        G6350   S          9052862 0
##  4 60652     60652   B5        G6350   S          129878… 0
##  5 60653     60653   B5        G6350   S          6041418 1696670
##  6 60654     60654   B5        G6350   S          1464813 113471
##  7 60655     60655   B5        G6350   S          114080… 0
##  8 60656     60656   B5        G6350   S          8465226 0
##  9 60657     60657   B5        G6350   S          5888324 2025836
## 10 60659     60659   B5        G6350   S          5251086 2818
## # … with 75 more rows, and 3 more variables: INTPTLAT10 <chr>,
## #   INTPTLON10 <chr>, geometry <POINT [US_survey_foot]>``````

For each zip code, this will calculate the centroid, and the output will be a point dataset.

Plot to double check that everything is ok. The `st_geometry()` function will once again just return the outline:

``````plot(st_geometry(chicagoZips))
plot(st_geometry(chicagoCentroids), add = TRUE, col = "red")``````

### 7.4.1 Visualize & Confirm

Once again, we can create an interactive map:

``````tm_shape(chicagoZips) +
tm_borders() +
tm_shape(chicagoCentroids) +
tm_dots()``````

## 7.5 Standardize CRS

If we immediately try to calculate the distance between the zip centroids and the locations of the resources using the `st_distance` function, we’ll get an error:

``st_distance(chicagoCentroids, methadoneSf, by_element = TRUE)``
``Error in st_distance(chicagoCentroids, methadoneSf, by_element = TRUE) : st_crs(x) == st_crs(y) is not TRUE``

Why is there an error? Because the projection of the centroids and the resource locations don’t match up. Let’s project the resource locations so that they match the projection of the centroids.

First, use the `st_crs` function to check that the coordinate reference system (or projection) is the same. They’re not, so we have to fix it.

``st_crs(chicagoCentroids)``
``````## Coordinate Reference System:
##   User input: EPSG:3435
##   wkt:
## PROJCRS["NAD83 / Illinois East (ftUS)",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4269]],
##     CONVERSION["SPCS83 Illinois East zone (US Survey feet)",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",36.6666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-88.3333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.999975,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",984250,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["USA - Illinois - SPCS - E"],
##         BBOX[37.06,-89.28,42.5,-87.02]],
##     ID["EPSG",3435]]``````
``st_crs(methadoneSf)``
``````## Coordinate Reference System:
##   User input: WGS 84
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]``````

We’ll take the CRS from the zip code centroids data, and use it as input to `st_transform` applied to the methadone clinics data.

``````newCRS <- st_crs(chicagoCentroids)
newCRS``````
``````## Coordinate Reference System:
##   User input: EPSG:3435
##   wkt:
## PROJCRS["NAD83 / Illinois East (ftUS)",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4269]],
##     CONVERSION["SPCS83 Illinois East zone (US Survey feet)",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",36.6666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-88.3333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.999975,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",984250,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["USA - Illinois - SPCS - E"],
##         BBOX[37.06,-89.28,42.5,-87.02]],
##     ID["EPSG",3435]]``````
``methadoneSf <- st_transform(methadoneSf, newCRS)``

If we check the CRS again, we now see that they match. Mismatched projections are a commonly made mistake in geospatial data processing.

``st_crs(chicagoCentroids)``
``````## Coordinate Reference System:
##   User input: EPSG:3435
##   wkt:
## PROJCRS["NAD83 / Illinois East (ftUS)",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4269]],
##     CONVERSION["SPCS83 Illinois East zone (US Survey feet)",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",36.6666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-88.3333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.999975,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",984250,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["USA - Illinois - SPCS - E"],
##         BBOX[37.06,-89.28,42.5,-87.02]],
##     ID["EPSG",3435]]``````
``st_crs(methadoneSf)``
``````## Coordinate Reference System:
##   User input: EPSG:3435
##   wkt:
## PROJCRS["NAD83 / Illinois East (ftUS)",
##         DATUM["North American Datum 1983",
##             ELLIPSOID["GRS 1980",6378137,298.257222101,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4269]],
##     CONVERSION["SPCS83 Illinois East zone (US Survey feet)",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",36.6666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",-88.3333333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",0.999975,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",984250,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["US survey foot",0.304800609601219],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["US survey foot",0.304800609601219]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["USA - Illinois - SPCS - E"],
##         BBOX[37.06,-89.28,42.5,-87.02]],
##     ID["EPSG",3435]]``````

Now we have the zip boundaries, the centroids of the zips, and the resource locations, as shown below. Next, we will calculate the distance to the nearest resource from each zip code centroid.

``````plot(st_geometry(chicagoZips))
plot(st_geometry(chicagoCentroids), col = "red", add = TRUE)
plot(st_geometry(methadoneSf), col = "blue", add = TRUE)``````

## 7.6 Calculate Distance

First, we’ll identify the resource that is the closest to a zip centroid using the `st_nearest_feature` function. (It will return the index of the object that is nearest, so we will subset the resources by the index to get the nearest object.)

``````nearestClinic_indexes <- st_nearest_feature(chicagoCentroids, methadoneSf)

nearestClinic``````
``````## Simple feature collection with 85 features and 8 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 1147259 ymin: 1829330 xmax: 1190680 ymax: 1930471
## projected CRS:  NAD83 / Illinois East (ftUS)
## First 10 features:
##       X
## 16   16
## 7     7
## 7.1   7
## 26   26
## 15   15
## 5     5
## 26.1 26
## 7.2   7
## 1     1
## 3     3
##                                                              Name
## 16                            Katherine Boone Robinson Foundation
## 7                                        A Rincon Family Services
## 7.1                                      A Rincon Family Services
## 26                              New Hope Community Service Center
## 15           HRDI- Grand Boulevard Professional Counseling Center
## 5                             Center for Addictive Problems, Inc.
## 26.1                            New Hope Community Service Center
## 7.2                                      A Rincon Family Services
## 1                   Chicago Treatment and Counseling Center, Inc.
## 3    Soft Landing Interventions/DBA Symetria Recovery of Lakeview
##                      Address    City State   Zip
## 16        4100 W. Ogden Ave. Chicago    IL 60623
## 7         3809 W. Grand Ave. Chicago    IL 60651
## 7.1       3809 W. Grand Ave. Chicago    IL 60651
## 26          2559 W. 79th St. Chicago    IL 60652
## 15           340 E. 51st St. Chicago    IL 60615
## 5           609 N. Wells St. Chicago    IL 60654
## 26.1        2559 W. 79th St. Chicago    IL 60652
## 7.2       3809 W. Grand Ave. Chicago    IL 60651
## 1    4453 North Broadway st. Chicago    IL 60640
## 3       3934 N. Lincoln Ave. Chicago    IL 60613
## 16        4100 W. Ogden Ave. Chicago IL 60623     census
## 7         3809 W. Grand Ave. Chicago IL 60651     census
## 7.1       3809 W. Grand Ave. Chicago IL 60651     census
## 26          2559 W. 79th St. Chicago IL 60652        osm
## 15           340 E. 51st St. Chicago IL 60615     census
## 5           609 N. Wells St. Chicago IL 60654     census
## 26.1        2559 W. 79th St. Chicago IL 60652        osm
## 7.2       3809 W. Grand Ave. Chicago IL 60651     census
## 1    4453 North Broadway st. Chicago IL 60640        osm
## 3       3934 N. Lincoln Ave. Chicago IL 60613     census
##                     geometry
## 16   POINT (1149437 1888620)
## 7    POINT (1150707 1908328)
## 7.1  POINT (1150707 1908328)
## 26   POINT (1160422 1852071)
## 15   POINT (1179319 1871281)
## 5    POINT (1174632 1904257)
## 26.1 POINT (1160422 1852071)
## 7.2  POINT (1150707 1908328)
## 1    POINT (1168556 1929911)
## 3    POINT (1162460 1926257)``````

Then, we will calculate the distance between the nearest resource and the zip code centroid with the `st_distance` function. As shown above, make sure both of your datasets are projected, and in the same projection, before you run `st_distance`.

``````minDist <- st_distance(chicagoCentroids, nearestClinic, by_element = TRUE)

minDist``````
``````## Units: [US_survey_foot]
##  [1] 36764.471 82821.757  5237.764  7419.050  7241.215   873.991
##  [7] 20515.835 38337.779  8516.342 15479.988  9823.172  4447.705
## [13] 33707.202 24111.173 24184.991 45224.230 31198.648 10113.372
## [19] 10890.143 13779.075 49874.400 32463.844 36633.404 21981.243
## [25] 13637.327 22151.171 63268.911  4242.008  3749.998  5116.415
## [31]  5530.724  8696.602  3965.649  4368.285  6998.725  6937.025
## [37]  3967.509 45944.187 32598.459 44548.952 58483.551  5416.487
## [43]  5728.851  2327.899  6646.411  5817.262   365.628 13619.270
## [49]  7021.182  5024.536  2829.488  1648.361  7709.690  2628.400
## [55]  1521.495  9423.848 16648.775  2185.237 11463.785 22537.644
## [61] 11872.133  7759.662 39665.409 11741.279 18689.070 27555.396
## [67]  5301.820  8808.542 25363.914 18986.590 26456.564  7562.917
## [73]  3496.941 11023.379  3128.000 16826.923  6172.874 25276.556
## [79] 17290.797 15198.524 25164.072 18438.696 19972.965 33211.229
## [85] 22518.257``````

This is in US feet. To change to a more meaningful unit, such as miles, we can use the `set_units()` function:

``````minDist_mi <- set_units(minDist, "mi")

minDist_mi``````
``````## Units: [mi]
##  [1]  6.96298198 15.68597011  0.99200275  1.40512596  1.37144503
##  [6]  0.16552893  3.88558262  7.26095754  1.61294677  2.93182182
## [11]  1.86045293  0.84236999  6.38395252  4.56651915  4.58049994
## [16]  8.56521226  5.90884663  1.91541520  2.06253115  2.60967860
## [21]  9.44592795  6.14846754  6.93815858  4.16312245  2.58283224
## [26]  4.19530595 11.98277234  0.80341215  0.71022833  0.96901998
## [31]  1.04748778  1.64708702  0.75107138  0.82732834  1.32551867
## [36]  1.31383303  0.75142366  8.70156788  6.17396283  8.43731838
## [41] 11.07645233  1.02585189  1.08501190  0.44089088  1.25879245
## [46]  1.10175645  0.06924787  2.57941243  1.32977202  0.95161851
## [51]  0.53588899  0.31219014  1.46017146  0.49780409  0.28816254
## [56]  1.78482335  3.15318335  0.41387148  2.17117583  4.26850167
## [61]  2.24851459  1.46963585  7.51240307  2.22373150  3.53960367
## [66]  5.21883539  1.00413457  1.66828778  4.80378114  3.59595236
## [71]  5.01072284  1.43237347  0.66230081  2.08776534  0.59242542
## [76]  3.18692368  1.16910723  4.78723618  3.27477863  2.87851400
## [81]  4.76593224  3.49218425  3.78276617  6.29001804  4.26482998``````

### 7.6.1 Merge Data

We then rejoin the minimum distances to the zip code data, by column binding `minDist_mi` to the original `chicagoZips` data.

``````minDistSf <- cbind(chicagoZips, minDist_mi)
minDistSf``````
``````## Simple feature collection with 85 features and 10 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 1058388 ymin: 1791133 xmax: 1205317 ymax: 1966816
## projected CRS:  NAD83 / Illinois East (ftUS)
## First 10 features:
##    ZCTA5CE10 GEOID10 CLASSFP10 MTFCC10 FUNCSTAT10  ALAND10
## 1      60501   60501        B5   G6350          S 12532295
## 2      60007   60007        B5   G6350          S 36493383
## 3      60651   60651        B5   G6350          S  9052862
## 4      60652   60652        B5   G6350          S 12987857
## 5      60653   60653        B5   G6350          S  6041418
## 6      60654   60654        B5   G6350          S  1464813
## 7      60655   60655        B5   G6350          S 11408010
## 8      60656   60656        B5   G6350          S  8465226
## 9      60657   60657        B5   G6350          S  5888324
## 10     60659   60659        B5   G6350          S  5251086
##    AWATER10  INTPTLAT10   INTPTLON10      minDist_mi
## 1    974360 +41.7802209 -087.8232440  6.9629820 [mi]
## 2    917560 +42.0086000 -087.9973398 15.6859701 [mi]
## 3         0 +41.9020934 -087.7408565  0.9920027 [mi]
## 4         0 +41.7479319 -087.7147951  1.4051260 [mi]
## 5   1696670 +41.8199645 -087.6059654  1.3714450 [mi]
## 6    113471 +41.8918225 -087.6383036  0.1655289 [mi]
## 7         0 +41.6947762 -087.7037764  3.8855826 [mi]
## 8         0 +41.9742800 -087.8271283  7.2609575 [mi]
## 9   2025836 +41.9402931 -087.6468569  1.6129468 [mi]
## 10     2818 +41.9914885 -087.7039859  2.9318218 [mi]
##                          geometry
## 1  MULTIPOLYGON (((1112613 185...
## 2  MULTIPOLYGON (((1058389 194...
## 3  MULTIPOLYGON (((1136069 190...
## 4  MULTIPOLYGON (((1145542 185...
## 5  MULTIPOLYGON (((1177007 187...
## 6  MULTIPOLYGON (((1170904 190...
## 7  MULTIPOLYGON (((1146378 183...
## 8  MULTIPOLYGON (((1110359 193...
## 9  MULTIPOLYGON (((1162394 192...
## 10 MULTIPOLYGON (((1148555 194...``````

### 7.6.2 Visualize & Confirm

We can now visualize the zip-level access to methadone clinics using our new access metric, using the `tmap` package. We’ll use quantile bins with five intervals.

``tmap_mode("plot")``
``## tmap mode set to plotting``
``````tm_shape(minDistSf) +
tm_polygons("minDist_mi", style = 'quantile', n=5,
title = "Minimum Distance (mi)") +
tm_layout(main.title = "Minimum Distance from Zip Centroid\n to Methadone Clinic",
main.title.position = "center",
main.title.size = 1)``````

Access by zip code can also be combined with locations of resources:

``````tm_shape(minDistSf) +
tm_polygons("minDist_mi", style = 'quantile', n=5,
title = "Minimum Distance (mi)") +
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)``````

``head(chicagoZips)``
``````## Simple feature collection with 6 features and 9 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 1058388 ymin: 1846408 xmax: 1189038 ymax: 1957548
## projected CRS:  NAD83 / Illinois East (ftUS)
## # A tibble: 6 x 10
##   ZCTA5CE10 GEOID10 CLASSFP10 MTFCC10 FUNCSTAT10 ALAND10 AWATER10
##   <chr>     <chr>   <chr>     <chr>   <chr>      <chr>   <chr>
## 1 60501     60501   B5        G6350   S          125322… 974360
## 2 60007     60007   B5        G6350   S          364933… 917560
## 3 60651     60651   B5        G6350   S          9052862 0
## 4 60652     60652   B5        G6350   S          129878… 0
## 5 60653     60653   B5        G6350   S          6041418 1696670
## 6 60654     60654   B5        G6350   S          1464813 113471
## # … with 3 more variables: INTPTLAT10 <chr>, INTPTLON10 <chr>,
## #   geometry <MULTIPOLYGON [US_survey_foot]>``````

## 7.7 Save Data

To save our final result to a CSV:

``write_sf(minDistSf, "chizips_Access.csv")``

We can also write out this data to a shapefile format:

``write_sf(minDistSf, "chizips_Access.shp")``