---
title: "Lecture 16" 
subtitle: "Spatial Data and Cartography (Part 2)"
date: "10/31/2018"
fontsize: 11pt
output: 
  beamer_presentation:
    theme: metropolis
    highlight: pygments
    fig_caption: false
    keep_tex: true
    latex_engine: xelatex
    includes:
      in_header: ../settings.tex
---

```{r setup, include=FALSE}
library(sf)
library(raster)

library(dplyr)
library(ggplot2)
library(patchwork)

set.seed(20180322)
```

```{r config, include=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  fig.width=7,
  fig.height=4.5,
  out.width="\\textwidth",
  fig.align="center",
  echo=TRUE,
  warning=FALSE
)

ggplot2::theme_set(ggplot2::theme_bw())

source("../util.R")
```


# Plotting

## Example Data - NC SIDS {.t}

\scriptoutput

```{r}
nc = st_read(system.file("shape/nc.shp", package="sf"), quiet = TRUE) %>% 
  select(-(AREA:CNTY_ID), -(FIPS:CRESS_ID))

tbl_df(nc)
```

## Base Plots {.t}

\scriptoutput

```{r out.width="0.8\\textwidth"}
plot(nc)
```

## Geometry Plot {.t}

\scriptoutput

```{r}
plot(st_geometry(nc), axes=TRUE)
```

## Graticules {.t}

\scriptoutput

```{r}
plot(nc[,"SID79"], axes=TRUE)
```

## Graticules {.t}

\scriptoutput

```{r}
plot(nc[,"SID79"], graticule=st_crs(nc), axes=TRUE)
```

## Graticules (EPSG:3631) {.t}

\scriptoutput

```{r}
plot(st_transform(nc[,"SID79"], 3631), axes=TRUE)
```

## Graticules (EPSG:3631) {.t}

\scriptoutput

```{r}
plot(st_transform(nc[,"SID79"], 3631), graticule=st_crs(nc), axes=TRUE)
```

## Graticules (EPSG:3631) {.t}

\scriptoutput

```{r}
plot(st_transform(nc[,"SID79"], 3631), graticule=st_crs(3631), axes=TRUE)
```


## ggplot2 {.t}

\scriptoutput

```{r}
ggplot(nc) + 
  geom_sf(aes(fill=SID79))
```

## ggplot2 + projections {.t}

\scriptoutput

```{r}
ggplot(st_transform(nc, 3631)) + 
  geom_sf(aes(fill=SID79 / BIR79))
```

## ggplot2 + viridis {.t}

\scriptoutput

```{r}
ggplot(st_transform(nc, 3631)) + 
  geom_sf(aes(fill=SID79 / BIR79)) +
  scale_fill_viridis_c()
```


## Example Data - Meuse

\scriptoutput

```{r}
data(meuse, meuse.riv, package="sp")

meuse = st_as_sf(meuse, coords=c("x", "y"), crs=28992)
meuse_riv = st_polygon(list(meuse.riv)) %>% st_sfc() %>% st_set_crs(28992)

tbl_df(meuse)
```

## Meuse {.t}

\scriptoutput

```{r}
plot(meuse, pch=16)
```



## Layering plots {.t}

\scriptoutput

```{r}
plot(meuse_riv, col=adjustcolor("lightblue", alpha.f=1), border = NA, 
     axes=TRUE, graticule=st_crs(4326), xlim=st_bbox(meuse)[c(1,3)], ylim=st_bbox(meuse)[c(2,4)])
plot(meuse[,"lead"], pch=16, add=TRUE)
```


## ggplot2 {.t}

\scriptoutput

```{r out.width="0.33\\textwidth"}
ggplot() +
  geom_sf(data=st_sf(meuse_riv), fill="lightblue", color=NA) +
  geom_sf(data=meuse, aes(color=lead), size=1)
```

## ggplot2 - axis limits {.t}

\scriptoutput

```{r out.width="0.65\\textwidth"}
ggplot() +
  geom_sf(data=st_sf(meuse_riv), fill="lightblue", color=NA) +
  geom_sf(data=meuse, aes(color=lead), size=1) +
  ylim(50.95, 50.99)
```

## ggplot2 - axis limits {.t}

\scriptoutput

