Working with Geospatial Data

class: center, middle, inverse, title-slide

# Working with Geospatial Data
## Joining spatial data
### Ben Baumer
### SDS 192</br>April 15th, 2020</br>(<a href="http://beanumber.github.io/sds192/lectures/mdsr_geo_09-joins.html" class="uri">http://beanumber.github.io/sds192/lectures/mdsr_geo_09-joins.html</a>)

---

class: center, middle, inverse
# [Spatial joins](http://wiki.gis.com/wiki/index.php/Spatial_Join)

---
## Two ideas

.pull-left[

- Joining based on **attributes**

- `sf` object + `data.frame` object

- add new variables

- `inner_join()`, `left_join()`, etc.

]

.pull-left[

- Joining based on **geometries**

- `sf` object + `sf` object

- make new geometries

- `st_join()`

- `left = TRUE/FALSE`

- `join` predicates:
    - `st_within`, `st_crosses`, `st_intersects`, etc.

.footnote[http://wiki.gis.com/wiki/index.php/Spatial_Join]

]

---
## Make a choropleth of the Northeast

```r
library(tigris)
library(sf)

us <- states(cb = TRUE, class = "sf")
```

```r
state_region_lkup <- tibble(
  postal_code = state.abb, 
  region = state.region
)
```

```r
class(us)
```

```
## [1] "sf"         "data.frame"
```

```r
class(state_region_lkup)
```

```
## [1] "tbl_df"     "tbl"        "data.frame"
```

---
## Attribute (non-spatial) join

```r
nrow(us)
```

```
## [1] 56
```

```r
nrow(state_region_lkup)
```

```
## [1] 50
```

```r
ne <- us %>%
  left_join(state_region_lkup, by = c("STUSPS" = "postal_code")) %>%
  filter(region == "Northeast")

nrow(ne)
```

```
## [1] 9
```

---
## Land area in the Northeast

.pull-left[

```r
ne_plot <- ggplot(ne) +
  geom_sf(
    aes(fill = ALAND / 1e6)
  ) + 
  coord_sf(crs = st_crs(2163)) + 
  scale_fill_distiller(
    "Land Area (square km)", 
    palette = "Greens", 
    direction = 1
  )
```

.footnote[https://epsg.io/2163]

]

.pull-right[

![](figures/northeast-land-1.png)

]

---
## Spatial join

```r
library(macleish)
macleish <- macleish_layers %>%
  pluck("boundary")
```

```r
ne %>%
* st_join(macleish)
```

```
## Error in st_geos_binop("intersects", x, y, sparse = sparse, prepared = prepared): st_crs(x) == st_crs(y) is not TRUE
```

```r
ne %>%
  st_crs()
```

```
## Coordinate Reference System:
##   User input: 4269 
##   wkt:
## GEOGCS["GCS_North_American_1983",
##     DATUM["North_American_Datum_1983",
##         SPHEROID["GRS_1980",6378137,298.257222101]],
##     PRIMEM["Greenwich",0],
##     UNIT["Degree",0.017453292519943295],
##     AUTHORITY["EPSG","4269"]]
```

---
## Convert to common CRS first!

```r
ne %>%
  st_transform(4326) %>%
* st_join(macleish, join = st_contains, left = FALSE)
```

```
## although coordinates are longitude/latitude, st_contains assumes that they are planar
```

```
## Simple feature collection with 1 feature and 13 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: -73.50814 ymin: 41.23796 xmax: -69.92839 ymax: 42.88659
## CRS:            EPSG:4326
##   STATEFP  STATENS    AFFGEOID GEOID STUSPS          NAME LSAD       ALAND
## 1      25 00606926 0400000US25    25     MA Massachusetts   00 20205125364
##       AWATER    region OBJECTID Shape_Leng Shape_Area
## 1 7129925486 Northeast        1   5893.578    1033988
##                         geometry
## 1 MULTIPOLYGON (((-70.23405 4...
```

---
## Can change join predicate!

```r
ne %>%
  st_transform(4326) %>%
* st_join(macleish, join = st_within, left = FALSE)
```

```
## although coordinates are longitude/latitude, st_within assumes that they are planar
```

```
## Simple feature collection with 0 features and 13 fields
## bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
## CRS:            EPSG:4326
##  [1] STATEFP    STATENS    AFFGEOID   GEOID      STUSPS     NAME      
##  [7] LSAD       ALAND      AWATER     region     OBJECTID   Shape_Leng
## [13] Shape_Area geometry  
## <0 rows> (or 0-length row.names)
```

.footnote[https://r-spatial.github.io/sf/reference/st_join.html]