Session4c.Rmd

---
source: Rmd
title: "Intersecting raster and vector data"
author: "Clémentine Cottineau"
teaching: 40
exercises: 20
---


```{r, include=FALSE}
source("bin/chunk-options.R")
knitr_fig_path("04c-")
```

# Question: "How can I crop raster objects to vector objects, and extract the summary of raster pixels?"

## Objectives:
- Crop a raster to the extent of a vector layer.
- Extract values from a raster that correspond to a vector file overlay.

## Keypoints:
- Use the `crop()` function to crop a raster object.
- Use the `extract()` function to extract pixels from a raster object that fall within a particular extent boundary.
- Use the `extent()` function to define an extent.

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

```{r lib}
library(tidyverse)
library(sf)
library(rgdal)
library(raster)
library(here)
```

---

# Introduction

**copied from the carpentry lesson [Manipulating Raster Data](https://datacarpentry.org/r-raster-vector-geospatial/11-vector-raster-integration/index.html)).**

We often work with spatial layers that have different spatial extents. The spatial extent of a shapefile or R spatial object represents the geographic “edge” or location that is the furthest north, south east and west. Thus is represents the overall geographic coverage of the spatial object.

![Image Source: National Ecological Observatory Network (NEON)](fig/spatial_extent.png)
The graphic below illustrates the extent of several of the spatial layers that we have worked with in this workshop:

- Area of interest (AOI) – blue
- Roads and trails – purple
- Vegetation plot locations (marked with white dots)– black
- A canopy height model (CHM) in GeoTIFF format – green

![Image Source: DCC](fig/rmd-11-compare-data-extents-1.png)

Frequent use cases of cropping a raster file include reducing file size and creating maps. Sometimes we have a raster file that is much larger than our study area or area of interest. It is often more efficient to crop the raster to the extent of our study area to reduce file sizes as we process our data. Cropping a raster can also be useful when creating pretty maps so that the raster layer matches the extent of the desired vector layers.

## Import the raster 


We import a *DSM* (Digital Surface Model) and a *DTM* (Digital Terrain Model) of TU-Delft.

The difference between the two altitudes gives us a *CHM*, a Canopy Height Model.


```{r}
DSM_TUD <- raster(here("data","tud-dsm.tif"))
DTM_TUD <- raster(here("data","tud-dtm.tif"))
CHM_TUD <- DSM_TUD - DTM_TUD

CHM_TUD_df <- as.data.frame(CHM_TUD, xy = TRUE)

oai_boundary_tudlib <- st_as_sfc(st_bbox(raster(here("data","tudlib-rgb.tif"))))
```



# Crop a Raster Using Vector Extent

We can use the `crop()` function to crop a raster to the extent of another spatial object. To do this, we need to specify the raster to be cropped and the spatial object that will be used to crop the raster. R will use the extent of the spatial object as the cropping boundary.

To illustrate this, we will crop the Canopy Height Model (CHM) to only include the area of interest (AOI). Let’s start by plotting the full extent of the CHM data and overlay where the AOI falls within it. The boundaries of the AOI will be colored blue, and we use `fill = NA` to make the area transparent.



```{r from-ep2}

ggplot() +
  geom_raster(data = CHM_TUD_df, aes(x = x, y = y, fill = layer)) + 
  scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) +
  geom_sf(data = oai_boundary_tudlib, color = "blue", fill = NA) +
  coord_sf()

```

Now that we have visualized the area of the CHM we want to subset, we can perform the cropping operation. We are going to use the `crop()` function from the `raster` package to create a new object with only the portion of the CHM data that falls within the boundaries of the AOI.


```{r}
CHM_TUD_Cropped <- crop(x = CHM_TUD, y = st_as_sf(oai_boundary_tudlib))
```

Now we can plot the cropped CHM data, along with a boundary box showing the full CHM extent. However, remember, since this is raster data, we need to convert to a data frame in order to plot using ggplot. To get the boundary box from CHM, the `st_bbox()` will extract the 4 corners of the rectangle that encompass all the features contained in this object. The `st_as_sfc()` converts these 4 coordinates into a polygon that we can plot:



```{r}
CHM_TUD_Cropped_df <- as.data.frame(CHM_TUD_Cropped, xy = TRUE)

ggplot() +
  geom_sf(data = st_as_sfc(st_bbox(CHM_TUD)), fill = "green",
          color = "green", alpha = .2) +  
  geom_raster(data = CHM_TUD_Cropped_df,
              aes(x = x, y = y, fill = layer)) + 
  scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) + 
  coord_sf()
```

The plot above shows that the full CHM extent (plotted in green) is much larger than the resulting cropped raster. Our new cropped CHM now has the same extent as the oai_boundary_tudlib object that was used as a crop extent (blue border below).

```{r}
ggplot() +
  geom_raster(data = CHM_TUD_Cropped_df,
              aes(x = x, y = y, fill = layer)) + 
  geom_sf(data = oai_boundary_tudlib, color = "blue", fill = NA) + 
  scale_fill_gradientn(name = "Canopy Height", colors = terrain.colors(10)) + 
  coord_sf()
```

We can look at the extent of all of our other objects for this field site.

```{r}
st_bbox(CHM_TUD)
```

```{r}
st_bbox(CHM_TUD_Cropped)
```

```{r}
st_bbox(oai_boundary_tudlib)
```

```{r}
leisure_locations_selection <- st_read(here("data", "delft-leisure.shp")) |> 
  filter(leisure %in% c("playground", "picnic_table"))

st_bbox(leisure_locations_selection)
```

Our plot location extent is not the largest but is larger than the AOI Boundary. It would be nice to see our vegetation plot locations plotted on top of the Canopy Height Model information.


# Extract Raster Pixels Values Using Vector Polygons

Often we want to extract values from a raster layer for particular locations - for example, plot locations that we are sampling on the ground. We can extract all pixel values within 20m of our x,y point of interest. These can then be summarized into some value of interest (e.g. mean, maximum, total).


![Image Source: National Ecological Observatory Network (NEON)](fig/BufferSquare.png)
To do this in R, we use the `extract()` function. The `extract()` function requires:

The raster that we wish to extract values from,
The vector layer containing the polygons that we wish to use as a boundary or boundaries,
we can tell it to store the output values in a data frame using `df = TRUE`. (This is optional, the default is to return a list, NOT a data frame.) .
We will begin by extracting all canopy height pixel values located within our `oai_boundary_tudlib` polygon.

```{r}
tree_height <- extract(x = CHM_TUD, y = st_as_sf(oai_boundary_tudlib), df = TRUE)

str(tree_height)
```

When we use the `extract()` function, R extracts the value for each pixel located within the boundary of the polygon being used to perform the extraction - in this case the `aoi_boundary_HARV` object (a single polygon). Here, the function extracted values from 621,642 pixels.

We can create a histogram of tree height values within the boundary to better understand the structure or height distribution of trees at our site. We will use the column `layer` from our data frame as our x values, as this column represents the tree heights for each pixel.

```{r}
ggplot() + 
  geom_histogram(data = tree_height, aes(x = layer)) +
  ggtitle("Histogram of CHM Height Values (m)") +
  xlab("Tree Height") + 
  ylab("Frequency of Pixels")
```

We can also use the `summary()` function to view descriptive statistics including min, max, and mean height values. These values help us better understand vegetation at our field site.

```{r}
summary(tree_height$layer)
```

# Summarize Extracted Raster Values

We often want to extract summary values from a raster. We can tell R the type of summary statistic we are interested in using the `fun = argument`. Let’s extract a mean height value for our AOI. Because we are extracting only a single number, we will not use the `df = TRUE` argument.

```{r}
mean_tree_height_AOI <- extract(x = CHM_TUD, y = st_as_sf(oai_boundary_tudlib), fun = mean)

head(mean_tree_height_AOI)
```

It appears that the mean height value, extracted from our LiDAR data derived canopy height model is 4.3 meters.

# Extract Data using x,y Locations

We can also extract pixel values from a raster by defining a buffer or area surrounding individual point locations using the `extract()` function. To do this we define the summary argument (`fun = mean`) and the buffer distance (`buffer = 20`) which represents the radius of a circular region around each point. By default, the units of the buffer are the same units as the data’s CRS. All pixels that are touched by the buffer region are included in the extract.

![Image Source:National Ecological Observatory Network (NEON)](BufferCircular.png)


Let’s put this into practice by figuring out the mean tree height in the 20m around the tower location (`point_Delft`). Because we are extracting only a single number, we will not use the `df = TRUE` argument.

```{r}
point_Delft <- st_read(here("data", "delft-leisure.shp"))

mean_tree_height_tower <- extract(x = CHM_TUD,
                                  y = point_Delft,
                                  buffer = 20,
                                  fun = mean)

mean_tree_height_tower
```

# Challenge: Extract Raster Height Values For Plot Locations

- Use the leisure locations object (`point_Delft`) to extract an average tree height for the area within 20m of each playground and picnic table in the study area. Because there are multiple playgrounds and picnic tables, there will be multiple averages returned, so the `df = TRUE` argument should be used.

- Create a plot showing the mean tree height of each area.


---

## Solution
```{r}
leisure_locations_selection <- st_read(here("data", "delft-leisure.shp")) |> 
  filter(leisure %in% c("playground", "picnic_table"))

# extract data at each plot location
mean_tree_height_plots_TUD <- extract(x = CHM_TUD,
                                       y = leisure_locations_selection,
                                       buffer = 20,
                                       fun = mean,
                                       df = TRUE)

# view data
head(mean_tree_height_plots_TUD)
```

```{r}
# plot data
ggplot(data = mean_tree_height_plots_TUD, aes(ID, layer)) + 
  geom_col() + 
  ggtitle("Mean Tree Height at each Plot") + 
  xlab("Plot ID") + 
  ylab("Tree Height (m)")
```

---

# Summary and keypoints.

We have seen how to crop a raster to the extent of a vector layer and how to extract values from a raster that correspond to a vector file overlay.

In short:

- Use the `crop()` function to crop a raster object.
- Use the `extract()` function to extract pixels from a raster object that fall within a particular extent boundary.