index.Rmd

---
title: "Geographic Data Science in R"
subtitle: "Katie Jolly"
author: "Slides: katiejolly.io/rnorth-19"
date: "August 16, 2019"
institution: "Slides: katiejolly.io/rnorth-19"
output:
  xaringan::moon_reader:
    lib_dir: libs
    css: xaringan-themer.css
    nature:
      highlightStyle: github
      highlightLines: true
      countIncrementalSlides: false
---



```{r setup, include=FALSE}
options(htmltools.dir.version = FALSE)
library(tidyverse)
library(sf)

knitr::opts_chunk$set(fig.align="center", fig.width=5, fig.height=5, warning = FALSE, message = FALSE)
```

```{r xaringan-themer, include = FALSE}
library(xaringanthemer)
duo_accent(
  primary_color = "#93E1D8",
  secondary_color = "#310A31",
  header_font_google = google_font("Roboto", "400"),
  text_font_google   = google_font("Lato", "300"),
  code_font_family = "Fira Code",
  code_font_url = "https://cdn.rawgit.com/tonsky/FiraCode/1.204/distr/fira_code.css",
  header_color = "#f54278",
  title_slide_text_color = "#354a66",
  title_slide_background_image = "images/mpls.png"
)
```

### Framing Geographic Data Science (GDS)

<br> <br>

As opposed to traditional Geographic Information Science (GIS), GDS is:

<br> <br>

<ul>
<li> Interdisciplinary with Geography, Computing, and Statistics <br><br>
<li> Code-based workflow <br><br>
<li> Maximally reproducible
</ul>

