4  Spatial Stuff

The most essential packages for doing GIS work (with ESRI products) are:

Honorable mention:

library(sf)
Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE
library(mapview)
library(arcgis)
Attaching core arcgis packages:
→ arcgisutils v0.1.1.9000
→ arcgislayers v0.1.0
library(tidyverse)
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
✔ ggplot2 3.4.0     ✔ purrr   0.3.4
✔ tibble  3.1.8     ✔ dplyr   1.1.0
✔ tidyr   1.2.1     ✔ stringr 1.5.0
✔ readr   2.1.3     ✔ forcats 0.5.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ purrr::%||%()    masks arcgisutils::%||%()
✖ purrr::compact() masks arcgisutils::compact()
✖ dplyr::filter()  masks stats::filter()
✖ dplyr::lag()     masks stats::lag()

4.1 ArcGIS REST API

The {arcgislayers} package allows users to read and write data from and to the ArcGIS REST API.

4.1.1 Reading Layers

sf_libraries_url <- "https://services.arcgis.com/Zs2aNLFN00jrS4gG/arcgis/rest/services/SF_Libraries/FeatureServer"
# arc_open can read a FeatureServer or a FeatureLayer directly
(sf_libraries_fs <- arc_open(sf_libraries_url))
<FeatureServer <2 layers, 0 tables>>
CRS: 3857
Capabilities: Query
  0: Libraries (esriGeometryPoint)
  1: Libraries with Air Conditioning (esriGeometryPoint)
(sf_libraries_lyr <- get_layer(sf_libraries_fs, name = "Libraries"))
<FeatureLayer>
Name: Libraries
Geometry Type: esriGeometryPoint
CRS: 3857
Capabilities: Query
sf_libraries <- arc_select(sf_libraries_lyr)
Registered S3 method overwritten by 'jsonify':
  method     from    
  print.json jsonlite
