In this lab, we will learn how to use the macleish package to retrieve weather and spatial data, and we will learn how to make interactive maps with leaflet.

library(tidyverse)

Goal: by the end of this lab, you will be able to place spatial geometries on an interactive map

The sf package

The sf package provides support for simple features in R. sf is a new, cutting-edge package developed by the r-spatial group group to be a complete, tidyverse-compatible suite of tools for working with geospatial data in R. It is under active development.

sf is designed to supersede the older spatial package sp. It is developed by the original developers of sp, along with Hadley Wickham and some other folks at RStudio.

In addition to supporting R packages, sf requires external dependencies. Similar to how you may have had to install git, you may have to install these dependencies at the operating system level. Installation instructions are here, and will vary across operating systems. In short, sf requires four external libraries:

  • GDAL: Geospatial Data Abstraction Library, which provides basic functionality for working with spatial data. GDAL is connected to R through the rgdal package.
  • GEOS provides geometric functions for computing area, length, intersections, etc. It is connected to R through the rgeos package.
  • Proj.4 is a library for computing projections. It is connected to R through several packages.
  • udunits, for converting one set of units to another.

In order to get sf to work, you need recent versions of all four of these libraries installed. It may take some trial-and-error to get all of these libraries installed properly. BE PATIENT!

You should install sf the way you would any other package, but pay careful attention to any error messages that crop up!

install.packages("sf")
  1. Install sf. Do not give up until it is successfully installed!

The macleish package

Much of what you need for #mp4 is in the macleish package for R. You can install this package directly from CRAN. You should have at least version 0.3.6.

# will install version 0.3.6
install.packages("macleish")
library(macleish)

The macleish package contains two main types of data: weather data from the two monitor stations, and spatial data in the form of shapefiles (actually sf data frames).

Spatial data

The spatial layers are stored in a list called macleish_layers. Each element in the list contains a different set of spatial data.

names(macleish_layers)

