2 Set up the workspace

# Clean up the workspace
rm(list = ls())

# Load in the libraries
require(tidyverse)
require(sf)

# set the default ggplot theme to theme_bw() and center ggplot titles
theme_set(theme_bw())
theme_update(plot.title = element_text(hjust = 0.5))


# Read in the Lake Winnipeg Data using read_csv()
LWpg = readr::read_csv("LakeWpgOutline.csv") %>% 
  dplyr::select(Long, Lat)

# Read in islands CSVs
com = readr::read_csv("Commissioner Isle.csv") 
deer = readr::read_csv("Deer Isle.csv") 
reindeer = readr::read_csv("Reindeer Isle.csv") 
black = readr::read_csv("Black Isle.csv")

3 Spatial Data with R

3.1 Spatial Points

First, we’ll look at a plot of Lake Winnipeg without assigning a coordinate reference system.

# We convert the LWpg data frame into an sf object with st_as_sf()
LWpg_sf <- st_as_sf(LWpg, coords = c("Long", "Lat"))

# Note the additional geometry list-column which now holds the simple 
# feature collection with the coordinates of all the points.

# To make it a complete geographical object we assign the WGS84 projection,
# which has the EPSG code 4326:

ggplot() +
  geom_sf(data = LWpg_sf) +
  ggtitle("No coordinate reference system")

Now we’ll assign a CRS to the simple feature dataframe. Notice how the scaling automatically updates between the 2 plots.

st_crs(LWpg_sf) <- 4326 # we can use EPSG as numeric here
st_crs(LWpg_sf) # Coordinate system set
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
ggplot() +
  geom_sf(data = LWpg_sf) +
  ggtitle("CRS Assigned (EPSG: 4326)")

We’ve plotted our outline of Lake Winnipeg, but it’s plotted as points which is not ideal. Net we will convert the spatial points to a spatial polygon so that we plot the shoreline.

3.2 Spatial Polygon

Plotting a lake shoreline isn’t really useful as spatial points. We want to display the data as a spatial polygon. If you have a CSV of your data locations and you are trying to plot a polygon, the polygon has to be closed (the first point has to be the same as the last point). We can bind the first row to the end of the matrix to close the polygon. Not ideal, but a fairly easy fix.

If your data isn’t closed the st_polygon() function will throw an error.

# Run a test to convert the raw csv coords into a polygon
# xx = as.matrix(head(LWpg))
# st_polygon(list(xx))

# Need the polygon to be closed, so we bind the first data point to
# the end of the matrix to close the polygon. Not ideal, but fairly easy fix
xx = as.matrix(head(LWpg))
xx = rbind(xx, xx[1,])
st_polygon(list(xx))

rm(xx)

Now that we know how to convert raw data into a spatial polygon and plot it using geom_sf()

# Create simple feature geometry list column
LWpg_poly = st_sfc(st_polygon(list(as.matrix(rbind(LWpg, LWpg[1,])))))


# Set the coordinate system (same as what we did before)
st_crs(LWpg_poly)
## Coordinate Reference System: NA
st_crs(LWpg_poly) <- 4326 # we can use EPSG as numeric here
st_crs(LWpg_poly)
## Coordinate Reference System:
##   EPSG: 4326 
##   proj4string: "+proj=longlat +datum=WGS84 +no_defs"
# Not really neccesary, but could convert to sf as we did for the points
LWpg_poly = st_as_sf(LWpg_poly)


ggplot() +
  geom_sf(data = LWpg_poly, fill = "lightblue") +
  ggtitle("Map of Lake Winnipeg Boundary")

3.2.1 Add Point to Map

If you want to add points to the map, they don’t need to be a spatial object (ie. Spatial Points). You can just add coordinates and they will get added

# Easy to add points to polygon - they don't have to be spatial objects
ggplot() +
  geom_sf(data = LWpg_poly, fill = "lightblue") +
  ggtitle("Map of Lake Winnipeg Boundary") +
  geom_point(data = data.frame(x = c(-98.5, -96.75), y = c(53.5, 50.9)),
             aes(x, y))

