Exercises: Spatial data wrangling with sf

Packages

library(tidyverse)
library(sf)
library(janitor)
library(tidygeocoder)

Data

Access

Access data

Create a folder data/.

NC PMTW Streams:

download.file(paste0("https://opendata.arcgis.com/datasets/", 
                    "0cc135e6c6244c9e9646b45ee3cb6c1e_0.zip", 
                    "?outSR=%7B%22latestWkid%22%3A3857%2C%22",
                    "wkid%22%3A102100%7D"),
              destfile = "data/pmtw.zip")

unzip("data/pmtw.zip", exdir = "data/")
file.remove("data/pmtw.zip")

NC Public Schools:

download.file(paste0("https://opendata.arcgis.com/datasets/", 
                    "dea6ff0e8b4743a0ba361e13a85a4c70_3.zip", 
                    "?outSR=%7B%22latestWkid%22%3A32119%2C%22",
                    "wkid%22%3A32119%7D"),
              destfile = "data/schools.zip")

unzip("data/schools.zip", exdir = "data/")
file.remove("data/schools.zip")

Read

Read in the shapefiles

nc_counties <- st_read(system.file("shape/nc.shp", package = "sf"), 
                       quiet = T) %>% 
  clean_names() %>% 
  select(name, area)

nc_pmtw <- st_read("data/PMTW_Streams_2020.shp", quiet = T) %>%
  clean_names() %>% 
  select(name = displ_name, designation = first_wrc, length = shape_len)

nc_schools <- st_read("data/Public_Schools.shp", quiet = T) %>% 
  clean_names()

Align CRS

nc_counties <- st_transform(nc_counties, crs = st_crs(nc_schools))
nc_pmtw <- st_transform(nc_pmtw, crs = st_crs(nc_schools))

st_crs(nc_counties)[1]

#> $input
#> [1] "NAD83 / North Carolina"

st_crs(nc_pmtw)[1]

#> $input
#> [1] "NAD83 / North Carolina"

st_crs(nc_schools)[1]

#> $input
#> [1] "NAD83 / North Carolina"

Exercise 1

Problem

Use st_intersects() and st_contains() to plot all the public schools in Durham and the counties that neighbor Durham.

Solution

durham <- nc_counties %>% 
  filter(str_detect(name, "^Durham"))

durham_counties <- durham %>% 
  st_intersects(nc_counties, sparse = FALSE) %>% 
  filter(.data = nc_counties, .)

durham_counties %>% 
  st_contains(nc_schools, sparse = FALSE) %>%
  apply(MARGIN = 2, FUN = any) %>% 
  filter(.data = nc_schools, .) %>% 
  ggplot() +
  geom_sf(color = "blue", alpha = .5, size = 2) +
  geom_sf(data = durham_counties, alpha = 0) +
  labs(title = "Durham Area Public Schools") +
  theme_void()

Exercise 2

Problem

Use st_buffer() to place a 500 meter buffer around the trout waters in Ashe county North Carolina.

Solution

ashe <- nc_counties %>% 
  filter(name == "Ashe")

ashe_trout <- ashe %>% 
  st_contains(nc_pmtw, sparse = FALSE) %>% 
  filter(.data = nc_pmtw, .)

ggplot(ashe) +
  geom_sf(alpha = 0.3) +
  geom_sf(data = st_buffer(ashe_trout, dist = 500), color = "pink", 
          alpha = 0.5) +
  geom_sf(data = ashe_trout, color = "blue") +
  labs(title = "Ashe county PMTW") +
  theme_void()

Exercise 3

Problem

Check the accuracy of the geocoding where method = census. Only use a sample of the schools. You can gauge the accuracy visually or by computing distances with st_distance(). Don’t forget to align the CRS!

Solution

Create a new sf object based on the geocode results. We’ll take a sample of 10 schools.

set.seed(3244)

nc_schools_sample <- nc_schools %>% 
  slice_sample(n = 10)

nc_schools_geocode <- nc_schools_sample %>% 
  st_drop_geometry() %>% 
  select(phys_addr, phys_city, phys_zip) %>% 
  mutate(state = "NC") %>% 
  geocode(street = phys_addr, city = phys_city, state = state, postalcode = phys_zip,
          method = "census") %>%
  filter(!is.na(lat)) %>% 
  st_as_sf(coords = c("long", "lat"), crs = "NAD27") %>% 
  st_transform(crs = st_crs(nc_schools))

Visualize how close these are.

ggplot(nc_schools_sample) + 
  geom_sf(color = "blue", alpha = .8, size = 5) + 
  geom_sf(data = nc_schools_geocode, color = "red", 
          shape = 21, size = 10) +
  geom_sf(data = nc_counties, alpha = 0) +
  theme_bw()