---
title: "Spatial"
author: "K Arnold"
date: "11/27/2021"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, comment = "")
library(conflicted)
conflict_prefer("filter", "dplyr")
conflict_prefer("select", "dplyr")
library(tidyverse)
library(sf)
library(ggspatial)
# library(raster)
library(mapview)
library(leafem)
```

Objectives:

- Identify the basic structure of geospatial data
- Identify the basic types of geospatial visualizations
- Apply data wrangling to geospatial data

References:

- [MDSR2e chapter](https://mdsr-book.github.io/mdsr2e/ch-spatial.html)
- [sf package](https://r-spatial.github.io/sf/index.html)
- [mapview package](https://r-spatial.github.io/mapview/)
- [geocomputation with R book](https://geocompr.robinlovelace.net/location.html)

## Cholera

```{r}
# dsn <- "data/SnowGIS_SHP/"
dsn <- fs::path(fs::path_wd(), "data", "SnowGIS_SHP")
if (!dir.exists(dsn))
  usethis::use_zip("http://rtwilson.com/downloads/SnowGIS_SHP.zip", destdir = "data/", cleanup = TRUE)
st_layers(dsn)
```

```{r}
cholera_deaths <-  st_read(dsn, layer = "Cholera_Deaths")
cholera_pumps <- st_read(dsn, layer = "Pumps")
```

```{r}
st_bbox(cholera_deaths)
```

What are those?

```{r}
ggplot(cholera_deaths) +
  annotation_map_tile(type = "osm", zoomin = 0) +
  geom_sf(aes(size = Count), alpha = 0.7)
```

Oops this is incorrectly projected.

```{r}
st_crs(cholera_deaths)
```

Project this into the EPSG 4326 / WGS84 coordinate system (standard for GPS systems and Google Earth).

```{r}
cholera_4326 <- cholera_deaths %>% st_transform(4326)
st_bbox(cholera_4326)
```

But that was still wrong: we need datum. The original still was missing some data.

```{r}
st_crs(cholera_deaths)$proj4string
st_crs(27700)$proj4string
```

```{r}
cholera_latlong <- cholera_deaths %>% 
  st_set_crs(27700) %>% 
  st_transform(4326)
pumps_latlong <- cholera_pumps %>% 
  st_set_crs(27700) %>% 
  st_transform(4326)
```



```{r}
ggplot(cholera_latlong) +
  #coord_sf(lims_method = "geometry_bbox") +
  annotation_map_tile(type = "osm", zoomin = 0) +
  geom_sf(aes(size = Count), alpha = 0.7) +
  geom_sf(data = pumps_latlong, size = 3, color = "red")
```

Mapview just takes care of all of this (but it's not as clear how to control scales etc.)

```{r}
mapview(cholera_pumps, col.regions = "red", homebutton = FALSE) +
  mapview(cholera_deaths, cex = 'Count', homebutton = TRUE)
```

## NC gerrymandering

```{r}
if (!nzchar(system.file(package = "fec12")))
  remotes::install_github("baumer-lab/fec12")
library(fec12)
```

```{r}
district_elections <- results_house %>%
  mutate(district = parse_number(district_id)) %>%
  group_by(state, district) %>%
  summarize(
    N = n(), 
    total_votes = sum(general_votes, na.rm = TRUE),
    d_votes = sum(ifelse(party == "D", general_votes, 0), na.rm = TRUE),
    r_votes = sum(ifelse(party == "R", general_votes, 0), na.rm = TRUE)
  ) %>%
  mutate(
    other_votes = total_votes - d_votes - r_votes,
    r_prop = r_votes / total_votes, 
    d_prop = d_votes / total_votes,
    winner = ifelse(r_votes > d_votes, "Republican", "Democrat")
  ) %>% 
  ungroup()
nc_results <- district_elections %>% 
  filter(state == "NC")
nc_results %>%
  select(district, r_prop, winner) %>% 
  arrange(desc(r_prop))
```

```{r}
results_house
```

Data source: <https://cdmaps.polisci.ucla.edu/>


```{r}
# src <- "http://cdmaps.polisci.ucla.edu/shp/districts113.zip"
# dsn_districts <- "data/districtShapes/"
dsn_districts <- fs::path(fs::path_wd(), "data", "districtShapes")
#st_layers(dsn_districts)
districts <- st_read(dsn_districts, layer = "districts113") %>%
  mutate(DISTRICT = parse_number(as.character(DISTRICT))) %>% 
  janitor::clean_names()
```

```{r}
districts$statename %>% unique()
```


```{r}
nc_shp <- districts %>%
  filter(statename == "North Carolina")
mi_shp <- districts %>% filter(statename == "Michigan")
```

```{r}
mapview(mi_shp)
```


```{r}
mi_shp
```

```{r}
mi_shp %>% 
  inner_join(district_elections %>% filter(state == "MI")) %>% 
  mapview(zcol = "r_prop")
```


```{r}
mapview(nc_shp)
```

```{r}
nc_merged <- nc_shp %>%
  #st_transform(4326) %>%
  inner_join(nc_results, by = c("district"))
```

```{r}
districts %>% filter(statename == "Michigan") %>% 
  mapview()
```


```{r}
#pal <- leaflet::colorNumeric(palette = "RdBu", domain = c(.2, .8))
mapview(nc_merged, 
        col.regions = RColorBrewer::brewer.pal(9, "RdBu"),
        at = seq(.2, .8, length.out = 9),
        zcol = "d_prop")
```


## GR

<https://www.accesskent.com/pdfs/GIS_DataDownloadLicense.pdf>

```{r eval=FALSE}
dsn <- "~/Downloads/City_of_Grand_Rapids_Parcels/"
parcels <- st_read(dsn, layer = "City_of_Grand_Rapids_Parcels")
parcels
```

```{r eval=FALSE}
parcels %>% 
  head(50) %>% 
  ggplot() +
    annotation_map_tile() +
    geom_sf() 
```


```{r eval=FALSE}
st_geometry(parcels) %>% plot()
```


```{r eval=FALSE}
parcels %>% 
  head(10) %>% 
  plot()
```


```{r eval=FALSE}
parcels %>% nrow()
parcels %>% distinct(Address, City, Zip_Code) %>% nrow()
```


## Geocoding

[tidygeocoder](https://jessecambon.github.io/tidygeocoder/)

```{r}
addr <- tibble(address = c(
  "3201 Burton Street SE, Grand Rapids, MI, 49546", # Calvin
  "2470 Burton St SE, Grand Rapids, MI, 49546" # Fresh Thyme
  )) %>% 
  tidygeocoder::geocode(address, method = "osm")
```

```{r}
addr
```


```{r}
addr %>% st_as_sf(coords = c("long", "lat"), crs = 4326) %>% 
  ggplot() +
    annotation_map_tile() +
    geom_sf()
```

