---
title: "Solution _Spatial data_"
author: "Maximilian H.K. Hesselbarth"
date: 2022/10/24
editor_options: 
  chunk_output_type: console
---

```{r, setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE
)

library(downloadthis)
```

```{r tidyfigure, echo = FALSE, fig.align = "center", out.width = '65%'}
knitr::include_graphics("img/sf.png", auto_pdf = FALSE)
```

```{r download, echo = FALSE}
download_link(link = "https://raw.githubusercontent.com/mhesselbarth/advanced-r-workshop/main/solution-spatial.Rmd",
              button_label = "Download .Rmd file", 
              button_type = "danger")
```

<br><br>

Make sure you can install and load all packages. This includes `terra` and `sf`, but also the `tidyverse`.

```{r load_libs}
library(tidyverse)

library(sf)
library(terra)
```

Next, go to [https://www.naturalearthdata.com](https://www.naturalearthdata.com) and download the "Small scale data, 1:110m" > "Cultural" > "Admin 1 – States, Provinces" data set. Additionally, download the 
"NLCD 2019 Land Cover (CONUS)" data set from [https://www.mrlc.gov](https://www.mrlc.gov). 

Once you downloaded all the data, read it into your R Session using the corresponding packages.

```{r read}
states <- sf::read_sf("slides/spatial/data/provinces.shp")

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

Make sure both the vector and the raster data have the same CRS (Hint: It's often faster to project vectors instead of raster. If projecting the raster, have a look at the 'method' argument).

```{r crs}
states <- sf::st_transform(x = states, crs = terra::crs(nlcd))
```

Next, remove Alaska and Hawaii from the states vector because there is no NLCD data for these states. Next select only the 5 largest states in area

```{r states_size}
states_large <- dplyr::filter(states, !name %in% c("Alaska", "Hawaii")) |> 
  dplyr::mutate(area_m = as.numeric(sf::st_area(geometry))) |> 
  dplyr::arrange(-area_m) |> 
  head(n = 5)
```

First plot the NLCD data and add the largest states to the map. Try to use the region as shape fill.

```{r plot}
plot(nlcd)

plot(states_large["region"], add = TRUE)
```

Now, pick one state (your home state, a state you recently visited, a state you want to visit, ...) and get the NLCD data for that state only.

```{r crop}
kansas <- dplyr::filter(states, name == "Kansas")

nlcd_kansas <- terra::crop(x = nlcd, y = kansas) |> 
  terra::mask(mask = kansas)
```

Next, get all values of the cropped NLCD data and remove all `NA` and `NaN` values. Calculate the relative amount of all remaining values. Which one is the most dominant land-cover class in your state?

```{r rel_class}
values_kansas <- terra::values(nlcd_kansas)

values_kansas <- values_kansas[is.finite(values_kansas)]

rel_classes <- sort((table(values_kansas) / length(values_kansas)) * 100, 
                    decreasing = TRUE)
```

Last, try to reclassify the raster into less classes (e.g., use the bigger classification found at[NLCD classes](https://www.mrlc.gov/data/legends/national-land-cover-database-class-legend-and-description))

```{r classification}
classes_present <- names(rel_classes)

class_mat <- cbind(as.numeric(classes_present), 
                   as.numeric(stringr::str_sub(classes_present, start = 1, end = 1)))

kansas_reclass <- terra::classify(x = nlcd_kansas, rcl = class_mat)
```