spatial.Rmd

```{r include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
```

# (PART) Spatial {-}

# Spatial Data and Maps

We'll explore the basics of simple features (sf) for building spatial datasets, then some common mapping methods:

- ggplot2
- tmap
- leaflet
- the base plot system occasionally
- something else

## Spatial Data 

To work with spatial data requires extending R to deal with it using packages.  Many have been developed, but the field is starting to mature using international open GIS standards.

**`sp`**  (until recently, the dominant library of spatial tools)

- Includes functions for working with spatial data
- Includes `spplot` to create maps
- Also needs `rgdal` package for `readOGR` – reads spatial data frames.  

**`sf`** (Simple Features)

- ISO 19125 standard for GIS geometries
- Also has functions for working with spatial data, but clearer to use.
- Doesn't need many additional packages, though you may still need `rgdal` installed for some tools you want to use.
- Replacing `sp` and `spplot` though you'll still find them in code. We'll give it a try...
- Works with ggplot2 and tmap for nice looking maps.

Cheat sheet: https://github.com/rstudio/cheatsheets/raw/master/sf.pdf

#### simple feature geometry sfg and simple feature column sfc



### Examples of simple geometry building in sf 

sf functions have the pattern st_* 

st means "space and time"

See Geocomputation with R at https://geocompr.robinlovelace.net/ or  https://r-spatial.github.io/sf/
	for more details, but here's an example of manual feature creation of sf geometries (sfg):

```{r message=FALSE}
library(tidyverse)
library(sf)
library(iGIScData)
```

[As usual, go to the relevant project, in this case **`generic_methods`**]
```{r fig.cap="Building simple geometries in sf"}
library(sf)
eyes <- st_multipoint(rbind(c(1,5), c(3,5)))
nose <- st_point(c(2,4))
mouth <- st_linestring(rbind(c(1,3),c(3, 3)))
border <- st_polygon(list(rbind(c(0,5), c(1,2), c(2,1), c(3,2), 
                              c(4,5), c(3,7), c(1,7), c(0,5))))
face <- st_sfc(eyes, nose, mouth, border)  # sfc = sf column 
plot(face)
```

The face was a simple feature column (sfc) built from the list of sfgs. 
An sfc just has the one column, so is not quite like a shapefile.

- But it can have a coordinate referencing system CRS, and so can be mapped.
- Kind of like a shapefile with no other attributes than shape

[**`westUS`**]

### Building a mappable sfc from scratch

```{r fig.cap="A simple map built from scratch with hard-coded data as simple feature columns"}
CA_matrix <- rbind(c(-124,42),c(-120,42),c(-120,39),c(-114.5,35),
  c(-114.1,34.3),c(-114.6,32.7),c(-117,32.5),c(-118.5,34),c(-120.5,34.5),
  c(-122,36.5),c(-121.8,36.8),c(-122,37),c(-122.4,37.3),c(-122.5,37.8),
  c(-123,38),c(-123.7,39),c(-124,40),c(-124.4,40.5),c(-124,41),c(-124,42))
NV_matrix <- rbind(c(-120,42),c(-114,42),c(-114,36),c(-114.5,36),
  c(-114.5,35),c(-120,39),c(-120,42))
CA_list <- list(CA_matrix);       NV_list <- list(NV_matrix)
CA_poly <- st_polygon(CA_list);   NV_poly <- st_polygon(NV_list)
sfc_2states <- st_sfc(CA_poly,NV_poly,crs=4326)  # crs=4326 specifies GCS
st_geometry_type(sfc_2states)
library(tidyverse)
ggplot() + geom_sf(data = sfc_2states)

```

**sf class**

Is like a shapefile:  has attributes to which geometry is added, and can be used like a data frame.

```{r fig.cap="Using an sf class to build a map, displaying an attribute"}
attributes <- bind_rows(c(abb="CA", area=423970, pop=39.56e6),
                        c(abb="NV", area=286382, pop=3.03e6))
twostates <- st_sf(attributes, geometry = sfc_2states)
ggplot(twostates) + geom_sf() + geom_sf_text(aes(label = abb))
```

### Creating features from shapefiles or tables

**sf's `st_read` reads shapefiles**

- shapefile is an open GIS format for points, polylines, polygons

You would normally have shapefiles (and all the files that go with them -- .shx, etc.)
stored on your computer, but we'll access one from the iGIScData external data folder [**`sierra`**]:

```{r}
library(iGIScData)
library(sf)
shpPath <- system.file("extdata","CA_counties.shp", package="iGIScData")
CA_counties <- st_read(shpPath)
plot(CA_counties)
```

**`st_as_sf` converts data frames**

- using coordinates read from x and y variables, with crs set to coordinate system (4326 for GCS)

```{r}
sierraFebpts <- st_as_sf(sierraFeb, coords = c("LONGITUDE", "LATITUDE"), crs=4326)
plot(sierraFebpts)
```

[**`air_quality`**]
```{r, message=FALSE, warning=FALSE, fig.cap="ggplot map of Bay Area TRI sites, census centroids, freeways"}
library(tidyverse)
library(sf)
library(iGIScData)
censusCentroids <- st_centroid(BayAreaTracts)
TRI_sp <- st_as_sf(TRI_2017_CA, coords = c("LONGITUDE", "LATITUDE"), 
        crs=4326) # simple way to specify coordinate reference
bnd <- st_bbox(censusCentroids)
ggplot() +
  geom_sf(data = BayAreaCounties, aes(fill = NAME)) +
  geom_sf(data = censusCentroids) +
  geom_sf(data = CAfreeways, color = "grey") +
  geom_sf(data = TRI_sp, color = "yellow") +
  coord_sf(xlim = c(bnd[1], bnd[3]), ylim = c(bnd[2], bnd[4])) +
  labs(title="Bay Area Counties, Freeways and Census Tract Centroids")
```

### Coordinate Referencing System

Say you have data you need to make spatial with a spatial reference

`sierra <- read_csv("sierraClimate.csv")`

EPSG or CRS codes are an easy way to provide coordinate referencing.  

Two ways of doing the same thing. 

1. Spell it out:
```
GCS <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
wsta = st_as_sf(sierra, coords = c("LONGITUDE","LATITUDE"), crs=GCS)
```

2. Google to find the code you need and assign it to the crs parameter:

`wsta <- st_as_sf(sierra, coords = c("LONGITUDE","LATITUDE"), crs=4326)`

#### *Removing* Geometry

There are many instances where you want to remove geometry from a sf data frame

- Some R functions run into problems with geometry and produce confusing error messages, like "non-numeric argument"

- You're wanting to work with an sf data frame in a non-spatial way

One way to remove geometry:

`myNonSFdf <- mySFdf %>% st_set_geometry(NULL)`

### Spatial join `st_join`

A spatial join with st_join
joins data from census where TRI points occur [**`air_quality`**]

```{r, message=FALSE}
TRI_sp <- st_as_sf(TRI_2017_CA, coords = c("LONGITUDE", "LATITUDE"), crs=4326) %>%
  st_join(BayAreaTracts) %>%
  filter(CNTY_FIPS %in% c("013", "095"))
```


### Plotting maps in the base plot system

There are various programs for creating maps from spatial data, and we'll look at a few after we've looked at rasters. As usual, the base plot system often does something useful when you give it data.

```{r}
plot(BayAreaCounties)
```

And with just one variable:

```{r}
plot(BayAreaCounties["POP_SQMI"])
```

There's a lot more we could do with the base plot system, but we'll mostly focus on
some better options in ggplot2 and tmap.


## Raster GIS in R

Simple *Features* are feature-based, of course, so it's not surprising that sf doesn't have support for rasters. So we'll want to use the raster package.

We can start by building one from scratch:

```{r}
library(raster)
new_ras <- raster(nrows = 10, ncols = 10, 
                      xmn = 0, xmx = 100, ymn = 0, ymx = 100,
                      vals = 1:100)
plot(new_ras)

```



**A bit of raster reading and map algebra with Marble Mountains elevation data [`marbles`]**

```{r message=FALSE}
library(raster)
rasPath <- system.file("extdata","elev.tif", package="iGIScData")
elev <- raster(rasPath)
slope <- terrain(elev, opt="slope")
aspect <- terrain(elev, opt="aspect")
slopeclasses <-matrix(c(0,0.2,1, 0.2,0.4,2, 0.4,0.6,3,
                        0.6,0.8,4, 0.8,1,5), ncol=3, byrow=TRUE)
slopeclass <- reclassify(slope, rcl = slopeclasses)

plot(elev)
plot(slope)
plot(slopeclass)
plot(aspect)
```

**Sinking Cove, Tennessee** is a karst valley system carved into the Cumberland Plateau, a nice place to see the use of a hillshade raster created from a digital elevation model using raster functions for slope, aspect, and hillshade:  

