---
title: "Spatial Data and Cartography (Part 2)"
subtitle: "Lecture 18"
author: "Dr. Colin Rundel"
footer: "Sta 344 - Fall 2022"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
execute:
  echo: true
  warning: false
  collapse: true
---

```{r setup, include=FALSE}
library(raster)
library(tidyverse)
library(patchwork)
library(sta344)

library(sf)
library(stringr)

knitr::opts_chunk$set(
  fig.align = "center"
)

options(width=60)

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

set.seed(20221102)
```
# Plotting

## Example Data - NC SIDS

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

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

## Geometry Plot

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

## Graticules

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

## Graticules

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

## Graticules (EPSG:3631)

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

## Graticules (EPSG:3631)

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

## Graticules (EPSG:3631)

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

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

## ggplot2 + projections

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

## ggplot2 + viridis

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

```{r}
data(meuse, meuse.riv, package="sp")
(meuse = st_as_sf(meuse, coords=c("x", "y"), crs=28992) |>
  as_tibble() |> st_as_sf())
```

##

```{r}
( meuse_riv = st_polygon(list(meuse.riv)) |>
    st_sfc() |>
    st_set_crs(28992) |>
    st_as_sf()
)
```

## Meuse

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

## Layering plots

::: {.small}
```{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

```{r}
ggplot() +
  geom_sf(data=meuse_riv, fill="lightblue", color=NA) +
  geom_sf(data=meuse, aes(color=lead), size=1)
```


## ggplot2 - axis limits

```{r}
ggplot() +
  geom_sf(data=meuse_riv, fill="lightblue", color=NA) +
  geom_sf(data=meuse, aes(color=lead), size=1) +
  ylim(50.95, 50.99)
```


## ggplot2 - axis limits

```{r}
ggplot() +
  geom_sf(data=meuse_riv, fill="lightblue", color=NA) +
  geom_sf(data=meuse, aes(color=lead), size=1) +
  ylim(329714, 333611)
```


## ggplot2 - bounding box{.t}

```{r}
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 Predicates

```{r include=FALSE}
ncc  = read_sf("data/gis/nc_counties/", quiet=TRUE) |> select(COUNTY)
air = read_sf("data/gis/airports/", quiet=TRUE) |> st_transform(st_crs(ncc))
hwy = read_sf("data/gis/us_interstates/", quiet=TRUE) |> st_transform(st_crs(ncc))
```

## DE-9IM

![](imgs/de_9im.png){fig-align="center" width="100%"}



## Spatial predicates 

![](imgs/predicates.png){fig-align="center" width="100%"}

::: {.small}
`st_within(a,b)`:

$$
\begin{bmatrix} 
T & * & F \\ 
* & * & F \\ 
* & * & * 
\end{bmatrix}
$$ 

 `st_touches(a,b)`:


$$
\begin{bmatrix} F & T & * \\ * & * & * \\ * & * & * \end{bmatrix} \cup 
\begin{bmatrix} F & * & * \\ T & * & * \\ * & * & * \end{bmatrix} \cup 
\begin{bmatrix} F & * & * \\ * & T & * \\ * & * & * \end{bmatrix}
$$
:::

## Sparse vs Full Results

::: {.small}
```{r}
st_intersects(ncc[20:30,], air) %>% str()
```
:::

. . .

::: {.small}
```{r}
st_intersects(ncc, air, sparse=FALSE) %>% str()
```
:::

## Examples

* Which counties have an airport?

* Which counties are adjacent to Durham County?

* Which counties have more than 4 neighbors?

::: {.small}
```{r}
ncc  = read_sf("data/gis/nc_counties/", quiet=TRUE)
air = read_sf("data/gis/airports/", quiet=TRUE) |> st_transform(st_crs(ncc))
hwy = read_sf("data/gis/us_interstates/", quiet=TRUE) |> st_transform(st_crs(ncc))
```
:::


## Which counties have an airport?

::: {.small}
```{r}
ncc |>
  mutate(
    airports = st_intersects(ncc, air) |> strip_attrs(),
    n_airports = purrr::map_int(airports, length)
  ) |>
  filter(n_airports > 0)
```
:::

## Which counties are adjacent to Durham County?

::: {.small}
```{r}
ncc |>
  mutate(
    touch_durham = st_touches(ncc, ncc |> filter(COUNTY == "Durham County")) |> strip_attrs(),
    n_touches = map_int(touch_durham, length)
  ) |>
  filter(n_touches > 0)
```
:::

## Which counties have more than 4 neighbors?

::: {.small}
```{r}
ncc |>
  mutate(
    neighbors = st_touches(ncc) |> strip_attrs(),
    n_neighbors = map_int(neighbors, length)
  ) |>
  ggplot(aes(fill = n_neighbors > 4)) +
    geom_sf()
```
:::



