spatial

class: inverse, center, middle
background-image: url("img/aerial_lsm.png")
background-position: 50% 50%
background-size: 75%

# Spatial data

### Maximilian H.K. Hesselbarth

#### University of Michigan (EEB)

2022/10/24

---

# Spatial data

- All data that contains **spatial information** about location

- More complex than just two extra columns (*x,y* coordinates)

- Two major data types, namely **raster** and **vector** data

.pull-left[

]

.pull-right[

<img src="img/sf.png" width="50%" style="display: block; margin: auto;" />
]

---

background-image: url("img/terra.png")
background-position: 95% 10%
background-size: 10%

# raster data

- **Cells** of equal size represent data

- Closely related to **satellite imagery**

- `terra` packages provides `SpatRaster`, `SpatVector` and `SpatExtent` classes

.yellow[(successor of formerly `raster` packages)]

---

background-image: url("img/terra.png")
background-position: 95% 10%
background-size: 10%

# `SpatRaster` data

- Multi-layer raster data

- Stores number of **columns/rows**, **spatial extent**, and **coordinate reference system**

- Values are held in RAM or on disk (if too large for memory)

```r
(x <- terra::rast(nrow = 50, ncol = 50, xmin = -50, xmax = 50, ymin = -50, ymax = 50, 
                  crs = NA))
```

```
## class       : SpatRaster 
## dimensions  : 50, 50, 1  (nrow, ncol, nlyr)
## resolution  : 2, 2  (x, y)
## extent      : -50, 50, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. :
```

```r
terra::ext(x)
## SpatExtent : -50, 50, -50, 50 (xmin, xmax, ymin, ymax)
terra::crs(x)
## [1] ""
terra::hasValues(x)
## [1] FALSE
```

---

background-image: url("img/terra.png")
background-position: 95% 10%
background-size: 10%

# `SpatRaster` data

.pull-left[

- Adding **values** and **information**

```r
vals <- runif(n = terra::nrow(x) * 
                terra::ncol(x))

*terra::values(x) <- vals
*x$lyr.2 <- vals * -1

*terra::crs(x) <- "epsg:4326"
```

```
## class       : SpatRaster 
## dimensions  : 50, 50, 2  (nrow, ncol, nlyr)
## resolution  : 2, 2  (x, y)
## extent      : -50, 50, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : memory  
##               memory  
## names       :        lyr.1,         lyr.2 
## min values  : 0.0006893259, -0.9999481156 
## max values  : 0.9999481156, -0.0006893259
```
]

.pull-right[

```r
plot(x)
```

<img src="spatial_files/figure-html/plot_rast-1.png" style="display: block; margin: auto;" />
]

---

background-image: url("img/terra.png")
background-position: 95% 10%
background-size: 10%

# `SpatRaster` data

.pull-left[

- Access data as **matrix** or **data.frame**

```r
x[]
terra::as.data.frame(x, xy = TRUE)
```

```
##          lyr.1      lyr.2
## [1,] 0.3201632 -0.3201632
## [2,] 0.5526989 -0.5526989
## [3,] 0.1201288 -0.1201288
## [4,] 0.7382276 -0.7382276
## [5,] 0.8010045 -0.8010045
## [6,] 0.1125737 -0.1125737
```

```
##     x  y     lyr.1      lyr.2
## 1 -49 49 0.3201632 -0.3201632
## 2 -47 49 0.5526989 -0.5526989
## 3 -45 49 0.1201288 -0.1201288
## 4 -43 49 0.7382276 -0.7382276
## 5 -41 49 0.8010045 -0.8010045
## 6 -39 49 0.1125737 -0.1125737
```

]

.pull-right[

- Use **mathematical** calculations with raster cells

```r
x2 <- x["lyr.1"] + x["lyr.2"]

c(x, x2)
```

```
## class       : SpatRaster 
## dimensions  : 50, 50, 3  (nrow, ncol, nlyr)
## resolution  : 2, 2  (x, y)
## extent      : -50, 50, -50, 50  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## sources     : memory  
##               memory  
##               memory  
## names       :        lyr.1,         lyr.2, lyr.1 
## min values  : 0.0006893259, -0.9999481156,     0 
## max values  : 0.9999481156, -0.0006893259,     0
```

]

---

background-image: url("img/terra.png")
background-position: 95% 10%
background-size: 10%

# National Land Cover Dataset