glimpse(sf_libraries)
Rows: 33
Columns: 23
$ objectid         <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16…
$ gross_sq_f       <chr> "8500", "6465", "6096", "8536", "7633", "6100", "8000…
$ block_lot        <chr> "8708110", "3564095", "6539034", "2919031", "0469001"…
$ zip_code         <chr> "94158", "94114", "94114", "94127", "94123", "94112",…
$ facility_i       <chr> "1113", "648", "1184", "1853", "1053", "896", "1858",…
$ city             <chr> "San Francisco", "San Francisco", "San Francisco", "S…
$ latitude         <chr> "37.775369728", "37.76406037", "37.750228042", "37.74…
$ department       <chr> "Public Library", "Public Library", "Public Library",…
$ longitude        <chr> "-122.393097384", "-122.431881717", "-122.435090242",…
$ dept_id          <chr> "48", "48", "48", "48", "48", "48", "48", "48", "48",…
$ common_nam       <chr> "Mission Bay Library", "Eureka Valley Branch Library/…
$ address          <chr> "960 04th St", "1 Jose Sarria Ct", "451 Jersey St", "…
$ supervisor       <chr> "6", "8", "8", "7", "2", "7", "5", "6", "10", "8", "9…
$ city_tenan       <chr> " ", " ", " ", " ", " ", " ", " ", " ", " ", " ", " "…
$ owned_leas       <chr> "Own", "Own", "Own", "Own", "Own", "Own", "Own", "Own…
$ globalid         <chr> "755632fc-18d9-4c0e-b65d-ec53d18a132b", "195a2a2b-302…
$ created_user     <chr> "nancy.milholland_sfdem", "nancy.milholland_sfdem", "…
$ created_date     <dttm> 2019-06-07 21:44:34, 2019-06-07 21:44:34, 2019-06-07…
$ last_edited_user <chr> "nancy.milholland_sfdem", "nancy.milholland_sfdem", "…
$ last_edited_date <dttm> 2019-06-07 21:44:34, 2019-06-07 21:44:34, 2019-06-07…
$ eas_id           <lgl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N…
$ air_conditioning <chr> "Yes", "No", "No", "No", "No", "No", "No", "No", "Sec…
$ geometry         <POINT [m]> POINT (-13624737 4547742), POINT (-13629055 454…
# You can also specify which columns to select, e.g.:
# arc_select(
#   sf_libraries_lyr,
#   fields = c("common_nam", "gross_sq_f", "address"),
#   where = "gross_sq_f < 8000"
# )

# With pipes:
# sf_libraries_url %>% 
#   arc_open() %>% 
#   get_layer(name = "Libraries") %>% 
#   arc_select()

# With pipes and tidyverse:
# (if the url points to a FeatureLayer instead of a FeatureServer)
# sf_libraries_url %>% 
#   arc_open() %>% 
#   select(common_nam, gross_sq_f, address) %>% 
#   filter(gross_sq_f < 8000) %>% 
#   collect()
mapview(sf_libraries, zcol = "air_conditioning")

It is also convenient to wrap this up into the body of a single function:

get_arcgis_layer <- function(lyr_name) {
  url <- glue::glue("https://services.arcgis.com/Zs2aNLFN00jrS4gG/arcgis/rest/services/{lyr_name}/FeatureServer/0")
  out <- arcgislayers::arc_select(arcgislayers::arc_open(url))
  return(out)
}

libraries <- get_arcgis_layer("SF_Libraries")

4.1.2 Writing Layers

If you have an sfgov.maps.arcgis.com account, you can write layers directly to your content. Read the authorization page for more information on credentials and tokens.

nc <- st_read(system.file("shape/nc.shp", package = "sf"))
tkn <- auth_code()
set_arc_token(tkn)

publish_res <- publish_layer(nc, "North Carolina SIDS sample")

4.2 ArcGIS Pro

If you have an ArcGIS Pro license, you can write directly to geodatabases within Projects using the {arcgisbinding} package.

Reading layer `OGRGeoJSON' from data source 
  `https://data.sfgov.org/api/geospatial/f2zs-jevy?accessType=DOWNLOAD&method=export&format=GeoJSON' 
  using driver `GeoJSON'
Simple feature collection with 11 features and 7 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -123.1738 ymin: 37.63983 xmax: -122.3279 ymax: 37.8632
Geodetic CRS:  WGS 84
library(arcgisbinding)
# arc.check_product()

# Get Supervisor Districts from DataSF:
sup_dists <- st_read("https://data.sfgov.org/api/geospatial/f2zs-jevy?accessType=DOWNLOAD&method=export&format=GeoJSON")

# Write to ArcGIS Pro project geodatabase
proj_path <- "...<full_path>.../ArcGIS/Projects/Test R Project/Test R Project.gdb/sup_dist"
arc.write(path = proj_path, data = sup_dists)

4.3 Spatial Joins

Use spatial joins to determine which points are ‘within’ which polygon:

nhoods <- st_read("https://data.sfgov.org/resource/j2bu-swwd.geojson")
Reading layer `j2bu-swwd' from data source 
  `https://data.sfgov.org/resource/j2bu-swwd.geojson' using driver `GeoJSON'
Simple feature collection with 41 features and 1 field
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -122.5149 ymin: 37.70813 xmax: -122.357 ymax: 37.8333
Geodetic CRS:  WGS 84
# The coordinate reference systems must match
st_crs(sf_libraries) == st_crs(sup_dists)
[1] FALSE
sf_libraries %>% 
  select(common_nam) %>% 
  st_transform(st_crs(sup_dists)) %>% 
  st_join(sup_dists %>% select(sup_dist), join = st_within) %>% 
  st_join(nhoods, join = st_within) %>% 
  st_drop_geometry() %>% 
  mutate(common_nam = gsub(" Branch| Library", "", common_nam)) %>% 
  head()
                           common_nam sup_dist               nhood
1                         Mission Bay        6         Mission Bay
2 Eureka Valley/ Harvey Milk Memorial        8 Castro/Upper Market
3                          Noe Valley        8          Noe Valley
4                         West Portal        7  West of Twin Peaks
5                              Marina        2              Marina
6                           Ingleside        7  West of Twin Peaks

4.4 Removing Farallon Islands from Supervisor Districts

d4 <- sup_dists %>% 
  filter(sup_dist == 4) %>% 
  st_cast("POLYGON") %>% 
  slice(1) %>% 
  st_cast("MULTIPOLYGON")
Warning in st_cast.sf(., "POLYGON"): repeating attributes for all
sub-geometries for which they may not be constant
sup_dists_no_farallon <- sup_dists %>% 
  filter(sup_dist != 4) %>% 
  bind_rows(d4)

mapview(sup_dists_no_farallon)

4.5 Census Data

The {tidycensus package} is fantastic, and the documentation is full of helpful examples.

library(tidycensus)

sf <- get_acs(
  state = "CA",
  county = "San Francisco",
  geography = "tract",
  variables = "B19013_001",
  geometry = TRUE,
  year = 2020
) %>% 
  st_transform(3857)
Getting data from the 2016-2020 5-year ACS
Downloading feature geometry from the Census website.  To cache shapefiles for use in future sessions, set `options(tigris_use_cache = TRUE)`.
sf_bbox <- libraries %>% 
  drop_na(city) %>% 
  st_buffer(3500) %>% 
  st_bbox()
  
sf %>%
  ggplot(aes(fill = estimate)) + 
  geom_sf(color = NA) + 
  labs(title = "Median Household Income, 2020", fill = NULL) +
  coord_sf(xlim = sf_bbox[c("xmin", "xmax")], ylim = sf_bbox[c("ymin", "ymax")], expand = TRUE) +
  scale_fill_viridis_c(option = "magma", labels = scales::dollar) +
  theme_bw()