```{r}
library(sf); library(tidyverse); library(tmap)
library(raster)
tmap_mode("plot")
DEMpath <- system.file("extdata/SinkingCove","DEM_SinkingCoveUTM.tif",package="iGIScData")
DEM <- raster(DEMpath)
slope <- terrain(DEM, opt='slope')
aspect <- terrain(DEM, opt='aspect')
hillsh <- hillShade(slope, aspect, 40, 330)
#
# Need to crop a bit since grid north != true north 
bbox0 <- st_bbox(DEM)
xrange <- bbox0$xmax - bbox0$xmin
yrange <- bbox0$ymax - bbox0$ymin
bbox1 <- bbox0
crop <- 0.05
bbox1[1] <- bbox0[1] + crop * xrange # xmin
bbox1[3] <- bbox0[3] - crop * xrange # xmax
bbox1[2] <- bbox0[2] + crop * yrange # ymin
bbox1[4] <- bbox0[4] - crop * yrange # ymax
bboxPoly <- bbox1 %>% st_as_sfc() # makes a polygon
#
tm_shape(hillsh, bbox=bboxPoly) +
  tm_raster(palette="-Greys",legend.show=F,n=20) +
  tm_shape(DEM) +
  tm_raster(palette=terrain.colors(24), alpha=0.5) +
  tm_graticules(lines=F)
```

*See ?raster to learn more about the rich array of raster GIS operations.*


## ggplot2 for maps

The Grammar of Graphics is the gg of ggplot.

- Key concept is separating aesthetics from data
- Aesthetics can come from variables (using aes()setting) or be constant for the graph

Mapping tools that follow this lead

- ggplot, as we have seen, and it continues to be enhanced
- tmap (Thematic Maps) https://github.com/mtennekes/tmap
Tennekes, M., 2018, tmap: Thematic Maps in R, *Journal of Statistical Software* 84(6), 1-39

```{r}
ggplot(CA_counties) + geom_sf()

```


Try `?geom_sf` and you'll find that its first parameters is mapping with `aes()` by default. The data property is inherited from the ggplot call, but commonly you'll want to specify data=something in your geom_sf call.

**Another simple ggplot, with labels**

```{r}
ggplot(CA_counties) + geom_sf() +
  geom_sf_text(aes(label = NAME), size = 1.5)

```

**and now with fill color**

```{r}
ggplot(CA_counties) + geom_sf(aes(fill = MED_AGE)) +
  geom_sf_text(aes(label = NAME), col="white", size=1.5)
```

**Repositioned legend, no "x" or "y" labels**

```{r warning=FALSE}
ggplot(CA_counties) + geom_sf(aes(fill=MED_AGE)) +
  geom_sf_text(aes(label = NAME), col="white", size=1.5) +
  theme(legend.position = c(0.8, 0.8)) +
  labs(x="",y="")

```


**Map in ggplot2, zoomed into two counties [`air_quality`]:** [@TRI]

```{r warning=FALSE}
library(tidyverse); library(sf); library(iGIScData)
census <- BayAreaTracts %>%
   filter(CNTY_FIPS %in% c("013", "095"))
TRI <- TRI_2017_CA %>%
  st_as_sf(coords = c("LONGITUDE", "LATITUDE"), crs=4326) %>%
  st_join(census) %>%
  filter(CNTY_FIPS %in% c("013", "095"),
         (`5.1_FUGITIVE_AIR` + `5.2_STACK_AIR`) > 0)
bnd = st_bbox(census)
ggplot() +
  geom_sf(data = BayAreaCounties, aes(fill = NAME)) +
  geom_sf(data = census, color="grey40", fill = NA) +
  geom_sf(data = TRI) +
  coord_sf(xlim = c(bnd[1], bnd[3]), ylim = c(bnd[2], bnd[4])) +
  labs(title="Census Tracts and TRI air-release sites") +
  theme(legend.position = "none")

```

### Rasters in ggplot2

Raster display in ggplot2 is currently a little awkward, as are rasters in general in the feature-dominated GIS world.