# Geometry Manipulation


## Casting

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

##

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

## Casting - POINT

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

##

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

## Casting - LINESTRING

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

##

```{r}
st_cast(nc, "MULTILINESTRING") |> st_geometry() |> plot()
```
## Grouping Features

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

## More Grouping

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

##

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

## dplyr and sf

```{r}
nc_cut |> 
  group_by(region) |> 
  summarize() |> 
  ggplot() + 
    geom_sf(aes(fill=region))
```
## Affine Transfomations

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

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

## Scaling Size

```{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

```{r out.height="0.5\\textheight"}
hwy = read_sf("data/us_interstates.gpkg", quiet = TRUE) |>
  st_transform(st_crs(nc))

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

## NC Interstate Highways

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

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

## Counties near the interstate (Projection)

```{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)

```{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)

```{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

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

::: {.medium}
```{r}
( nc_house = read_sf("data/nc_districts112.gpkg", quiet = TRUE) |>
    select(ID, DISTRICT) )
```
:::

##

```{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

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

::: {.small}
```{r}
circs = nc_house |> 
  lwgeom::st_minimum_bounding_circle()

plot(circs |> filter(DISTRICT == 1) |> st_geometry(), axes=TRUE)
plot(nc_house |> select(DISTRICT) |> filter(DISTRICT == 1), add=TRUE)
```
:::

##

::: {.small}
```{r}
ggplot(mapping = aes(fill=DISTRICT)) +
  geom_sf(data=nc_house) +
  geom_sf(data=circs, alpha=0.15) +
  guides(color="none", fill="none")
```
:::

## Calculating Reock

```{r}
nc_house |>
  mutate(reock = (st_area(nc_house) / st_area(circs)) |> as.numeric()) |>
  ggplot(aes(fill = reock)) +
    geom_sf()

```

##

::: {.small}
```{r}
nc_house |>
  mutate(reock = st_area(nc_house) / st_area(circs)) |>
  arrange(reock) |>
  print(n=13)
```
:::

# Raster Data (stars)

## Example data - Meuse

```{r}
( meuse_rast = stars::read_stars(
    system.file("external/test.grd", package="raster")
  ) |>
    st_transform(st_crs(meuse_riv))
 )
```

## stars class

::: {.small}
```{r}
str(meuse_rast)
```
:::


## Plotting

```{r}
plot(meuse_rast)
```


##

::: {.small}
```{r}
plot(
  meuse_riv, 
  col=adjustcolor("lightblue",alpha.f = 0.5), border=NA,
  ylim = c(329500, 333611), axes=TRUE
)
plot(meuse_rast, add=TRUE)
```
:::


## ggplot

```{r}
ggplot() +
  stars::geom_stars(data=meuse_rast)
```

##

```{r}
ggplot() +
  stars::geom_stars(data=meuse_rast) +
  geom_sf(data=meuse_riv, fill="lightblue", color=NA, alpha=0.5) +
  ylim(329500, 333611)
```

##

```{r}
ggplot() +
  stars::geom_stars(data=meuse_rast) +
  geom_sf(data=meuse_riv, fill="lightblue", color=NA, alpha=0.5) +
  ylim(329500, 333611) +
  scale_fill_gradient(na.value = NA)
```

## Rasters and Projections

::: {.small}
```{r}
meuse_rast_ll = st_transform(meuse_rast, "+proj=longlat +datum=NAD83 +no_defs")
```
:::

:::: {.columns .small}
::: {.column width='50%'}
```{r}
plot(meuse_rast, axes=TRUE)
```
:::

::: {.column width='50%'}
```{r message=FALSE}
plot(meuse_rast_ll, axes=TRUE)
```
:::
::::




##

:::: {.columns .small}
::: {.column width='50%'}
```{r}
meuse_rast
```

:::

::: {.column width='50%'}
```{r}
meuse_rast_ll
```
:::
::::

## Cropping

```{r}
meuse_rast_riv = meuse_rast[ meuse_riv ]
plot(meuse_rast_riv, axes=TRUE)
```

## Segmentation

```{r}
meuse_rast |>
  mutate(bins = cut(test.grd, 3) ) |>
  select(bins) |>
  plot()
```

## Polygonization

::: {.small}
```{r}
meuse_rast_poly = meuse_rast |>
  mutate(bins = cut(test.grd, 3) ) |>
  select(bins) |>
  st_as_sf()
plot(meuse_rast_poly)
```
:::


##

::: {.small}
```{r}
meuse_rast_poly |>
  group_by(bins) |>
  summarize() |>
  plot()
```
:::


##

::: {.small}
```{r}
meuse_rast_poly |>
  group_by(bins) |>
  summarize() |>
  mutate(area = st_area(geometry))
```
:::