You can access a particular layer using the [[ operator or the pluck() function. The sf package does most of the heavy-lifting for spatial objects.

library(sf)
macleish_layers %>%
  pluck("buildings")

Note the coordinates, the projection string, and the features. This buildings layers contains POLYGON geometries that contain the information to draw the shapes of the buildings at the field station. These layers are of class sf and data.frame. This means that you can work with them just like data.frames, because they are data.frames!

map(macleish_layers, class)

You can plot the spatial data in a context-free fashion, but that is not really very informative, is it?

macleish_layers %>%
  pluck("buildings") %>%
  plot()

Some of the other layers have POINT geometries. Here we examine the landmarks at MacLeish.

macleish_layers %>%
  pluck("landmarks")

Please see the help for more information.

?macleish_layers

There is also a vignette that contains additional helpful context, examples, and exploratory data analysis. Read it!

vignette("macleish")
  1. What kind of geometry does the streams layer contain?
# SAMPLE SOLUTION
macleish_layers %>%
  pluck("streams")
  1. Plot the streams layer in a context-free fashion.
# SAMPLE SOLUTION
macleish_layers %>%
  pluck("streams") %>%
  plot()

Weather data

There are two data frames bundled with macleish that contain weather data from 2015: whately_2015 contains measurements from the Whately tower, and orchard_2015 contains measurements from the Orchard tower.

These measurements are averaged over every ten minute intervals and sent in real-time to a server in Bass Hall. You can can retrieve the full set of weather readings using the etl framework implemented by the macleish package.

mac <- etl("macleish") %>%
  etl_update()
whately <- mac %>%
  tbl("whately") %>%
  collect(n = Inf)

Please see the documentation for etl_extract.etl_macleish() for more information about how this process works.

  1. Create a data.frame that contains all of the weather readings from the orchard weather station. How many observations are there?
# SAMPLE SOLUTION
orchard <- mac %>%
  tbl("orchard") %>%
  collect(n = Inf)
nrow(orchard)

Static maps with ggplot2

Newer version of ggplot2 have built-in support for sf objects through the geom_sf() function. Since all sf objects include a geometry variable that contains the spatial information, geom_sf() automatically binds the geometry column to the x and y aesthetic. Other than that wrinkle, everything else should work as expected.

macleish_layers %>%
  pluck("buildings") %>%
  ggplot() +
  geom_sf()

Of course, you can layer additional features onto your maps by calling geom_sf() more than once.

ggplot() +
  geom_sf(data = pluck(macleish_layers, "buildings")) +
  geom_sf(data = pluck(macleish_layers, "streams"), color = "blue")
  1. Make a ggplot with the buildings, trails, landmarks, and streams layers all depicted (in different colors).
# SAMPLE SOLUTION
ggplot() +
  geom_sf(data = macleish_layers[["buildings"]]) +
  geom_sf(data = macleish_layers[["trails"]], color = "brown") +
  geom_sf(data = macleish_layers[["landmarks"]], color = "gold") +
  geom_sf(data = macleish_layers[["streams"]], color = "blue")

Interactive maps with leaflet

The Bechtel Environmental Classroom is located here:

bechtel <- tibble(lat = 42.449167, lon = -72.679389)

Let’s put it on a basic plot in leaflet. Note that here, bechtel is a regular data.frame – it does not contain any spatial geometries.

library(leaflet)

base_plot <- leaflet() %>%
  addTiles() %>%
  addMarkers(
    lng = ~lon, lat = ~lat, data = bechtel,
    popup = "Bechtel Environmental Classroom"
  )
base_plot

Note that the syntax here is similar to that of ggplot2, but with the pipe (%>%). Please consult the Leaflet documentation for more details.

You can add multiple layers to your leaflet plot.

base_plot %>%
  addPolygons(
    data = macleish_layers[["buildings"]],
    weight = 1, popup = ~name
  ) %>%
  addPolygons(
    data = macleish_layers[["forests"]],
    weight = 1, fillOpacity = 0.2, popup = ~type
  )
  1. Make an interactive visualization of the MacLeish property. Include layers for buildings, streams, and trails.
# SAMPLE SOLUTION
require(leaflet)
m <- leaflet() %>%
  addTiles() %>%
  addPolygons(data = macleish_layers[["buildings"]], weight = 1, popup = ~name) %>%
  addPolylines(data = macleish_layers[["trails"]], weight = 1, color = "brown") %>%
  addPolylines(data = macleish_layers[["streams"]], weight = 2)
m
# SAMPLE SOLUTION
forest_pal <- colorFactor("Greens", macleish_layers[["forests"]]$type)

m <- leaflet() %>%
  # Base groups
  addTiles(group = "OpenStreetMap") %>%
  addProviderTiles("Esri.WorldTopoMap", group = "Topography") %>%
  addProviderTiles("Esri.WorldImagery", group = "Satellite") %>%
  addProviderTiles("Stamen.TonerLite", group = "Toner Lite") %>%
  # Boundaries
  addPolygons(data = macleish_layers[["boundary"]], weight = 1, fillOpacity = 0.01, group = "Boundaries") %>%
  #  addPolygons(data = macleish_layers[["property_boundary"]], weight = 1, fill = "red", fillOpacity = 0.1) %>%
  # Man-Made structures
  addPolygons(
    data = macleish_layers[["buildings"]],
    weight = 1, popup = ~name, group = "Structures"
  ) %>%
  #  addPolygons(data = macleish_layers[["reservoir"]],
  #              weight = 1, group = "Structures") %>%
  addMarkers(
    data = macleish_layers[["landmarks"]],
    popup = ~Label, group = "Structures"
  ) %>%
  addPolylines(
    data = macleish_layers[["trails"]],
    weight = 1, color = "brown",
    popup = ~name, group = "Structures"
  ) %>%
  # Natural elements
  addPolygons(
    data = macleish_layers[["forests"]],
    color = ~ forest_pal(type), weight = 0.1,
    fillOpacity = 0.2,
    popup = ~type, group = "Natural"
  ) %>%
  addPolygons(
    data = macleish_layers[["wetlands"]],
    weight = 1, group = "Natural"
  ) %>%
  addPolylines(
    data = macleish_layers[["streams"]],
    weight = 2, group = "Natural"
  ) %>%
  # Layers control
  addLayersControl(
    baseGroups = c("OpenStreetMap", "Topography", "Satellite", "Toner Lite"),
    overlayGroups = c("Boundaries", "Structures", "Natural"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>%
  setView(lat = bechtel$lat, lng = bechtel$lon, zoom = 15)
m

Your learning

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

Prompt: Where did you get stuck in this lab? What errors are you seeing, and where are they appearing?