---
title: "Spatial data with sf"
subtitle: "Programming for Statistical Science"
author: "Shawn Santo"
institute: ""
date: ""
output:
  xaringan::moon_reader:
    css: "slides.css"
    lib_dir: libs
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: false
editor_options: 
  chunk_output_type: console
---

```{r include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE,
                      comment = "#>", highlight = TRUE,
                      fig.align = "center")

library(sf)
library(tidyverse)
library(patchwork)
library(mapview)
```

## Supplementary materials

Full video lecture available in Zoom Cloud Recordings

Additional resources

- Simple Features for R [vignettes](https://r-spatial.github.io/sf/)
- [CRS in R](https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf) 
  by Melanie Frazier
- [Leaflet for R](https://rstudio.github.io/leaflet/)

---

class: inverse, center, middle

# Introduction

---

## Spatial data is different

Our **typical tidy data frame**:

.tiny[
```{r echo=FALSE}
nycflights13::flights
```
]

---

## Spatial data is different

A **simple features object**:

.tiny[
```{r echo=FALSE}
nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = T)
nc %>% 
  select(AREA:NAME)
```
]

---

Another simple features object:

.tiny[
```{r echo=FALSE}
nc_gamelands <- read_sf("data/Game_Lands_-_general.shp")
nc_gamelands["GML_HAB"]
```
]

---

## Spatial data plotting needs care

```{r echo=FALSE, out.width="100%", fig.height=5}

lcc <- "+proj=lcc +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
nc_lcc <- nc %>% st_transform(crs = lcc)
aea <- "+proj=aea +lat_1=20 +lat_2=60 +lat_0=40 +lon_0=-96 +x_0=0 +y_0=0 +ellps=GRS80 +datum=NAD83 +units=m +no_defs"
nc_aea <- nc %>% st_transform(crs = aea)
rob <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
nc_rob <- nc %>% st_transform(crs = rob)

par(mar=c(3,2,2,1),mfrow=c(2,2))
plot(st_geometry(nc), graticule = T, axes = T)
plot(st_geometry(nc_lcc), graticule = T, axes = T)
plot(st_geometry(nc_aea), graticule = T, axes = T)
plot(st_geometry(nc_rob), graticule = T, axes = T)
```

---

```{r echo=FALSE, out.width="100%", fig.height=6}
nc <- st_transform(nc, st_crs(nc_gamelands))
durham <- nc[nc$NAME == "Durham", ]
durham <- nc[durham, ]
durham_gamelands <- nc_gamelands[durham, ]
plot(st_geometry(durham), axes = T, col = "#add8e6",
     main = "Durham Area Counties")
```

---

```{r echo=FALSE, out.width="100%", fig.height=6}
durham_gamelands_2 <- durham_gamelands %>% st_transform(crs = "+proj=utm +zone=17 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs")
plot(st_geometry(durham_gamelands_2), axes = T, col = "#ff6700",
     main = "Durham Area Public Game Lands")
```

---

class: center, middle

## Can we combine the two plots?

---

```{r echo=FALSE, out.width="100%", fig.height=5.5}
plot(st_geometry(durham), axes = T, main = "Durham Area Public Game Lands")
plot(st_geometry(durham_gamelands_2), add = T, axes = T, col = "#ff6700")
```

Where are the game lands?

---

class: center, middle

## We can, but more care is needed.

---

```{r echo=FALSE, out.width="100%", fig.height=5.5}
plot(st_geometry(durham), axes = T, main = "Durham Area Public Game Lands")
plot(st_geometry(durham_gamelands), add = T, axes = T, col = "#ff6700")
```

---

## Spatial data challenges

1. Different data types exist.

2. Special attention must be given to the coordinate reference system (CRS).

3. Manipulating spatial data objects is similar but not identical to 
   manipulating data frame objects.
   
---

class: inverse, center, middle

# Spatial data and R

---

## Analysis of spatial data in R

.pull-left[

<br/>

- Package `raster` contains classes and tools for handling spatial raster data.

<br/><br/>

- Package `sf` combines the functionality of `sp`, `rgdal`, and `rgeos` 
  into a single package based on tidy simple features.

]

.pull-right[

![](images/vector_raster_comparison.png)


]

<br/>

Whether or not you use vector or raster data depends on the type of problem 
and the data source. Our focus will be on vector data and package `sf`. 

*Source:* https://commons.wikimedia.org/wiki/File:Raster_vector_tikz.png

---

## Installing package `sf`

From https://r-spatial.github.io/sf/index.html

**Windows**

Installing `sf` from source works under windows when Rtools is installed. 
This downloads the system requirements from rwinlib.

**MacOS**

```bash
brew install pkg-config
brew install gdal
```
Once gdal is installed, you will be able to install sf package from source in R.

**Linux**

For Unix-alikes, GDAL (>= 2.0.1), GEOS (>= 3.4.0) and Proj.4 (>= 4.8.0) are 
required.

---

## Features and simple features

- A **feature** is a thing or object in the real world: a house, a city, a park,
  a forest, etc.
<br/><br/>
- A **simple feature** as defined by OpenGIS Abstract specification is to have 
  both spatial and non-spatial attributes. Spatial attributes are geometry 
  valued, and simple features are based on 2D geometry with linear 
  interpolation between vertices.
<br/><br/>

.tiny[
```{r eval=FALSE}
Simple feature collection with 100 features and 1 field
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 123829.8 ymin: 14740.06 xmax: 930518.6 ymax: 318255.5
projected CRS:  NAD83 / North Carolina
First 10 features:
          NAME                       geometry
1         Ashe MULTIPOLYGON (((387344.7 27... #<<
2    Alleghany MULTIPOLYGON (((408601.4 29...
3        Surry MULTIPOLYGON (((478715.7 27...
4    Currituck MULTIPOLYGON (((878193.4 28...
5  Northampton MULTIPOLYGON (((769834.9 27...
6     Hertford MULTIPOLYGON (((812327.7 27...
7       Camden MULTIPOLYGON (((878193.4 28...
8        Gates MULTIPOLYGON (((828444.5 29...
9       Warren MULTIPOLYGON (((671746.3 27...
10      Stokes MULTIPOLYGON (((517435.1 27...
```
]

---

## Simple features examples

```{r echo=FALSE}
pt <- st_point(c(30, 10))
ls <- st_linestring(matrix(c(30, 10, 10, 30, 40, 40), byrow = TRUE, ncol =
                            2))
poly <- st_polygon(list(matrix(
  c(30, 10, 40, 40, 20, 40, 10, 20, 30, 10),
  ncol = 2,
  byrow = TRUE
)))
poly_hole <- st_polygon(list(matrix(
  c(35, 10, 45, 45, 15, 40, 10, 20, 35, 10),
  ncol = 2,
  byrow = TRUE
),
matrix(
  c(20, 30, 35, 35, 30, 20, 20, 30),
  ncol = 2,
  byrow = TRUE
)))

pts <- st_multipoint(matrix(
  c(10, 40, 40, 30, 20, 20, 30, 10),
  ncol = 2,
  byrow = TRUE
))
lss <- st_multilinestring(list(matrix(
  c(10, 10, 20, 20, 10, 40), ncol = 2, byrow = TRUE
),
matrix(
  c(40, 40, 30, 30, 40, 20, 30, 10),
  ncol = 2,
  byrow = TRUE
)))
polys <- st_multipolygon(list(list(matrix(
  c(30, 20, 45, 40, 10, 40, 30, 20),
  ncol = 2,
  byrow = TRUE
)),
list(matrix(
  c(15, 5, 40, 10, 10, 20, 5, 10, 15, 5),
  ncol = 2,
  byrow = TRUE
))))
polys_hole <- st_multipolygon(list(list(matrix(
  c(40, 40, 20, 45, 45, 30, 40, 40),
  ncol = 2,
  byrow = TRUE
)),
list(
  matrix(
    c(20, 35, 10, 30, 10, 10, 30, 5, 45, 20, 20, 35),
    ncol = 2,
    byrow = TRUE
  ),
  matrix(
    c(30, 20, 20, 15, 20, 25, 30, 20),
    ncol = 2,
    byrow = TRUE
  )
)))
```

```{r echo=FALSE, fig.width=8, fig.height=8, fig.align="center", out.width="100%"}
par(mar = c(1, 1, 2, 1), mfrow = c(4, 4))

plot(pt,
     axes = FALSE,
     main = "Point",
     pch = 16)
box()
plot(ls, axes = FALSE, main = "Linestring")
box()
plot(poly,
     axes = FALSE,
     col = "lightgrey",
     main = "Polygon")
box()
plot(poly_hole,
     axes = FALSE,
     col = "lightgrey",
     main = "Polygon w/ Hole(s)")
box()

plot(pts,
     axes = FALSE,
     main = "Multipoint",
     pch = 16)
box()
plot(lss, axes = FALSE, main = "Multilinestring")
box()
plot(polys,
     axes = FALSE,
     col = "lightgrey",
     main = "Multipolygon")
box()
plot(polys_hole,
     axes = FALSE,
     col = "lightgrey",
     main = "Multipolygon w/ Hole(s)")
box()
```

---

## `sf` objects

.tiny[
```{r}
nc <- st_read(system.file("shape/nc.shp", package = "sf"), quiet = TRUE)
nc
```
]

---

## Class and other attributes: `sf`

```{r}
class(nc)
```

```{r}
names(attributes(nc))
```


---

## `sfc` objects

```{r}
nc_polygons <- st_geometry(nc)
nc_polygons
```

---

## Class and other attributes: `sfc`

```{r}
class(nc_polygons)
```

```{r}
names(attributes(nc_polygons))
```

<br/>

--

We see that `nc` has a class attribute `sf`, and object `nc_polygons` has a
class attribute `sfc`. What methods are available?

---

```{r}
methods(class = "sf")
```

---

```{r}
methods(class = "sfc")
```

---

## Reading and writing spatial data

- `st_read()` / `st_write()`, Shapefile, GeoJSON, KML, ...

- `st_as_sfc()` 

- `st_as_text()`,  well-known text format

- `st_as_binary()`, well-known binary format

<br/><br/><br>

See https://r-spatial.github.io/sf/articles/sf2.html for the full
set of driver availability.

---

## Plotting with `plot()`

```{r fig.height=6}
plot(nc)
```

---

```{r fig.width=7, fig.height=6}
plot(nc["NAME"])
```

---

```{r fig.width=12, fig.height=6}
par(oma=c(0,2,0,0))
plot(nc["AREA"], graticule = TRUE, axes = TRUE, las = 1)
```

---

## What is happening with `[` and the `sf` object?

```{r}
nc["AREA"]
```

---

```{r}
nc$AREA
```

---

```{r fig.width=12, fig.height=6}
par(oma=c(0,2,0,0))
plot(nc["AREA"], col = "lightblue", graticule = TRUE, 
     axes = TRUE, las = 1)
```

---

```{r fig.width=12, fig.height=6}
par(oma=c(0,2,0,0))
plot(st_geometry(nc), graticule = TRUE, axes = TRUE, las = 1)
```

---

## Plotting with `ggplot()`

```{r fig.width=9, fig.height=6}
ggplot(nc) +
  geom_sf() +
  theme_bw(base_size = 16)
```

---

```{r fig.width=12, fig.height=8}
ggplot(nc) +
  geom_sf(aes(fill = AREA)) +
  scale_fill_gradient(low = "#fee8c8", high = "#7f0000") +
  theme_bw(base_size = 16)
```

---

class: center, middle

```{r echo=FALSE, fig.width=12, fig.height=8}
p1 <- ggplot(nc) +
  geom_sf(aes(fill = SID74)) +
  scale_fill_gradient(low = "#fff7f3", high = "#49006a") +
  theme_bw(base_size = 16)
  
p2 <- ggplot(nc) +
  geom_sf(aes(fill = SID79)) +
  scale_fill_gradient(low = "#fff7f3", high = "#49006a") +
  theme_bw(base_size = 16)

p1 / p2
```

---

```{r eval=FALSE, fig.width=12, fig.height=8}
p1 <- ggplot(nc) +
  geom_sf(aes(fill = SID74)) +
  scale_fill_gradient(low = "#fff7f3", high = "#49006a") +
  theme_bw(base_size = 16)
  
p2 <- ggplot(nc) +
  geom_sf(aes(fill = SID79)) +
  scale_fill_gradient(low = "#fff7f3", high = "#49006a") +
  theme_bw(base_size = 16)

p1 / p2
```

<br/><br/>

--

Visually, what is wrong with the last plot?

---

## Plotting with `mapview()`

```{r eval=FALSE}
mapview(nc)
```

--

```{r eval=FALSE}
mapviewOptions(legend.pos = "bottomright")
mapview(nc["SID74"], col.regions = sf.colors(10), 
        layer.name = "SID 1974")
```

<br/><br/>

*These should run in RStudio. There is an issue embedding this overlay in the slides*.

---

## Exercise

Use `ggplot` to create a choropleth map for the proportion of sudden infant 
deaths, for the period of July 1, 1974 to June 30, 1979.

```{r echo=FALSE, out.width="100%", fig.height=5}
nc %>% 
  select(BIR74, SID74) %>% 
  mutate(SID74_prop = SID74 / (BIR74 + SID74)) %>% 
  st_as_sf() %>% 
  ggplot() +
  geom_sf(aes(fill = SID74_prop)) +
  scale_fill_gradient(low = "#fff7f3", high = "#49006a") +
  labs(title = "July 1, 1974 to June 30, 1979", 
       fill  = "",
       subtitle = "Proportion of SID by county") +
  theme_bw(base_size = 16)
```

???

```{r eval=FALSE}
nc %>% 
  select(BIR74, SID74) %>% 
  mutate(SID74_prop = SID74 / (BIR74 + SID74)) %>% 
  st_as_sf() %>% 
  ggplot() +
  geom_sf(aes(fill = SID74_prop)) +
  scale_fill_gradient(low = "#fff7f3", high = "#49006a") +
  labs(title = "July 1, 1974 to June 30, 1979", 
       fill  = "",
       subtitle = "Proportion of SID by county") +
  theme_bw(base_size = 16)
```

---

class: inverse, center, middle

# Map layers

---

## Game Lands data

The North Carolina Department of Environment and Natural Resources, 
Wildlife Resources Commission and the NC Center for Geographic Information and 
Analysis has a shapefile data set available on all public Game Lands in NC.

https://www.nconemap.gov/datasets/e5ddff9b96204c6181be7c022e61d946_0

We can directly download and unzip the shapefile via

```{r eval=FALSE}
download.file("https://opendata.arcgis.com/datasets/e5ddff9b96204c6181be7c022e61d946_0.zip?outSR=%7B%22latestWkid%22%3A32119%2C%22wkid%22%3A32119%7D",
              destfile = "data/Gamelands.zip")
unzip("data/Gamelands.zip", exdir = "data/")
```

--

To see the available files

```{r}
list.files(path = "data/", pattern = "Game_Lands*")
```

---

## Read in the shapefile

```{r}
nc_gamelands <- st_read("data/Game_Lands_-_general.shp", quiet = TRUE)
```

.tiny[
```{r}
print(nc_gamelands, n = 5)
```
]

---

## Metadata for each `sf` object

`nc`:

```{r eval=FALSE}
Simple feature collection with 100 features and 14 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
geographic CRS: NAD27 #<<
```

`nc_gamelands`:

```{r eval=FALSE}
Simple feature collection with 94 features and 6 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: 127456.7 ymin: 26544.91 xmax: 923528.7 ymax: 318097.4
projected CRS:  NAD83 / North Carolina #<<
```

---

## Check the CRS

```{r eval=FALSE}
st_crs(nc)
```

```{r eval=FALSE}
Coordinate Reference System:
  User input: NAD27 
  wkt:
GEOGCRS["NAD27",
    DATUM["North American Datum 1927",
        ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["latitude",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["longitude",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4267]]
```


---

.tiny[
```{r eval=FALSE}
st_crs(nc_gamelands)
```

```{r eval=FALSE}
Coordinate Reference System:
  User input: NAD83 / North Carolina 
  wkt:
PROJCRS["NAD83 / North Carolina",
    BASEGEOGCRS["NAD83",
        DATUM["North American Datum 1983",
            ELLIPSOID["GRS 1980",6378137,298.257222101,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4269]],
    CONVERSION["SPCS83 North Carolina zone (meters)",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",33.75,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
                
                  ⋮    
          
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting (X)",east,
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["northing (Y)",north,
            ORDER[2],
            LENGTHUNIT["metre",1]],
    USAGE[
        SCOPE["unknown"],
        AREA["USA - North Carolina"],
        BBOX[33.83,-84.33,36.59,-75.38]],
    ID["EPSG",32119]]
```
]

---

## Coordinate reference systems (CRS)

- CRS provide a standardized way of describing locations.

- Different CRS arise from various ways data were gathered, the locations,
  and purposes of the data.
  
- A CRS is comprised of
    - an ellipsoid, to define the earth's shape;
    - a datum, to define the origin and orientation of coordinate axes;
    - a projection, to go from 3D to 2D.

- It is important that you transform your spatial data to a common CRS before
  plotting.

---

## Transform CRS

```{r}
nc_gamelands <- st_transform(nc_gamelands, crs = st_crs(nc))
```

--

Check they are equal:

```{r}
st_crs(nc) == st_crs(nc_gamelands)
```

---

## Map overlay with `plot()`

```{r fig.show='hide'}
plot(st_geometry(nc), axes = T, las = 1, main = "NC Public Game Lands",
     cex.main = 3, cex.lab = 2, cex.axis = 1.5)

plot(st_geometry(nc_gamelands), 
     add = T, #<<
     col = "#ff6700")

legend("bottomleft", legend = "WRC Game Lands", fill = "#ff6700")
```

---

## Map overlay with `plot()`

```{r echo=FALSE, fig.width=12, fig.height=8}
plot(st_geometry(nc), axes = T, las = 1, main = "NC Public Game Lands",
     cex.main = 3, cex.lab = 2, cex.axis = 1.5)
plot(st_geometry(nc_gamelands), add = T, col = "#ff6700")
legend("bottomleft", legend = "WRC Game Lands", fill = "#ff6700")
```

---

## Map overlay with `mapview()`

```{r eval=FALSE}
nc_mapview <- mapview(nc, alpha.regions = .2, alpha = .9,
                  label = nc[, "NAME", drop = T],
                  layer.name = "NC Counties")
```

```{r eval=FALSE}
nc_gamelands_mapview <- mapview(nc_gamelands, col.regions = "#ff6700",
        label = round(nc_gamelands[, "SUM_ACRES", drop = T], 2),
        layer.name = "NC Gamelands")
```

```{r eval=FALSE}
nc_mapview + nc_gamelands_mapview
```

<br/><br/><br/>

*These should run in RStudio. There is an issue embedding this overlay in the slides*.

---

## Exercise

Create a map that includes NC county boundaries, Game Lands, and hazardous waste
sites. Data for the hazardous waste sites is available at
https://www.nconemap.gov/datasets/hazardous-waste-sites

This data set represents the location of sites within North Carolina that 
are regulated by the hazardous waste portions of the 
Resource Conservation and Recovery Act (RCRA).

---

class: inverse, center, middle

# Manipulating `sf` type objects

---

## Change the CRS

We'll make a quick change to the CRS to better manipulate the geometries.

```{r}
nc <- st_transform(nc, st_crs(32119))
nc_gamelands <- st_transform(nc_gamelands, st_crs(32119))
```

*Source*: https://spatialreference.org/ref/epsg/32119/

<br/>

--

To make it easier to view the tibbles, we'll drop some of the fields.

```{r}
nc <- nc %>% 
  select(NAME)

nc_gamelands <- nc_gamelands %>% 
  select(GML_HAB)
```


---

## Intersects

```{r}
durham_county <- nc %>% 
  filter(NAME == "Durham")
durham_county
```

---

```{r}
nc[durham_county, ]
```

<br/>

--

**What is happening here? How can we verify this in the help?**

---

```{r}
st_intersects(nc, durham_county, sparse = F) %>% 
  nc[., ]
```

<br/>

--

Intersects finds if `nc` and `durham_county` geometries share any space.

---

## Touches

```{r}
st_touches(nc, durham_county, sparse = F) %>% 
  nc[., ]
```

<br/>

--

Touches identifies if `nc` and `durham_county` geometries share a common
point but their interiors do not intersect.

---

## Join

Suppose we want to plot all the game lands that intersect with Durham county or
one of its neighboring counties.

```{r}
durham_area_counties <- st_intersects(nc, durham_county, sparse = F) %>% 
  nc[., ]

durham_area_gamelands <- st_join(nc_gamelands, durham_area_counties,
                                 left = FALSE,
                                 join = st_intersects)
```

--

.pull-left[
```{r eval=FALSE}
nc_gamelands
```

```{r eval=FALSE}
          GML_HAB                       geometry
1           Alcoa MULTIPOLYGON (((512096.2 18...
2 Alligator River MULTIPOLYGON (((869633.1 24...
3      Angola Bay MULTIPOLYGON (((713079.4 11...
                      ⋮            
92      White Oak River MULTIPOLYGON (((779438 1142...
93 Whitehall Plantation MULTIPOLYGON (((669904.7 77...
94    William H. Silver MULTIPOLYGON (((233548.1 20...                                  
```
]

.pull-right[
```{r eval=FALSE}
durham_area_counties
```

```{r eval=FALSE}
        NAME                       geometry
13 Granville MULTIPOLYGON (((632225.8 25...
14    Person MULTIPOLYGON (((626993.2 27...
29    Orange MULTIPOLYGON (((607985.9 23...
30    Durham MULTIPOLYGON (((607985.9 23...
37      Wake MULTIPOLYGON (((616777.2 20...
48   Chatham MULTIPOLYGON (((559249.4 19...
```
]

---

```{r}
durham_area_gamelands
```

---

## Proximity

Suppose we want to find all the counties within 17,550 meters of Durham county.

```{r}
st_is_within_distance(durham_county, nc, dist = 17550, sparse = F) %>% 
  nc[., ]
```

---

## Exercise

Create a plot of North Carolina's Game Lands and all the waste sites within
100 meters of a Game Land area.

???

## Solution

.solution[
```{r eval=FALSE}
waste <- st_read("data/Hazardous_Waste_Sites.shp", quiet = T)
waste <- st_transform(waste, st_crs(nc_gamelands))
```

```{r eval=FALSE}
close_waste_lgl <- st_is_within_distance(waste, nc_gamelands, 
                                         dist = 100, sparse = F)
close_waste <- waste %>% 
  filter(apply(close_waste_lgl, 1, sum) > 0)
```

```{r eval=FALSE}
nc_gamelands_mapview <- mapview(nc_gamelands, col.regions = "#ff6700",
        label = nc_gamelands$GML_HAB,
        layer.name = "NC Gamelands")

nc_waste_mapview <- mapview(close_waste, col.regions = "#65ff00",
                            alpha = .3,
                            alpha.regions = .3,
                            label = waste[, "SITE_NAME", drop = T],
                            layer.name = "NC Waste Sites")

nc_gamelands_mapview + nc_waste_mapview
```

```{r eval=FALSE}
ggplot(nc) +
  geom_sf(alpha = .3) +
  geom_sf(data = close_waste, color = "#65ff00", size = 3) +
  geom_sf(data = nc_gamelands, fill = "#ff6700", alpha = .5) +
  theme_bw()
```

]

---

## References

1. Interactive Viewing of Spatial Data in R. (2020).         
   https://r-spatial.github.io/mapview/index.html.
   
2. Melanie Frazier. Coordinate Reference Systems in R. 
   https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf.
   
3. Simple Features for R. (2020). 
   https://r-spatial.github.io/sf/.