We can use a trick: converting rasters to a grid of points [**`marbles`**:

```{r}
library(tidyverse)
library(sf)
library(raster)
rasPath <- system.file("extdata","elev.tif", package="iGIScData")
elev <- raster(rasPath)
shpPath <- system.file("extdata","trails.shp", package="iGIScData")
trails <- st_read(shpPath)
elevpts = as.data.frame(rasterToPoints(elev))
ggplot() +
  geom_raster(data = elevpts, aes(x = x, y = y, fill = elev)) +
  geom_sf(data = trails)
```

## tmap

Basic building block is tm_shape(data) followed by various layer elements such as tm_fill()
shape can be features or raster.
See Geocomputation with R Chapter 8 "Making Maps with R" for more information.
https://geocompr.robinlovelace.net/adv-map.html

```{r}
library(spData)
library(tmap)
tm_shape(world) + tm_fill() + tm_borders()
```

**Color by variable [`air_quality`]**

```{r}
library(sf)
library(tmap)
tm_shape(BayAreaTracts) + tm_fill(col = "MED_AGE")

```

**tmap of sierraFeb with hillshade and point symbols [`sierra`]**

```{r warning=FALSE}
library(tmap)
library(sf)
library(raster)
library(iGIScData)
tmap_mode("plot")
tmap_options(max.categories = 8)
sierra <- st_as_sf(sierraFeb, coords = c("LONGITUDE", "LATITUDE"), crs = 4326)
rasPath <- system.file("extdata","ca_hillsh_WGS84.tif", package="iGIScData")
hillsh <- raster(rasPath)
bounds <- st_bbox(sierra)
tm_shape(hillsh,bbox=bounds)+
  tm_raster(palette="-Greys",legend.show=FALSE,n=10) + tm_shape(sierra) + tm_symbols(col="TEMPERATURE",
     palette=c("blue","red"), style="cont",n=8) +
  tm_legend() + 
  tm_layout(legend.position=c("RIGHT","TOP"))

```

*Note: "-Greys" needed to avoid negative image, since "Greys" go from light to dark, and to match reflectance as with b&w photography, they need to go from dark to light.*

**`UpperSinkingCoveKarst`**

From a hydrologic and geochemical study of a fluviokarstic valley system in Tennessee [@davis1993geomorphology]:

<img src="img/CaveCoveCaveEntrance.png">
<img src="img/SinkingCoveCave.png">
<img src="img/SinkingCoveCaveSteephead.png">

```{r}
library(sf); library(tidyverse); library(readxl); library(tmap)
wChemData <- read_excel(system.file("extdata/SinkingCove","SinkingCoveWaterChem.xlsx", package="iGIScData")) %>%
  mutate(siteLoc = str_sub(Site,start=1L, end=1L))
wChemTrunk <- wChemData %>% filter(siteLoc == "T") %>% mutate(siteType = "trunk")
wChemDrip <- wChemData %>% filter(siteLoc %in% c("D","S")) %>% mutate(siteType = "dripwater")
wChemTrib <- wChemData %>% filter(siteLoc %in% c("B", "F", "K", "W", "P")) %>% mutate(siteType = "tributary")
wChemData <- bind_rows(wChemTrunk, wChemDrip, wChemTrib)
sites <- read_csv(system.file("extdata/SinkingCove", "SinkingCoveSites.csv", package="iGIScData"))
wChem <- wChemData %>%
  left_join(sites, by = c("Site" = "site")) %>%
  st_as_sf(coords = c("longitude", "latitude"), crs = 4326)
library(raster)
tmap_mode("plot")
DEMpath <- system.file("extdata/SinkingCove","DEM_SinkingCoveUTM.tif",package="iGIScData")
DEM <- raster(DEMpath)
slope <- terrain(DEM, opt='slope')
aspect <- terrain(DEM, opt='aspect')
hillsh <- hillShade(slope, aspect, 40, 330)
bounds <- st_bbox(wChem)
xrange <- bounds$xmax - bounds$xmin
yrange <- bounds$ymax - bounds$ymin
xMIN <- as.numeric(bounds$xmin - xrange/10)
xMAX <- as.numeric(bounds$xmax + xrange/10)
yMIN <- as.numeric(bounds$ymin - yrange/10)
yMAX <- as.numeric(bounds$ymax + yrange/10)
#st_bbox(c(xmin = 16.1, xmax = 16.6, ymax = 48.6, ymin = 47.9), crs = st_crs(4326))
newbounds <- st_bbox(c(xmin=xMIN, xmax=xMAX, ymin=yMIN, ymax=yMAX), crs= st_crs(4326))
tm_shape(hillsh,bbox=newbounds) +
  tm_raster(palette="-Greys",legend.show=F,n=20) +
  tm_shape(DEM) + tm_raster(palette=terrain.colors(24), alpha=0.5,legend.show=F) +
  tm_shape(wChem) + tm_symbols(size="TH", col="Lithology", scale=2, shape="siteType") +
  #tm_legend(legend.outside = T) +
  tm_layout(legend.position = c("left", "bottom")) +
  tm_graticules(lines=F)
```

## Interactive Maps

The word "static" in "static maps" isn't something you would have heard in a cartography class 30 years ago, since essentially *all* maps then were static. Very important in designing maps is considering your audience, and one characteristic of the audience of those maps of yore were that they were printed and thus fixed on paper.  A lot of cartographic design relates to that property:  

- Figure-to-ground relationships assume "ground" is a white piece of paper (or possibly a standard white background in a pdf), so good cartographic color schemes tend to range from light for low values to dark for high values.
- Scale is fixed and there are no "tools" for changing scale, so a lot of attention must be paid to providing scale information.
- Similarly, without the ability to see the map at different scales, inset maps are often needed to provide context.

Interactive maps change the game in having tools for changing scale, and *always* being "printed" on a computer or device where the color of the background isn't necessarily white. We are increasingly used to using interactive maps on our phones or other devices, and often get frustrated not being able to zoom into a static map.

A widely used interactive mapping system is Leaflet, but we're going to use tmap to access Leaflet behind the scenes and allow us to create maps with one set of commands.  The key parameter needed is tmap_mode which must be set to "view" to create an interactive map. 

[**`UpperSinkingCoveKarst`**]

With an interactive map, we do have the advantage of a good choice of base maps and the ability to resize and explore the map, but symbology is more limited, mostly just color and size, with only one variable in a legend. 

```{r}
tmap_mode("view")
bounds <- st_bbox()
wChem2map <- filter(wChem, Month == 8)
minVal <- min(wChem2map$TH); maxVal <- max(wChem2map$TH)
tm_basemap(leaflet::providers$Esri.WorldTopoMap) + 
  tm_shape(wChem2map) + tm_symbols(col="siteType", size="TH", scale=2) +
  tm_layout(title=paste("Total Hardness ",as.character(minVal),"-",as.character(maxVal)," mg/L", sep=""))
tm_basemap(leaflet::providers$Esri.WorldTopoMap) + 
  tm_shape(wChem2map) + tm_symbols(col="Lithology", size="TH", scale=2)
```

[**`air_quality`**]

```{r}
tmap_mode("view")
tm_shape(BayAreaTracts) + tm_fill(col = "MED_AGE", alpha = 0.5)
```


[**`sierra`**]
```{r}
library(tmap)
library(sf)
tmap_mode("view")
tmap_options(max.categories = 8)
sierra <- st_as_sf(sierraFeb, coords = c("LONGITUDE", "LATITUDE"), crs = 4326)
bounds <- st_bbox(sierra)
tm_basemap(leaflet::providers$Esri.NatGeoWorldMap) +
  tm_shape(sierra) + tm_symbols(col="TEMPERATURE",
  palette=c("blue","red"), style="cont",n=8,size=0.2) +
  tm_legend() + 
  tm_layout(legend.position=c("RIGHT","TOP"))

```

[**`landslides`**]

slideCentroids.shp crs = 26910


#### Leaflet

Now that we've seen an app that used it, let's look briefly at Leaflet itself, and we'll see that even the Leaflet package in R actually uses JavaScript...

Leaflet is designed as "An open-source JavaScript library for mobile-friendly interactive maps"   https://leafletjs.com
"The **R** package **leaflet** is an interface to the JavaScript library **Leaflet** to create interactive web maps. It was developed on top of the htmlwidgets framework, which means the maps can be rendered in **RMarkdown** (v2) documents (which is why you can see it in this document), Shiny apps, and RStudio IDE / the R console." 

https://blog.rstudio.com/2015/06/24/leaflet-interactive-web-maps-with-r/

https://github.com/rstudio/cheatsheets/blob/master/leaflet.pdf

```{r}
library(leaflet)
m <- leaflet() %>%
  addTiles() %>%  # default OpenStreetMap tiles
  addMarkers(lng=174.768, lat=-36.852,
             popup="The birthplace of R")
m 
```

## Exercises

1. Using the method of building simple sf geometries, build a simple 1x1 square object and plot it. Remember that you have to close the polygon, so the first vertex is the same as the last (of 5) vertices. Provide your code only.

```{r include=F}
library(sf)
library(tidyverse)
s <- st_polygon(list(rbind(c(0,1),c(0,2),c(1,2),c(1,1),c(0,1))))
square <- st_sfc(s)
plot(square)
```

2. Build a map in ggplot of Colorado, Wyoming, and Utah with these boundary vertices in GCS. As with the square, remember to close each figure, and assign the crs to what is needed for GCS: 4326. Submit map as exported plot, and code in the submittal text block. [**`westUS`**]

- Colorado: (-109,41),(-102,41),(-102,37),(-109,37)
- Wyoming: (-111,45),(-104,45),(-104,41),(-111,41)
- Utah: (-114,42),(-111,42),(-111,41),(-109,41),(-109,37),(-114,37)
- Arizona: (-114,37),(-109,37),(-109,31.3),(-111,31.3),(-114.8,32.5),
           (-114.6,32.7),(-114.1,34.3),(-114.5,35),(-114.5,36),(-114,36)
- New Mexico: (-109,37),(-103,37),(-103,32),(-106.6,32),(-106.5,31.8),
           (-108.2,31.8),(-108.2,31.3),(-109,31.3)

```{r include=F}
CO <- st_polygon(list(rbind(c(-109,41),c(-102,41),c(-102,37),c(-109,37),c(-109,41))))
WY <- st_polygon(list(rbind(c(-111,45),c(-104,45),c(-104,41),c(-111,41),c(-111,45))))
UT <- st_polygon(list(rbind(c(-114,42),c(-111,42),c(-111,41),c(-109,41),c(-109,37),
                            c(-114,37),c(-114,42))))
AZ <- st_polygon(list(rbind(c(-114,37),c(-109,37),c(-109,31.3),c(-111,31.3),c(-114.8,32.5),c(-114.6,32.7),c(-114.1,34.3),c(-114.5,35),c(-114.5,36),c(-114,36),c(-114,37))))
NM <- st_polygon(list(rbind(c(-109,37),c(-103,37),c(-103,32),c(-106.6,32),c(-106.5,31.8),c(-108.2,31.8),c(-108.2,31.3),c(-109,31.3),c(-109,37))))
sfc5states <- st_sfc(CO,WY,UT,AZ,NM, crs=4326)
ggplot() + geom_sf(data=sfc5states)

```

3. Add in the code for CA and NV and create kind of a western US map...
```{r include=F}
CO <- st_polygon(list(rbind(c(-109,41),c(-102,41),c(-102,37),c(-109,37),c(-109,41))))
WY <- st_polygon(list(rbind(c(-111,45),c(-104,45),c(-104,41),c(-111,41),c(-111,45))))
UT <- st_polygon(list(rbind(c(-114,42),c(-111,42),c(-111,41),c(-109,41),c(-109,37),c(-114,37),c(-114,42))))
AZ <- st_polygon(list(rbind(c(-114,37),c(-109,37),c(-109,31.3),c(-111,31.3),c(-114.8,32.5),c(-114.6,32.7),c(-114.1,34.3),c(-114.5,35),c(-114.5,36),c(-114,36),c(-114,37))))
NM <- st_polygon(list(rbind(c(-109,37),c(-103,37),c(-103,32),c(-106.6,32),c(-106.5,31.8),c(-108.2,31.8),c(-108.2,31.3),c(-109,31.3),c(-109,37))))

CA <- st_polygon(list(rbind(c(-124,42),c(-120,42),c(-120,39),c(-114.5,35),
  c(-114.1,34.3),c(-114.6,32.7),c(-117,32.5),c(-118.5,34),c(-120.5,34.5),
  c(-122,36.5),c(-121.8,36.8),c(-122,37),c(-122.4,37.3),c(-122.5,37.8),
  c(-123,38),c(-123.7,39),c(-124,40),c(-124.4,40.5),c(-124,41),c(-124,42))))
NV <- st_polygon(list(rbind(c(-120,42),c(-114,42),c(-114,36),c(-114.5,36),
  c(-114.5,35),c(-120,39),c(-120,42))))

sfc7states <- st_sfc(CO,WY,UT,AZ,NM,CA,NV, crs=4326)
ggplot() + geom_sf(data=sfc7states)

```

4. Create an sf class from the seven states adding the fields `name`, `abb`, `area_sqkm`, and `population`, and create a map labeling with the name.

- Colorado, CO, 269837, 5758736
- Wyoming, WY, 253600, 578759
- Utah, UT, 84899, 3205958
- Arizona, AZ, 295234, 7278717
- New Mexico, NM, 314917, 2096829
- California, CA, 423970, 39368078
- Nevada, NV, 286382, 3080156

```{r include=F}
attributes <- bind_rows(c(name="Colorado", abb="CO", area=269837, pop=5758736),
                        c(name="Wyoming", abb="WY", area=253600, pop=578759),
                        c(name="Utah", abb="UT", area=84899, pop=3205958),
                        c(name="Arizona", abb="AZ", area=295234, pop=7278717),
                        c(name="New Mexico", abb="NM", area=314917, pop=2096829),
                        c(name="California", abb="CA", area=423970, pop=39368078),
                        c(name="Nevada", abb="NV", area=286382, pop=3080156))
SW_States <- st_sf(attributes, geometry = sfc7states)
ggplot(SW_States) + geom_sf(aes(fill=pop)) + geom_sf_text(aes(label = name))
```

5. Create a tibble for the highest peaks in the 7 states, with the following names, elevations in m, longitude and latitude, and add them to that map:

- Wheeler Peak, 4011, -105.4, 36.5
- Mt. Whitney, 4421, -118.2, 36.5
- Boundary Peak, 4007, -118.35, 37.9
- Kings Peak, 4120, -110.3, 40.8
- Gannett Peak, 4209, -109, 43.2
- Mt. Elbert, 4401, -106.4, 39.1
- Humphreys Peak, 3852, -111.7, 35.4

Note: the easiest way to do this is with the tribble function, starting with:
```
peaks <- tribble(
  ~peak, ~elev, ~longitude, ~latitude,
  "Wheeler Peak", 4011, -105.4, 36.5,
```

```{r include=F}
peaks <- tribble(
  ~peak, ~elev, ~longitude, ~latitude,
  "Wheeler Peak", 4011, -105.4, 36.5,
  "Mt. Whitney", 4421, -118.2, 36.5,
  "Boundary Peak", 4007, -118.35, 37.9,
  "Kings Peak", 4120, -110.3, 40.8,
  "Gannett Peak", 4209, -109, 43.2,
  "Mt. Elbert", 4401, -106.4, 39.1,
  "Humphreys Peak", 3852, -111.7, 35.4)
peaksp <- st_as_sf(peaks, coords=c("longitude", "latitude"), crs=4326)

ggplot(SW_States) + geom_sf(aes(fill=pop)) + geom_sf(data=peaksp) + geom_sf_label(data=peaksp, aes(label=peak))

  
```

6. Use a spatial join to add the points to the states to provide a new attribute maximum elevation, and display that using geom_sf_text() with the state polygons.

```{r include=F}
SW_States %>%
  st_join(peaksp) %>%
  ggplot() + geom_sf() + geom_sf_text(aes(label=elev))

```

7. From the CA_counties and CAfreeways feature data in iGIScData, make a simple map in ggplot, with freeways colored red.

```{r include=F}
ggplot(CA_counties) + geom_sf() + geom_sf(data=CAfreeways, col="red")
```

8. After adding the raster library, create a raster from the built-in `volcano` matrix of elevations from Auckland's Maunga Whau Volcano, and use plot() to display it.  We'd do more with that dataset but we don't know what the cell size is.

```{r include=F}
library(raster)
v <- raster(volcano)
plot(v)
```

9. Use tmap to create a simple map from the SW_States (polygons) and peaksp (points) data we created earlier.  Hints: you'll want to use tm_text with text set to "peak" to label the points, along with the parameter `auto.placement=TRUE`. [**`westUS`**]

```{r include=F}
library(tmap)
tmap_mode("plot")
tm_shape(SW_States) + tm_borders() +
  tm_shape(peaksp) + tm_symbols(col = "red") + tm_text(text="peak", auto.placement=T)
```

10. Change the map to the view mode, but don't use the state borders since the basemap will have them. Just before adding shapes, set the basemap to leaflet::providers$Esri.NatGeoWorldMap, then continue to the peaks after the + to see the peaks on a National Geographic basemap.

```{r include=F}
tmap_mode("view")
tm_basemap(leaflet::providers$Esri.NatGeoWorldMap) +
tm_shape(peaksp) + tm_symbols(col = "red") + tm_text(text="peak", auto.placement=T)

```