In this lab, we will learn how to compute various geometric properties of spatial objects.

library(tidyverse)
library(macleish)
library(sf)

Goal: by the end of this lab, you will be able to use sf to perform computations on spatial objects.

dplyr operations

Most dplyr operations will just work on sf objects, since they are data.frames.

  1. Use rename() to change the variable called Sheet1__Ve in the forests layer to forest_type.
# SAMPLE SOLUTION
forests <- macleish_layers[["forests"]] %>%
  rename(forest_type = Sheet1__Ve)
  1. Use mutate() and st_area() to compute the area of each forest.
# SAMPLE SOLUTION
forests <- forests %>%
  mutate(area = st_area(geometry))

Spatial operations

We can perform geometric operations on spatial objects. For example, the st_intersection() function returns the geometry covered by both spatial objects. For example, in order to find the parts of the streams that are entirely inside the MacLeish boundary, we intersect those layers.

boundary <- macleish_layers %>%
  pluck("boundary")

streams <- macleish_layers %>%
  pluck("streams")

streams_inside <- boundary %>%
  st_intersection(streams)

library(leaflet)
macleish_map <- leaflet() %>%
  addTiles() %>%
  addPolygons(data = boundary, weight = 1)
macleish_map %>%
  addPolylines(data = streams_inside)

To find the corresponding part of the streams that lay outside of the MacLeish boundary, we use st_difference().

streams_outside <- streams %>%
  st_difference(boundary)

macleish_map %>%
  addPolylines(data = streams_outside)
  1. Find all segments of trails that extend beyond the MacLeish boundary
# SAMPLE SOLUTION
trails <- macleish_layers %>%
  pluck("trails")

trails_outside <- trails %>%
  st_difference(boundary)

macleish_map %>%
  addPolylines(data = trails_outside)

Spatial joins

Since all sf objects are also data.frames, you can add additional attributes to a spatial object using by joining that data.frame to your sf object in the usual way.

A distinct idea is to join two spatial layers together. This can be done with the st_join() function.

Just like regular joins, spatial joins can be inner joins or left joins (the default). However, spatial joins require an additional predicate that determines the type of spatial join. Common values include st_intersects(), st_contains(), st_touches(), etc. Please see help(st_intersects) for the full list.

For example, a way to find the complete streams that are entirely inside the MacLeish boundary is to join those two layers together. In this case we use the st_within predicate and set left to FALSE.

streams_inside <- streams %>%
  st_join(boundary, join = st_within, left = FALSE)

macleish_map %>%
  addPolylines(data = streams_inside)
  1. Find the stream that cross the boundary of MacLeish.
# SAMPLE SOLUTION
streams_xing <- streams %>%
  st_join(boundary, join = st_crosses, left = FALSE)

macleish_map %>%
  addPolylines(data = streams_xing)

Spatial aggregation

The dplyr operations group_by() and summarize() can also be used to aggregate spatial features. This may change the geometric nature of the features. The default method of aggregation is st_union(). Thus, for example, a collection of POINT features will become MULTIPOINT.

Suppose that we wanted to walk directly from the lowest point at MacLeish to the highest point. Both points are contained in the landmarks layer.

run_points <- macleish_layers %>%
  pluck("landmarks") %>%
  filter(str_detect(Label, "Point")) %>%
  summarize(N = n())

Note that since summarize() always condenses things down to one row, run_points has just one feature, and it is of type MULTIPOINT.

Unfortunately, leaflet canโ€™t plot MULTIPOINT geometries. We want to to show a line that connects the high point and the low point. So we need to convert our MULTIPOINT feature to a LINESTRING feature. We do this using the st_cast() function.

run_linestring <- run_points %>%
  st_cast("LINESTRING")

macleish_map %>%
  addPolylines(data = run_linestring)
  1. How far is it from the lowest point to the highest point?
run_linestring %>%
  st_length()
  1. Use group_by() and summarize() to compute the total area of each type of forest.
# SAMPLE SOLUTION
forests_total <- forests %>%
  group_by(forest_type) %>%
  summarize(N = n(), total_area = sum(area))
  1. Which type of forest takes occupies the largest area? Which type of forest is the most common?
# SAMPLE SOLUTION
forests_total %>%
  arrange(desc(total_area))
  1. Put only the Red Oak-Hemlock Forest on a leaflet map. This should be one MULTIPOLYGON feature.
red_oak <- forests_total %>%
  filter(forest_type == "Red Oak-Hemlock Forest")

# SAMPLE SOLUTION

macleish_map %>%
  addPolygons(
    data = red_oak,
    weight = 1,
    popup = ~forest_type
  )

Interpolation

sf contains a function called st_interpolate_aw() that will compute area-weighted interpolations of various quantities. Please explore and use at your own risk.

Your learning

Please respond to the following prompt on Slack in the #mod-spatial channel.

Prompt: What questions do you still have about working with spatial data in R? What features would you like to know more about?