```{r out.width="0.65\\textwidth"}
ggplot() +
  geom_sf(data=st_sf(meuse_riv), fill="lightblue", color=NA) +
  geom_sf(data=meuse, aes(color=lead), size=1) +
  ylim(329714, 333611)
```

## ggplot2 - bounding box{.t}

\scriptoutput

```{r out.width="0.65\\textwidth"}
ggplot() +
  geom_sf(data=st_sf(meuse_riv), fill="lightblue", color=NA) +
  geom_sf(data=meuse, aes(color=lead), size=1) +
  ylim(st_bbox(meuse)["ymin"], st_bbox(meuse)["ymax"])
```

# Geometry Manipulation

## Casting {.t}

\scriptoutput

```{r}
nc_pts = st_cast(nc, "MULTIPOINT")

tbl_df(nc_pts)
```

## {.t}

\scriptoutput

```{r}
plot(st_geometry(nc), border='grey')
plot(st_geometry(nc_pts), pch=16, cex=0.5, add=TRUE)
```

## Casting - POINT {.t}

\scriptoutput

```{r}
st_cast(nc, "POINT")
```

## {.t}

\scriptoutput

```{r}
plot(st_geometry(nc), border='grey')
plot(st_geometry(st_cast(nc, "POINT")), pch=16, cex=0.5, add=TRUE)
```

## Casting - LINESTRING {.t}

\scriptoutput

```{r}
st_cast(nc, "MULTILINESTRING") %>% as_tibble()
```

## {.t}

\scriptoutput

```{r}
st_cast(nc, "MULTILINESTRING") %>% st_geometry() %>% plot()
```


## Grouping Features {.t}

\scriptoutput

```{r}
nc_state = st_union(nc)
tbl_df(nc_state)
plot(nc_state)
```

## More Grouping {.t}

\scriptoutput

```{r}
nc_cut = nc %>%
  mutate(X = st_centroid(nc) %>% st_coordinates() %>% .[,1]) %>%
  mutate(region = cut(X, breaks = 5))

tbl_df(nc_cut)
```

## {.t}

\scriptoutput

```{r}
ggplot(nc_cut) +
  geom_sf(aes(fill=region))
```

## dplyr and sf {.t}

\scriptoutput

```{r}
nc_cut %>% 
  group_by(region) %>% 
  summarize() %>% 
  ggplot() + 
    geom_sf(aes(fill=region))
```


## Affine Transfomations {.t}

\scriptoutput

```{r}
rotate = function(a) matrix(c(cos(a), sin(a), -sin(a), cos(a)), 2, 2)

ctrd = st_centroid(nc_state)
state_rotate = lwgeom::st_make_valid( (nc_state) * rotate(-pi/4) )
plot(state_rotate, axes=TRUE)
```

## Scaling Size

\scriptoutput

```{r}
ctrd = st_centroid(st_geometry(nc))
area = st_area(nc) %>% strip_attrs()

nc_rot = nc
st_geometry(nc_rot) = (st_geometry(nc) - ctrd) * 0.75 + ctrd

plot(nc_rot[,"SID79"])
```

# Highway Example

## Highways {.t}

\scriptoutput

```{r out.height="0.5\\textheight"}
hwy = st_read("https://www2.stat.duke.edu/~cr173/Sta444_Fa18/slides/data/us_interstates.gpkg", 
              quiet=TRUE, stringsAsFactors=FALSE) %>% st_transform(st_crs(nc))
ggplot() +
  geom_sf(data=nc) +
  geom_sf(data=hwy, col='red')
```

## NC Interstate Highways {.t}

\scriptoutput

```{r}
hwy_nc = st_intersection(hwy, nc)

ggplot() + 
  geom_sf(data=nc) +
  geom_sf(data=hwy_nc, col='red')
```

## Counties near the interstate (Projection) {.t}

\scriptoutput

```{r}
nc_utm  = st_transform(nc, "+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs")

ggplot() + 
  geom_sf(data=nc_utm) +
  geom_sf(data=hwy_nc, col='red')
```

## Counties near the interstate (Buffering) {.t}

\scriptoutput

```{r}
hwy_nc_buffer = hwy_nc %>%
  st_transform("+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs") %>%
  st_buffer(10000)

ggplot() + 
  geom_sf(data=nc_utm) +
  geom_sf(data=hwy_nc, color='red') +
  geom_sf(data=hwy_nc_buffer, fill='red', alpha=0.3)
```

## Counties near the interstate (Buffering + Union) {.t}