.footnote[Source: https://www.robinlovelace.net/2017/05/02/can-geographic-data-save-the-world/]

---
class: center, middle


## Geographic data is a large category! 

---

background-image: url("https://pvsmt99345.i.lithium.com/t5/image/serverpage/image-id/50441iC2825B417745932D?v=1.0")
## Raster data `r emo::ji("camera")`

Think: aerial imagery, elevation models, remote sensing

.footnote[Source: https://community.alteryx.com/t5/Data-Science-Blog/Vector-and-Raster-A-Tale-of-Two-Spatial-Data-Types/ba-p/336141]

---

background-image: url("https://pvsmt99345.i.lithium.com/t5/image/serverpage/image-id/49570i26EF3FAEACD21BD4/image-size/medium?v=1.0&px=400")

## Vector data `r emo::ji("round_pushpin")`

Think: Census data, points on a map, roads

.footnote[Source: https://community.alteryx.com/t5/Data-Science-Blog/Vector-and-Raster-A-Tale-of-Two-Spatial-Data-Types/ba-p/336141]

---
class: inverse, center, middle

### Quirks (fun parts) of working with geographic data

---

### Projections and coordinate reference systems

How do you **translate** your data from a 3D shape to a 2D map...

```{r echo=FALSE}
knitr::include_graphics("https://media.opennews.org/cache/06/37/0637aa2541b31f526ad44f7cb2db7b6c.jpg")
```

---

[West Wing video](https://youtu.be/OH1bZ0F3zVU?t=34)

![](images/ww.PNG)

---

### Two common projections for US data

When working with something around the size of a **state** or **smaller**, it can often be a good idea to use state plane or UTM projections. <br>

```{r echo = FALSE}
knitr::include_graphics("images/projections.png")
```

<br>

Many states and cities also have **custom** or **recommended** projections, so it's worth doing some research before picking something! <br>

.footnote[Projection resources: [ArcGIS help](http://desktop.arcgis.com/en/arcmap/10.3/guide-books/map-projections/what-are-map-projections.htm), [Geocomputation with R: Reprojecting data](https://geocompr.robinlovelace.net/reproj-geo-data.html)]

<br>



---

class: inverse, center, middle

## Honeybee Permits in Minneapolis `r emo::ji("honeybee")`

---

## I'll talk about...

* Reading in spatial data

* Reprojecting data

* Basic maps

* Spatial join

* Neighborhood definition

* Spatial clustering (Moran's I)

* Modifiable areal unit problem


---

## Reading in spatial data

### Shapefiles

The most common file format is the **shapefile**, which is actually a collection of files. It's important to keep all of these parts in one directory! 

```{r echo = FALSE}
knitr::include_graphics("images/honeybee-files.PNG")
```

But when you read in the data, it only looks like you're using the *.shp* file.

```{r results = 'hide', eval = FALSE}
library(sf)
honeybees <- st_read("data/honeybees/Honey_Bee_Permits_2017.shp")
```


---

## Reading in spatial data

### APIs

If you're getting data from somewhere like an open data portal, using the API endpoint can often be easier.

* Easier for people trying to run your file

* Easier for your file management

```{r results = 'hide'}
honeybees <- st_read("https://opendata.arcgis.com/datasets/f99ce43936d74f718e92a37a560ad875_0.geojson")
```

There are reasons for one way over another, but I prefer APIs when possible. Either way you'll get the same data.

---

## Reprojecting data

We should first check the current projection.


```{r}
st_crs(honeybees)
```





--
When I look at this, I notice +proj=<mark>longlat</mark>, which is a geographic, not projected, coordinate system.

--


I'll use a UTM projection.

```{r}
honeybees_utm <- honeybees %>%
  st_transform(26915) # UTM 15N zone
st_crs(honeybees_utm)
```

---

## Basic maps

`sf` objects have a `plot` function. 

```{r}
plot(honeybees_utm %>% dplyr::select(HiveType))
```

---

## Basic maps

But they also work nicely with `ggplot2`

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

---

## Basic maps

We can use the ggplot layering logic to add contextual data, like Minneapolis neighborhood boundaries.

```{r results='hide'}
neighborhoods <- st_read("https://opendata.arcgis.com/datasets/055ca54e5fcc47329f081c9ef51d038e_0.geojson") %>%
  st_transform(26915)
```

```{r echo = FALSE}
library(showtext)
font_add_google("Lato", regular.wt = 300, bold.wt = 400)

showtext_auto()
ggplot() +
  geom_sf(data = neighborhoods, fill = "#D8E2DC", color = "white") +
  geom_sf(data = honeybees, color = "#FCFC62", size = 2) +
  theme_minimal() +
  theme(axis.text = element_blank(),
        panel.grid.major = element_line("transparent"),
        text = element_text(family = "Lato")) +
  labs(title = "Honeybee Permits in Minneapolis", caption = "Source: Minneapolis Open Data \n Projection: UTM 15N")
```

---

## Basic maps

```{r echo = FALSE}
library(showtext)
font_add_google("Lato", regular.wt = 300, bold.wt = 400)

showtext_auto()
ggplot() +
  geom_sf(data = neighborhoods, fill = "#D8E2DC", color = "white") +
  geom_sf(data = honeybees, color = "#FCFC62", size = 2) +
  theme_minimal() +
  theme(axis.text = element_blank(),
        panel.grid.major = element_line("transparent"),
        text = element_text(family = "Lato")) +
  labs(title = "Honeybee Permits in Minneapolis", caption = "Source: Minneapolis Open Data \n Projection: UTM 15N")
```

<br> 

When I look at this map as a geographer, I look for patterns of clustering or dispersion. I'll now walk through how to quantify that. 

---

## Spatial join

One way to think about clustering is to ask whether or not the permits are clustered **by neighborhood**. 

```{r echo=FALSE, fig.width=5, fig.height=3}
knitr::include_graphics("https://i.stack.imgur.com/CVVSH.png")
```


---

## Spatial join

```{r}
permits_per_nb <- neighborhoods %>% 
  st_join(honeybees_utm) %>% # which neighborhood is each permit in?
  group_by(BDNAME) %>% # and when we sum by neighborhood
  summarise(permits = sum(!is.na(OBJECTID))) # how many permits total?

summary(permits_per_nb$permits)
```

<br>

On average, there are **1.034** honeybee permits per neighborhood. 

<br>

But is this equally likely everywhere in the city? Or are permits more likely to be in certain areas?

---

## Spatial join

```{r echo = FALSE, fig.width=6, fig.height=6}
ggplot(permits_per_nb) +
  geom_sf(aes(fill = permits), color = "gray50", lwd = 0.3) +
  theme_minimal() +
  theme(axis.text = element_blank(),
        panel.grid.major = element_line("transparent"),
        text = element_text(family = "Lato", size = 12)) +
  labs(title = "Honeybee Permits in Minneapolis Neighborhoods", caption = "Source: Minneapolis Open Data \n Projection: UTM 15N") +
  scale_fill_gradient(low = "#fff7f3", high = "#7a0177") +
  guides(fill = guide_colorbar(title = "Permits", ticks = FALSE, position = "bottom", nbin = 5))
```


---

## Neighbor definition


```{r echo=FALSE}
knitr::include_graphics("http://geohealthinnovations.org/wp-content/uploads/2013/01/toblerquote.png")
```

---

## Neighbor definition

One common and straightforward way to define neighbors is **queen contiguity**.

```{r echo=FALSE}
knitr::include_graphics("https://encrypted-tbn0.gstatic.com/images?q=tbn:ANd9GcRxG5jeEH-7MmKzlzOcooJMOeMAMzZKqrgePBFpNP43w9W8ACq35g")
```

```{r}
library(spdep)

# need to convert to SpatialPolygons object
honeybees_sp <- as(permits_per_nb, "Spatial") 

# define neighbor structure
nb_queen <- poly2nb(honeybees_sp, queen = TRUE)

# create a matrix of binary spatial weights 
# (connected or not connected)
weights <- nb2listw(nb_queen, style = "B") 
```

---

## Neighbor definiton

```{r echo = FALSE}
coords <- coordinates(honeybees_sp)

plot(honeybees_sp, border = "gray50", col = "white")
plot(nb_queen, coords, add = TRUE, col = "#7a0177")
```


On average, each of the neighborhoods in Minneapolis has **5.6** queen's case neighbors.

---

## Spatial clustering

<br>
<br>

```{r echo = FALSE, out.width="500px"}
knitr::include_graphics("https://encrypted-tbn0.gstatic.com/images?q=tbn:ANd9GcTlm-7-YfREIPTXXfloomGr0jdZk1GhPimm7WH8lZGWwVQIHzDA")
```

---

## Spatial clustering

### Moran's I

> The Moran’s I statistic is the correlation coefficient for the relationship between a variable [like honeybee permits] and its surrounding values (Gimond)


```{r echo = FALSE, out.width="400px"}
knitr::include_graphics("https://mgimond.github.io/Spatial/img/MoranI_scatter_plot.png")
```




---

## Spatial clustering

### Moran's I

```{r}
set.seed(123)

# 10000 simulations
moran.mc(honeybees_sp$permits, weights, nsim=9999) 
```

---

## Spatial clustering

### Moran's I

statistic = <mark>0.13366</mark>, observed rank = 9869, p-value = <mark>0.0131</mark>

Based on this, we reject our null hypothesis and say that the permits are **slightly clustered** across neighborhoods.

```{r, echo = FALSE, out.width="250px"}
m <- moran.mc(honeybees_sp$permits, weights, nsim=9999) 

bee.lag = lag.listw(weights, honeybees_sp$permits)

ggplot(mapping = aes(x = honeybees_sp$permits, y = bee.lag)) +
  geom_point() +
  theme_minimal() +
  labs(y = "Neighbor honeybee permits", x ="Honeybee permits") +
  geom_smooth(method = "lm", fill = "#dee1e3", color = "#f7e02d") 
```

Does this make sense?

---

## Spatial clustering

```{r echo = FALSE, fig.width=6, fig.height=6}
ggplot(permits_per_nb) +
  geom_sf(aes(fill = permits), color = "gray50", lwd = 0.3) +
  theme_minimal() +
  theme(axis.text = element_blank(),
        panel.grid.major = element_line("transparent"),
        text = element_text(family = "Lato", size = 12)) +
  labs(title = "Honeybee Permits in Minneapolis Neighborhoods", caption = "Source: Minneapolis Open Data \n Projection: UTM 15N") +
  scale_fill_gradient(low = "#fff7f3", high = "#7a0177")  +
  guides(fill = guide_colorbar(title = "Permits", ticks = FALSE, position = "bottom"))
```


---

## Modifiable areal unit problem

Grouping points by an areal unit may distort or exaggerate the actual data pattern!

```{r echo = FALSE, out.width="200px"}
knitr::include_graphics("http://4.bp.blogspot.com/-_pvMVmjCTQU/TsY9QaWcJaI/AAAAAAAAA30/TGhvGkcnMPY/s1600/Penn+State+Map.gif")
```

---

## Modifiable areal unit problem

What if we look at honeybee permits by community instead?

```{r results = 'hide', echo = FALSE}
bg_mpls <- st_read("https://opendata.arcgis.com/datasets/e0a3b6a4c23b4a03b988388553cb9638_0.geojson") %>% 
  st_transform(26915) 

bg_join <- st_join(bg_mpls, honeybees_utm)

bg <- bg_join %>%
  group_by(CommName) %>%
  summarise(permits = sum(!is.na(OBJECTID))) 
```


```{r echo = FALSE, fig.width=6, results = 'show'}
ggplot(bg) +
  geom_sf(aes(fill = permits), color = "gray50", lwd = 0.3) +
  theme_minimal() +
  theme(axis.text = element_blank(),
        panel.grid.major = element_line("transparent"),
        text = element_text(family = "Lato", size = 12)) +
  labs(title = "Honeybee Permits in Minneapolis Block Groups", caption = "Source: Minneapolis Open Data & US Census Bureau \n Projection: UTM 15N") +
  scale_fill_gradient(low = "#fff7f3", high = "#7a0177")  +
  guides(fill = guide_colorbar(title = "Permits", ticks = FALSE, position = "bottom"))

```

Does this look more or less clustered?

---

## Modifiable areal unit problem

According to Moran's I, it's spatially random. But we found that it was clustered by neighborhood?

<br>
<br>

```{r echo = FALSE}
bg_sp <- as(bg, "Spatial") # convert our data to sp data

bg_obj <- poly2nb(bg_sp) # create the neighborhood object

weights <- nb2listw(bg_obj, style = "B") # create a matrix of binary spatial weights (connected or not connected)

set.seed(123)
moran.mc(bg_sp$permits, weights, nsim=9999) # 10000 simulations
```

<br>
<br>

What could be some consequences of this?


---

## Modifiable areal unit problem

Gerrymandering, one of my research areas, is largely an application of the modifiable areal unit problem. How we carve up our space matters!

```{r echo = FALSE, out.width="500x"}
knitr::include_graphics("https://pbs.twimg.com/media/B-8ljgjU0AASq8g.jpg")
```

---

class: center, middle

### Thank you for listening and happy to answer any questions!

@katiejolly6

katiejolly6@gmail.com

katiejolly.io