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

## introduction
In this lab, we'll explore the basics of working with raster data, including attribute, spatial, and geometry operations. This lab follows chapters [3](https://geocompr.robinlovelace.net/attr.html), [4](https://geocompr.robinlovelace.net/spatial-operations.html), and [5](https://geocompr.robinlovelace.net/geometry-operations.html) of Geocomputation with R by Robin Lovelace.

## prerequisites

```{r install, include=TRUE}
library(terra)
library(dplyr)
library(spData)
library(spDataLarge)
library(tmap)
library(geodata)
```

## manipulating raster objects
Raster data represents continuous surfaces, as opposed to the discrete features represented in the vector data model. Here we'll learn how to create raster data objects from scratch and how to do basic data manipulations.

Let's create a **SpatRaster** object using a digitial elevation model for Zion National Park.

```{r include=TRUE}

```

We can also create rasters from scratch using the **rast()** function. Here we create 36 cells centerd around (0,0). By default the CRS is set to WGS84, but we could change this with the "crs" argument.
Because we are working in WGS84, the resolution is in units of degrees. **rast()** fills the values of the cells row-wise starting in the upper left corner.

```{r include=TRUE}

```

The SpatRaster class can also handle multiple layers.

```{r include=TRUE}

```

We can subset layers using either the layer number or name
```{r include=TRUE}

```

We can combine SpatRaster objects into one, using the c function
```{r include=TRUE}

```

Let's create an example raster for elevation

```{r include=TRUE}

```

Rasters can also hold categorical data. Let's create an example raster for soil types.

```{r include=TRUE}

```

## raster subsetting
We can index rasters using, row-column indexing, cell IDs, coordinates, other spatial objects.

```{r include=TRUE}
# row 1, column 1

```

If we had a two layered raster, subsetting would return the values in both layers.

```{r include=TRUE}

```

We can also modify/overwrite cell values.
```{r include=TRUE}

```

Replacing values in multi-layer rasters requires a matrix with as many columns as layers and rows as replaceable cells.
```{r include=TRUE}

```

## summarizing raster objects
We can get info on raster values just by typing the name or using the summary function.
```{r include=TRUE}

```

We can get global summaries, such as standard deviation.
```{r include=TRUE}

```

Or we can use **freq()** to get the counts with categories.
```{r include=TRUE}

```

## spatial subsetting
We can move from subsetting based on specific cell IDs to extract info based on spatial objects.

To use coordinates for subsetting, we can "translate" coordinates into a cell ID with the **terra** function **cellFromXY()** or **terra::extract()**.

```{r include=TRUE}

```


Raster objects can also subset with another raster object. Here we extract the values of our elevation raster that fall within the extent of a masking raster.

```{r include=TRUE}

```

In the previous example, we just got the values back. In some cases, we might want the output to be the raster cells themselves. We can do this use the "[" operator and setting "drop = FALSE"

This example returns the first 2 cells of the first row of the "elev" raster.

```{r include=TRUE}

```

Another common use of spatial subsetting is when we use one raster with the same extent and resolution to mask the another. In this case, the masking raster needs to be composed of logicals or NAs.

```{r include=TRUE}


```

We can also use a similar approach to replace values that we suspect are incorrect.

```{r include=TRUE}

```

## map algebra
Here we define map algebra as the set of operations that modify or summarize raster cell values with reference to surrounding cells, zones, or statistical functions that apply to every cell.

## local operations
Local operations are computed on each cell individually. We can use oridinary arithemetic or logical statements.

```{r include=TRUE}

```

We can also classify intervals of values into groups. For example, we could classify a DEM into low, middle, and high elevation cells

-   first we need to construct a reclassification matrix
-   the first column corresponds to the lower end of the class
-   the second column corresponds to the upper end of the class
-   the third column corresponds to the new value for the specified ranges in columns 1 and 2

```{r include=TRUE}

```

For more efficient processing, we can use a set of map algebra functions.
-   **app()** applies a function to each cell of a raster to summarize the values of multiple layers into one layer
-   **tapp()** is an extension of "app()" that allows us to apply on operation on a subset of layers
-   **lapp()** allows us to apply a function to each cell using layers as arguments

We can use the **lapp()** function to compute the Normalized Difference Vegetation Index (NDVI).
Let's calculate NDVI for Zion National Park using multispectral satellite data.

```{r include=TRUE}

```

We need to define a function to calculate NDVI.
```{r include=TRUE}

```

So now we can use **lapp()** to calculate NDVI in each raster cell. To do so, we just need the NIR and red bands.

```{r include=TRUE}

```


## focal operations
Local operations operate on one cell, though from multiple layers. Focal operations take into account a central (focal) cell and its neighbors. The neighborhood (or kernel, moving window, filter) can take any size or shape. A focal operation applies an aggregation function to all cells in the neighborhood and updates the value of the central cell before moving on to the next central cell

We can use the **focal()** function to perform spatial filtering. We define the size, shape, and weights of the moving window using a matrix. Here we find the minimum.

```{r include=TRUE}

```

## zonal operations
Similar to focal operations, zonal operations apply an aggregation function to multiple cells. However, instead of applying operations to neighbors, zonal operations aggregate based on "zones". Zones can are defined using a categorical raster and do not necessarily have to be neighbors

For example, we could find the average elevation for different soil grain sizes.
```{r include=TRUE}

```

## merging rasters
In some cases, data for a region will be stored in multiple, contiguous files. To use them as a single raster, we need to merge them.\

In this example, we download elevation data for Austria and Switzerland and merge the two rasters into one.

```{r include=TRUE}

```



## geometric operations
When merging or performing map algebra, rasters need to match in their resolution, projection, origin, and/or extent

In the simplest case, two images differ only in their extent. Let's start by increasing the extent of a elevation raster.

```{r include=TRUE}

```

Performing algebraic operations on objects with different extents doesn't work.
```{r include=TRUE, eval=FALSE}

```

We can align the extent of the 2 rasters using the extend() function. Here we extend the elev object to the extent of elev_2 by adding NAs.

```{r include=TRUE}

```

the origin function returns the coordinates of the cell corner closes to the coordinates (0,0). We can also manually change the origin.

```{r include=TRUE}

```

## aggregation and disaggregation
Faster datasets can also differ in their resolution to match resolutions we can decrease the resolution by aggregating or increase the resolution by disaggregating.

Let's start by changing the resolution of a DEM by a factor of 5, by taking the mean.

```{r include=TRUE}

```

We have some choices when increasing the resolution. Here we try the bilinear method.
```{r include=TRUE}

```

## resampling

Aggregation/disaggregation work when both rasters have the same origins. What do we do in the case where we have two or more rasters with different origins and resolutions? Resampling computes values for new pixel locations based on custom resolutions and origins.

In most cases, the target raster would be an object you are already working with, but here we define a target raster.

```{r include=TRUE}

```