4 Polygons and Hole

As far as I can tell, you don’t have to do anything special to plot holes (islands) in a polygon. If the polygons intersect the sf package seems to convert it for you. In other packages you might need to import the coords clockwise or counterclockwise to define a hole in a polygon

# convert the 2 dataframe to matrices
LWpg1 = as.matrix(rbind(LWpg, LWpg[1,]))
com1 = as.matrix(rbind(com, com[1,]))
deer1 = as.matrix(rbind(deer, deer[1,]))
reindeer1 = as.matrix(rbind(reindeer, reindeer[1,]))
black1 = as.matrix(rbind(black, black[1,]))

LWpg_poly = st_polygon(list(LWpg1, com1, deer1, reindeer1, black1))


# Create simple feature geometry list column
LWpg_poly = st_sfc(LWpg_poly)

st_crs(LWpg_poly) <- 4326 

ggplot() +
  geom_sf(data = LWpg_poly, fill = "lightblue") +
  ggtitle("Lake Winnipeg with Islands")

5 Combining sp and sf

5.1 Spatial Polygons: Part 1

This time, instead of having to close the polygons manually, use polygon functions in the sp package to create the polygons.

The first method creates a class “SpatialPolygons” and stores all the polygons together under a single id (1)

xy.sp <- sp::SpatialPolygons(list(
  sp::Polygons(list(sp::Polygon(LWpg),
                    sp::Polygon(com, hole = TRUE),
                    sp::Polygon(deer, hole = TRUE),
                    sp::Polygon(reindeer, hole = TRUE),
                    sp::Polygon(black, hole = TRUE)), "1")
))
# see str(xy.sp) for the interal structure

sf_poly <- as(xy.sp, "sf")
sf::st_crs(sf_poly) <- 4326

ggplot(sf_poly) +
  geom_sf(fill = "lightblue") +
  ggtitle("Lake Winnipeg")

5.2 Spatial Polygons: Part 2

The second method creates a class “SpatialPolygons” and stores all the polygons separately with IDs 1 to 5.

I prefer this second method because it allows more flexibility for plotting

xy.sp <- sp::SpatialPolygons(list(
  sp::Polygons(list(sp::Polygon(LWpg)),"1"),
  sp::Polygons(list(sp::Polygon(com, hole = TRUE)),"2"),
  sp::Polygons(list(sp::Polygon(deer, hole = TRUE)),"3"),
  sp::Polygons(list(sp::Polygon(reindeer, hole = TRUE)),"4"),
  sp::Polygons(list(sp::Polygon(black, hole = TRUE)),"5")
))
# see str(xy.sp) for the interal structure

sf_poly <- as(xy.sp, "sf")
sf::st_crs(sf_poly) <- 4326
sf_poly$id <- c(1, 2, 2, 2, 2)


ggplot(sf_poly) +
  geom_sf(aes(fill = factor(id))) +
  scale_fill_manual(values = c("1" = "lightblue", "2" = "burlywood")) +
  ggtitle("Lake Winnipeg") +
  theme(legend.position = "none")

6 Cropping

Say you only want to plot the south basin and part of the narrows. With geom_sf() you use coord_sf() - not coord_cartesian().

The beauty about using the sf package is that it creates plots that are to scale so you don’t have to manually size the outputs.

ggplot(sf_poly) +
  geom_sf(aes(fill = factor(id))) +
  scale_fill_manual(values = c("1" = "lightblue", "2" = "antiquewhite")) +
  theme(legend.position = "none") + 
  coord_sf(ylim = c(50.3, 51.5), xlim = c(-97,-96.25)) + 
  theme(panel.grid = element_blank(),
        panel.background = element_rect(fill = "antiquewhite")) + 
  scale_x_continuous(breaks = c(-96.25, -96.5, -96.75, -97))