# Overview
Phenology is the timing of life history events. Important phenological events for plants involve the growth of leaves, flowering, and senescence (death of leaves). Plants species adapt the timing of these events to local climate conditions to ensure successful reproduction. Subsequently, animal species often adapt their phenology to take advantage of food availability. As the climate shifts this synchronization is being thrown out of whack. Shifts in phenology are therefore a common yardstick of understanding how and if ecosystems are adjusting to climate change.\

Plant species may employ the following phenological strategies:\

-   winter deciduous: lose leaves in the winter, grow new leaves in the spring\
-   drought deciduous: lose leaves in the summer when water is limited\
-   evergreen: maintain leaves yearround\

credit: this lab is based on a materials developed by Chris Kibler.

# Task
In this lab we are analyzing plant phenology near the Santa Clara River which flows from Santa Clarita to Ventura. We will investigate the phenology of the following plant communities:\

-   riparian forests: grow along the river, dominated by winter deciduous cottonwood and willow trees\
-   grasslands: grow in openspaces, dominated by drought deciduous grasses\
-   chaparral shrublands: grow in more arid habitats, dominated by evergreen shrubs\

To investigate the phenology of these plant communities we will a time series of Landsat imagery and polygons identifying the locations of study sites within each plant community.\

## Data

**Landsat Operational Land Imager (OLI sensor)**\

-   8 pre-processed scenes\
    -   Level 2 surface reflectance products\
    -   erroneous values were set to NA\
    -   scale factor set to 100\
    -   bands 2-7\
    -   dates in filenname\

**Study sites**

-   polygons representing sites\
    - study_site: character string with plant type
    
## Summary of approach
-   Convert spectral relfectance into a measure of vegetation productivity (NDVI)
-   Calulate NDVI throughout the year
-   Summarize NDVI values within vegetation communities
-   Visualize changes in NDVI within vegetation communities

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

```{r include=TRUE, message=FALSE}
library(terra)
library(sf)
library(dplyr)
library(tidyr)
library(stringr)
library(ggplot2)
library(here)
library(tmap)
library(cowplot)

rm(list = ls())
here::i_am("labs/week8.Rmd")
```

### create NDVI function
Let's start by defining a function to compute the Normalized Difference Vegetation Index (NDVI). NDVI computes the difference in reflectance in the near infrared and red bands, normalized by their sum. 

```{r ndvi-function}

}
```

### computing NDVI for a single scene
We have 8 scenes collected by the Landsat OLI sensor on 8 different days throughout the year.
Let's start by loading in the first scene collected on June 12, 2018
```{r}

```

Now let's update the names of the layers to match the spectral bands they correspond to
```{r}

```

Now we can apply the NDVI function we created to compute NDVI for this scene using the `lapp()` function. The `lapp()` function applies a function to each cell using layers as arguments. Therefore, we need to tell `lapp()` which layers (or bands) to pass into the function. The NIR band is the 4th layer and the red band is the 3rd layer in our raster. In this case, because we defined the NIR band as the first argument and the red band as the second argument in our function, we tell `lapp()` to use the 4th layer first and 3rd layer second. 
```{r}

```


### computing NDVI for all scenes - attempt 2

The first attempt was pretty clunky and required a lot of copy/pasting. Because we're performing the same operations over and over again, this is a good opportunity to generalize our workflow into a function!

#### starting from scratch
Let's start over and see how we could do this more efficiently. We'll clear our environment and redefine our function for NDVI.
```{r}
rm(list = ls())
here::i_am("labs/week8.Rmd")

ndvi_fun = function(nir, red){
  (nir - red) / (nir + red)
}
```

#### outlining our function
Let's first sketch out what operations we want to perform so we can figure out what our function needs.
```{r eval=FALSE}

# note: this code is not meant to run! we're just outlining the function we want to create


# what do we notice as what we need to pass into our function?

```

#### creating a list of files/scenes
We want a list of the scenes so that we can tell our function to compute NDVI for each. To do that we look in our data folder for the relevant file.

```{r}
# here we're asking for the names of all the files in the week8 folder
# the "pattern" option returns names that end in .tif (the file extension for the landsat scences)
# the "full.names" option returns the full file path for each scene


```

#### updating our function
Now let's update our function to work with list of file names we created
```{r}

# now we're passing our function a number that will correspond to the index in the list of file names


# let's test our function by asking it to read in the first file 

```

#### run function on all scenes!
Now we can use our function to create a NDVI layer for each scene and stack them into a single rasterstack.
```{r}
# create NDVI layer for each scene (day) and stack into a single rasterstack


# update layer names to match date

```


### Compare NDVI across vegetation communities

Now that we have computed NDVI for each of our scenes (days) we want to compare changes in NDVI values across different vegetation communities.

#### read in study sites
First, we'll read in a shapefile of study sites.
```{r include=TRUE}


# plot study sites on a single NDVI layer

```

#### extract NDVI at study sites
Here we find the average NDVI within each study site. The output of extract is a data frame with rows that match the study site dataset, so we bind the results to the original dataset. 
```{r include=TRUE}

```

#### clean results
We're done! Except our data is very untidy... Let's tidy it up!\

-   convert to data frame
-   turn from wide to long format
-   turn layer names into date format
```{r include=TRUE}

```

#### plot results
```{r include=TRUE}

```