\scriptoutput

```{r}
hwy_nc_buffer = hwy_nc %>%
  st_transform("+proj=utm +zone=17 +datum=NAD83 +units=m +no_defs") %>%
  st_buffer(10000) %>%
  st_union() %>%
  st_sf()
  
ggplot() + 
  geom_sf(data=nc_utm) +
  geom_sf(data=hwy_nc, color='red') +
  geom_sf(data=hwy_nc_buffer, fill='red', alpha=0.3)
```

## Example {.t}

How many counties in North Carolina are within 5, 10, 20, or 50 km of an interstate highway?


# Gerrymandering Example

## NC House Districts - 112th Congress

\footnoteoutput

```{r}
nc_house = st_read(
  "https://www2.stat.duke.edu/~cr173/Sta444_Fa18/slides/data/nc_districts112.gpkg", 
  stringsAsFactors = FALSE, quiet = TRUE) %>%
  select(ID, DISTRICT)

tbl_df(nc_house)
```

##

```{r}
nc_house = nc_house %>%
  st_transform("+proj=utm +zone=17 +datum=NAD83 +units=km +no_defs")
plot(nc_house[,"DISTRICT"], axes=TRUE)
```

## Measuring Compactness - Reock Score {.t}

The Reock score the a measure of compactness that is calculated as the the ratio area of a shape to the area of its minimum bounding circle.

```{r out.width="0.4\\textwidth"}
circs = nc_house %>% st_geometry() %>% lwgeom::st_minimum_bounding_circle()

sub = nc_house$DISTRICT == 1
plot(circs[sub])
plot(nc_house[sub,"DISTRICT"], add=TRUE)
```

##

```{r}
plot(nc_house[,"DISTRICT"])
plot(circs,add=TRUE)
```

## Calculating Reock

```{r out.width="0.8\\textwidth"}
nc_house = nc_house %>%
  mutate(reock = st_area(nc_house) / st_area(circs))
plot(nc_house[,"reock"])
```

##

```{r}
tbl_df(nc_house) %>% 
  arrange(reock) %>%
  print(n=13)
```


# Raster Data

## Example data - Meuse {.t}

\scriptoutput

```{r}
meuse_rast = raster(system.file("external/test.grd", package="raster"))

meuse_rast
```

## {.t}

\scriptoutput

```{r}
plot(meuse_rast)
plot(meuse_riv, add=TRUE, col=adjustcolor("lightblue",alpha.f = 0.5), border=NA)
```

## raster class

\scriptoutput
```{r}
str(meuse_rast)
```

## raster features {.t}

\scriptoutput

```{r}
extent(meuse_rast)

dim(meuse_rast)

res(meuse_rast)

projection(meuse_rast)

meuse_rast[20,]
```

## Rasters and Projections {.t}

\scriptoutput

```{r message=FALSE}
library(rgdal)
meuse_rast_ll = projectRaster(meuse_rast, crs="+proj=longlat +datum=NAD27 +no_defs")

par(mfrow=c(1,2))
plot(meuse_rast)
plot(meuse_rast_ll)
```

## {.t}

\scriptsize

```{r}
meuse_rast

meuse_rast_ll
```



## Simple Features $\longleftrightarrow$ Rasters {.t}

\scriptoutput

```{r error=TRUE}
meuse_riv_rast = rasterize(meuse_riv, meuse_rast)

meuse_riv_rast = rasterize(as(meuse_riv, "Spatial"), meuse_rast)
plot(meuse_riv_rast)
```


## Rasters and Spatial Models {.t}

\scriptoutput

```{r}
head(meuse)

head(st_coordinates(meuse))
```

## {.t}

\scriptoutput

```{r message=FALSE}
library(fields)

tps = Tps(x = st_coordinates(meuse), Y=meuse$elev)
pred_grid = xyFromCell(meuse_rast, seq_along(meuse_rast))

meuse_elev_pred = meuse_rast
meuse_elev_pred[] = predict(tps, pred_grid)

plot(meuse_elev_pred)
```

## ggplot and rasters

```{r fig.height=4}
p = rasterToPolygons(meuse_elev_pred) %>% st_as_sf()

(ggplot() + geom_sf(data=meuse, aes(color=elev), size=1)) + 
(ggplot() + geom_sf(data=p, aes(fill=test), color=NA))
```