- Freely available at [https://www.usgs.gov](https://www.usgs.gov)

.pull-left[

```r
(nlcd <- terra::rast("data/nlcd_2019.img"))
```

```
## class       : SpatRaster 
## dimensions  : 104424, 161190, 1  (nrow, ncol, nlyr)
## resolution  : 30, 30  (x, y)
## extent      : -2493045, 2342655, 177285, 3310005  (xmin, xmax, ymin, ymax)
## coord. ref. : Albers Conical Equal Area 
## source      : nlcd_2019.img 
## color table : 1 
## categories  : NLCD Land Cover Class, Histogram, Red, Green, Blue, Opacity 
## name        :        NLCD Land Cover Class 
## min value   :                 Unclassified 
## max value   : Emergent Herbaceous Wetlands
```

]

.pull-right[

]

---

background-image: url("img/terra.png")
background-position: 95% 10%
background-size: 10%

# Cropping raster

- Crop to **area of interest**

```r
bbox <- terra::ext(-84.023438, -83.320313, 42.032974, 42.553080)

raster_crop <- terra::project(terra::rast(bbox), terra::crs(nlcd))

*nlcd_crop <- terra::crop(x = nlcd, y = raster_crop) |>
  terra::project(y = "epsg:4326")
```

---

background-image: url("img/terra.png")
background-position: 95% 10%
background-size: 10%

# Reclassification

- **Re-classify** values using classification matrix

```r
class_matrix <- cbind(c(21, 22, 23, 41, 42, 43), 
                      rep(x = c(1, 2), each = 3))

*nlcd_class <- terra::classify(nlcd_crop, rcl = class_matrix, others = NA)
```

---

background-image: url("img/sf.png")
background-position: 95% 10%
background-size: 10%

# Vector data

- Based on **points**, **lines**, and **polygons**

- Often allows realistic representation of spatial objects

- `sf` packages provides main classes: `POINT`, `LINESTRING`, `POLYGON` (and `MULTI*`)

.yellow[(successor of formerly `sp` packages)]

---

background-image: url("img/sf.png")
background-position: 95% 10%
background-size: 10%

# Read simple features

```r
cities <- sf::st_read("data/cities-mi.shp") |> 
  sf::st_transform(crs = terra::crs(nlcd))
```

```
## Reading layer `cities-mi' from data source 
##   `/Users/mhessel/Dropbox (University of Michigan)/06-Talks/Advanced-R-Workshop/slides/spatial/data/cities-mi.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1773 features and 7 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: -90.15955 ymin: 41.71626 xmax: -82.43374 ymax: 48.09634
## Geodetic CRS:  WGS84(DD)
```

```r
states <- sf::st_read("data/provinces.shp") |>
  sf::st_transform(crs = terra::crs(nlcd))
```

```
## Reading layer `provinces' from data source 
##   `/Users/mhessel/Dropbox (University of Michigan)/06-Talks/Advanced-R-Workshop/slides/spatial/data/provinces.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 51 features and 83 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -171.7911 ymin: 18.91619 xmax: -66.96466 ymax: 71.35776
## Geodetic CRS:  WGS 84
```

---

background-image: url("img/sf.png")
background-position: 95% 10%
background-size: 10%

# simple features

```
## Simple feature collection with 1773 features and 7 fields
## Geometry type: POINT
## Dimension:     XY
## Bounding box:  xmin: 448831.5 ymin: 2121120 xmax: 1096453 ymax: 2810620
## Projected CRS: Albers Conical Equal Area
## First 10 features:
##          name     type     county          gazetteer_ squaremile latitude
## 1     Merritt Township        Bay    Merritt township     31.653 43.51838
## 2      Sciota Township Shiawassee     Sciota township     26.770 42.90760
## 3    Kingston Township    Tuscola   Kingston township     35.843 43.44955
## 4       Grant Township      Mason      Grant township     48.830 44.13276
## 5  Bay De Noc Township      Delta Bay de Noc township     68.145 45.77043
## 6    Sheridan Township      Mason   Sheridan township     35.853 44.05481
## 7     Gilford Township    Tuscola    Gilford township     34.795 43.51942
## 8      Vernon  Village Shiawassee      Vernon village      1.022 42.93992
## 9      Ensign Township      Delta     Ensign township     58.441 45.87256
## 10    Cornell Township      Delta    Cornell township     60.023 45.91859
##    longitude                 geometry
## 1  -83.74868 POINT (983267.6 2344054)
## 2  -84.32242 POINT (945817.9 2270564)
## 3  -83.17502  POINT (1030069 2342515)
## 4  -86.36132 POINT (767460.5 2388076)
## 5  -86.91345 POINT (706236.7 2565161)
## 6  -86.09039 POINT (789859.7 2381650)
## 7  -83.63167 POINT (992590.2 2345386)
## 8  -84.03009 POINT (968917.9 2277091)
## 9  -86.86951 POINT (708553.8 2576777)
## 10 -87.29777 POINT (674947.3 2578747)
```

---

background-image: url("img/sf.png")
background-position: 95% 10%
background-size: 10%

# Work with simple features

- Allow to use `tidyverse` syntax

```r
(regions <- dplyr::filter(states, !name %in% c("Alaska", "Hawaii")) |>
  dplyr::select(name, region) |>
  dplyr::mutate(area_m = as.numeric(sf::st_area(geometry))) |>
  dplyr::group_by(region) |>
* dplyr::summarise(area_m = sum(area_m), area_km = area_m / 1000000) |>
  dplyr::arrange(-area_km))
```

```
## Simple feature collection with 4 features and 3 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -2357181 ymin: 311395.7 xmax: 2257301 ymax: 3167163
## Projected CRS: Albers Conical Equal Area
## # A tibble: 4 × 4
##   region     area_m  area_km                                            geometry
##   <chr>       <dbl>    <dbl>                                       <POLYGON [m]>
## 1 West      3.08e12 3075966. ((-1617641 1143814, -1746332 1223770, -1733511 124…
## 2 South     2.33e12 2326982. ((726613.2 841097.8, 654204.3 827934.5, 615792.6 8…
## 3 Midwest   2.13e12 2130036. ((-511400 1906969, -504303.7 2015504, -670478 2027…
## 4 Northeast 4.52e11  451826. ((1714981 2054464, 1708595 2047918, 1704096 203790…
```

---

background-image: url("img/sf.png")
background-position: 95% 10%
background-size: 10%

# Filter simple features

- Easy to **spatially filter** objects

```r
*row_id <- sf::st_intersects(x = states, y = dplyr::filter(cities, type == "City"),
                            sparse = FALSE)

*state_mi <- states[row_id, ]
```

---

background-image: url("img/sf.png")
background-position: 95% 10%
background-size: 10%

# Plotting `sf` objects

.pull-left[

]

.pull-right[

- Can be plotted with `ggplot2` **geom**

```r
ggplot() + 
  geom_sf(data = regions, aes(fill = area_km), alpha = 0.5) +
  geom_sf(data = state_mi, fill = NA, color = "red") +
  scale_fill_viridis() +
  theme_classic()
```

]
---

# raster-vector interactions

- **raster** and **vector** data can be used together

```r
states_mainland <- dplyr::filter(states, !name %in% c("Alaska", "Hawaii"))
plot(nlcd)
plot(states_mainland$geometry, add = TRUE)
```

---

# raster-vector interactions (I)

- `crop` cut outs parts using extent

- `mask` changes all values outside the mask

```r
nlcd_mi <- terra::crop(x = nlcd, y = state_mi) |> 
  terra::mask(mask = state_mi)
```

---

# raster-vector interactions (II)

- `extract` returns raster values at locations of vector object

```r
nlcd_cities <- terra::extract(x = nlcd_mi, 
                              y = dplyr::filter(cities, type == "Village"), 
                              cells = TRUE)
head(nlcd_cities)
```

```
##   ID        NLCD Land Cover Class      cell
## 1  1     Developed, Low Intensity 421920227
## 2  2  Developed, Medium Intensity 417355557
## 3  3     Developed, Low Intensity 370187955
## 4  4 Emergent Herbaceous Wetlands 361982561
## 5  5     Developed, Low Intensity 368848346
## 6  6     Developed, Low Intensity 443220896
```

---

class: inverse

## Thank you for your attention

### Questions?

.pull-left[

Further resources: [https://mhesselbarth.github.io/advanced-r-workshop/resources](https://mhesselbarth.github.io/advanced-r-workshop/resources)

Exercise: [https://mhesselbarth.github.io/advanced-r-workshop/exercise-spatial](https://mhesselbarth.github.io/advanced-r-workshop/exercise-spatial)

]

.pull-right[

]