In thk686/keitt.ssi.2019: Learning materials for Geospatial Data Analysis in R, part of the UT Summer Statistics Institute

knitr::opts_chunk$set(warning = F, message = F)

# echo activity solutions?
solutions <- FALSE

Geospatial Data Analysis in R

Load some packages.

Note: You may need to set your working directory to the same as this file.

library(sf)
library(raster)
library(dplyr)
library(ggplot2)
library(RStoolbox) # provides fortify.raster

The st_join function

The sf library has an st_join function that performs a spatial join. For example, if you have a points dataset and you want to know which polygon in a polygons dataset goes with each point, then you use the st_join command.

busstops <- st_read(dsn = "ncshape", layer = "busstopsall", quiet = T)
censusblk <- st_read(dsn = "ncshape", layer = "censusblk_swwake", quiet = T)
join <- st_join(busstops, censusblk) # rows of censusblk matched to bus stops
head(join[, 1:7]) # 'cat' is the bus stop id

Activity: Hospital Elevation

Try joining the lakes layer (lhs) to the swwake census blocks layer (rhs). You will get a lot of NA values where there were not matches. You can filter these out.

Functions: st_read, st_join, is.na, filter, ggplot, geom_sf

censusblk <- st_read(dsn = "ncshape", layer = "censusblk_swwake", quiet = TRUE)
lakes <- st_read(dsn = "ncshape", layer = "lakes", quiet = TRUE)
join <- st_join(lakes, censusblk) %>% filter(!is.na(cat.y)) # remove non-matching rows
ggplot(join) +
  geom_sf(aes(fill = BLOCK_ID))

The extract function

The function extract in the raster package will choose pixel values in a raster based on a the points in an SpatialPoints object. The first argument of extract is a raster object. The second argument can be a wide variety of objects specifying points: 2-column matrix, 2-column data frame, a SpatialPoints object or an vector of indices. Methods are available for interpolating from neighboring cells and buffering the selection region. Here's an example. One niggling details is that you have to strip the Z-value from the coordinates of your SpatialPoints because extract expects only x, y values, not x, y, z.

zipc <- raster("ncrast/zipcodes.tif")
firest <- st_read(dsn = "ncshape", layer = "firestations", quiet = TRUE) %>% st_zm() # remove z-values
fszip <- extract(zipc, st_coordinates(firest)) # raster does not know about sf objects
res <- st_sf(LABEL = firest$LABEL, zipcode = fszip, geometry = firest$geometry) %>%
  filter(is.finite(zipcode))
print(res)
ggplot() +
  geom_sf(data = res, aes(color = as.factor(zipcode), size = 2), alapha = 0.1) +
  geom_sf(data = res, aes(size = 0.5)) +
  geom_raster(data = fortify(zipc), aes(x = x, y = y, fill = as.factor(zipcodes)), alpha = 0.5) +
  guides(fill = FALSE, color = FALSE, size = FALSE)

Activity: Hospital Elevation

Use the hospitals point layer and the elev_state_500m raster to extract the elevation of North Carolina hospitals.

Functions: raster, st_read, st_zm, extract, rgb, plot

elev.state <- raster("ncrast/elev_state_500m.tif")
hospitals <- st_read(dsn = "ncshape", layer = "hospitals", quiet = TRUE) %>% st_zm()
nc.state <- st_read(dsn = "ncshape", layer = "nc_state", quiet = TRUE) %>% st_zm()
hospitals$elevation <- extract(elev.state, st_coordinates(hospitals))
ggplot(nc.state) +
  geom_sf() +
  geom_raster(data = fortify(elev.state), aes(x = x, y = y, fill = elev_state_500m)) +
  scale_fill_gradientn(colours = terrain.colors(100, alpha = 0.2)) +
  scale_color_gradientn(colours = terrain.colors(100)) +
  geom_sf(data = hospitals, aes(color = elevation, size = 2), alpha = 0.5) +
  geom_sf(data = hospitals, aes(size = 0.5)) +
  guides(color = FALSE, fill = FALSE, size = FALSE) +
  theme_bw()

thk686/keitt.ssi.2019 documentation built on June 28, 2019, 1:28 p.m.