Lab4

Introduction

In this lab, we will answer the following questions:
In the Southern Rockies ecoregion,
Q1: How long does forest take to ‘recover’ post-fire?
Q2: Does time to recovery vary as a function of pre-fire beetle infestation?

1. Set-up

1.1 Create folders

1.1.1 Create a Lab4 folder nested within the folder you use for this class. Do not use spaces in the name of the folder (i.e., not “Lab 4”)
1.1.2. Create a data subfolder within your Lab4 folder

1.2 Set up your R Script

1.2.1 Open RStudio
1.2.2 Create an R Script file. Click File -> NewFile -> R Script
1.2.3 Save the R Script to your Lab4 folder. Click File > Save As

1.3 Set your working Directory and start coding

We need to set a working directory, which is the default location of any files you read into R or save out of R. We’ll do this by specifying the path name to the Lab4 folder in your new R Script.

1.3.1. Copy and past the text in the gray rectangle below into your R Script
1.3.2. MODIFY it to your own path name
1.3.3. Run the line of code by selecting anywhere in that line and clicking the Run icon. The shortcut is Control + Enter
1.3.4. If it ran correctly, you should see the line of code appear in the Console with no errors
1.3.5. Going forward in this tutorial, anything in a gray box is code that you should copy and paste into your script, modify as appropriate, and run. 1.3.6 Remember that you can use hashtages to comment out code. These are like notes to you that R skips over. 1.3.7 Save your R Script frequently.

setwd("/Users/megancattau/Dropbox/0_BoiseState/Courses/2024_spr_Ecol_Dist/Labs/Lab4_Recovery")

1.4. Install, load, and call the packages below and their dependencies.

1.4.1. Install the packages
If you do not already have these packages installed, run this line:   install.packages(“lme4”, type = “source”)

Remember that you only need to install packages once, even if you shut down your R session or start working on another Script. However, you will need to set your working directory and load the packages every time you start a new R session.

1.4.2. Load the packages

library(lme4)

Note: Anytime you close an R session, you would need to rerun your code. The script is saved, so the code is available to you, but the working directory will no longer be set, the libraries will no longer be loaded, objects you had previously created will no longer be available, etc.

1.5 Import Data

In this script, we will use data that has been preprocessed by (1) importing or generating annual 250m resolution rasters for disturbances as detected by the Aerial Detection Survey (ADS), fire (MTBS), MODIS Vegetation Continuous Fields (VCF), and the National Land Cover Dataset (NLCD) for EPA L3 Ecoregion Southern Rockies; (2) sample these data at regualrly spaced points (250m), (3) export them as a csv, and (4) further format them to capture VCF pre- and post-fire in forested areas.

The MODIS Vegetation Continuous Fields (VCF) data (MOD44B Version 6.1 data product) will be our forest response metric. They represent the percent tree cover per 250m pixel as detected by MODIS. These data can be found here: https://lpdaac.usgs.gov/products/mod44bv061/ Data Pool link: https://e4ftl01.cr.usgs.gov/MOLT/MOD44B.061/ Data are generated annually from JD 065 (early March) imagery.

1.5.1 Download the ‘samples_fire.csv’ data from our Canvas site   1.5.2 Place it in a data subfolder in your Lab4 folder   1.5.3 Import it and look at the column headings

samples<-read.csv("data/samples_fire.csv")
samples<-samples[,-1]
names(samples)

Each row in the data represents a pixel in our AOI. The attributes include the x and y coordinates, a unique ID (FID), an indicator of whether beetle infestation occurred prior to fire, a VCF variable name (you can ignore this), a VCF value, and how many years since fire that value was taken from. Note above what VCF is measuring. Note that time ‘0’ represents an average value for all years pre-fire for which data were available.

2. Analysis

We will evaluate Recovery since fire by answering the following questions:
In the Southern Rockies ecoregion,
Q1: At what age post-fire is forest ‘recovered’?
Q2: Does recovery vary as a function of pre-fire beetle infestation?

2.1 Q1: How long does forest take to ‘recover’ post-fire?

We will evaluate how many years after fire post-fire VCF resembles pre-fire VCF

2.1.1 Fit a model to the data

# Fit quadratic equation to data (excluding year 0)
data<-samples[samples$time!=0,]
glm.fit<-glm(VCF~ poly(time,2), data=data)
summary(glm.fit)

# Predict recovery over 1-32 years
time_pred=seq(1, 32, len=10000)
VCF_pred = predict(glm.fit, newdata=data.frame(time=time_pred))
predictions_VCF<-data.frame(time_pred, VCF_pred)

2.1.2 Determine age of ‘recovery’

# Calculate mean and sd pre-fire and post-fire VCF over time
VCF_post_fire_mean<-aggregate(samples$VCF, by=list(samples$time), FUN=mean, na.rm=TRUE)
VCF_post_fire_sd<-aggregate(samples$VCF, by=list(samples$time), FUN=sd, na.rm=TRUE)
VCF_post_fire_plus_sd<-cbind(VCF_post_fire_mean$Group.1, (VCF_post_fire_mean$x+VCF_post_fire_sd$x))
VCF_post_fire_minus_sd<-cbind(VCF_post_fire_mean$Group.1, (VCF_post_fire_mean$x-VCF_post_fire_sd$x))

# Determine pre-fire VCF mean and sd
VCF_post_fire_mean[1,2]
VCF_post_fire_sd[1,2]
# pre-fire mean is VCF 33.8 +/- 13.4

# Solve for pre-fire value, meaning and what timestep are post-fire mean VCF values predicted to be equivalent to pre-fire mean values?
recovery = mean(predictions_VCF[predictions_VCF$VCF_pred>33.7&predictions_VCF$VCF_pred<33.8,]$time_pred)
recovery
# Forest "recovers" 29.3 years after fire

2.1.3 Visualize ‘recovery’

plot(VCF_post_fire_mean, ylim=range(0:50), xlab="Time since fire (years)", ylab="Percent tree cover (VCF)", pch=20)
lines(VCF_post_fire_plus_sd, lty=3, col="darkgrey")
lines(VCF_post_fire_minus_sd, lty=3, col="darkgrey")
lines(time_pred, VCF_pred, lty=1, col="lightgreen")
abline(h=VCF_post_fire_mean[1,2], col="darkgreen", lty=4)

legend(2, 50, legend=c("Pre-fire forest state", "Recovery"), lty=c(4,1), lwd=c(2.5, 2.5), col=c("darkgreen", "lightgreen"))
text(recovery-5, VCF_post_fire_mean[1,2]+3, labels=c("Recovery=29.3 years"))

2.1.4 Sidenote: Perhaps you are interested in the period of lagged mortality you see here and want to know when does forest start to recover post-fire.

# get the minimum predicted value for VCF
min(VCF_pred)
# determine at what timestep that minimum value occurs
end_mortality = predictions_VCF[predictions_VCF$VCF_pred==min(predictions_VCF$VCF_pred),]$time_pred
end_mortality
# Forest mortality ends 9.6 years after fire 

2.2 Q2: Does time to recovery vary as a function of pre-fire beetle infestation?

Does recovery vary as a function of pre-fire infestation? We will evaluate how many years after fire post-fire VCF resembles pre-fire VCF for areas that have experienced a previous infestation versus those that have not.

2.2.1 Fit a model to the data

# Fit quadratic equation to data (excluding year 0). VCF ~ time+ previous infestation
data$infest_bin<-factor(data$infest_bin)
fit_fire_beetle<-lmer(VCF~ poly(time,2)+(1|infest_bin), data=data)
summary(fit_fire_beetle)

# Predict recovery over 1-32 years ~ infest or not
newdata_infest=data.frame(time=time_pred, infest_bin=factor(1))
newdata_noinfest=data.frame(time=time_pred, infest_bin=factor(0))
infest_pred = predict(fit_fire_beetle, newdata=newdata_infest)
noinfest_pred = predict(fit_fire_beetle, newdata=newdata_noinfest)
predictions_infest<-data.frame(time_pred, infest_pred)
predictions_noinfest<-data.frame(time_pred, noinfest_pred)

2.2.2 Determine age of ‘recovery’

# subset data by infested and not
VCF_infest_infest<-samples[samples$infest_bin==1,]
VCF_infest_noinfest<-samples[samples$infest_bin==0,]

# Calculate mean and sd pre-fire and post-fire VCF over time
VCF_post_fire_infest_mean<-aggregate(VCF_infest_infest$VCF, by=list(VCF_infest_infest$time), FUN=mean, na.rm=TRUE)
VCF_post_fire_infest_sd<-aggregate(VCF_infest_infest$VCF, by=list(VCF_infest_infest$time), FUN=sd, na.rm=TRUE)
VCF_post_fire_infest_plus_sd<-cbind(VCF_post_fire_infest_mean$Group.1, (VCF_post_fire_infest_mean$x+VCF_post_fire_infest_sd$x))
VCF_post_fire_infest_minus_sd<-cbind(VCF_post_fire_infest_mean$Group.1, (VCF_post_fire_infest_mean$x-VCF_post_fire_infest_sd$x))

VCF_post_fire_noinfest_mean<-aggregate(VCF_infest_noinfest$VCF, by=list(VCF_infest_noinfest$time), FUN=mean, na.rm=TRUE)
VCF_post_fire_noinfest_sd<-aggregate(VCF_infest_noinfest$VCF, by=list(VCF_infest_noinfest$time), FUN=sd, na.rm=TRUE)
VCF_post_fire_noinfest_plus_sd<-cbind(VCF_post_fire_noinfest_mean$Group.1, (VCF_post_fire_noinfest_mean$x+VCF_post_fire_noinfest_sd$x))
VCF_post_fire_noinfest_minus_sd<-cbind(VCF_post_fire_noinfest_mean$Group.1, (VCF_post_fire_noinfest_mean$x-VCF_post_fire_noinfest_sd$x))

# Determine pre-fire VCF mean and sd
VCF_post_fire_infest_mean[1,2]
# For areas with prior infestation, VCF pre-fire mean is 38.4
VCF_post_fire_noinfest_mean[1,2]
# For areas without prior infestation, VCF pre-fire mean is 33.0

# Solve for pre-fire value, meaning and what timestep are post-fire mean VCF values predicted to be equivalent to pre-fire mean values?
recovery_infest = mean(predictions_infest[predictions_infest$infest_pred>38.4 & predictions_infest$infest_pred<38.5,]$time_pred, na.rm=TRUE)
recovery_infest
# Forest "recovers" 31.8 years after fire for areas with previous infestation

recovery_noinfest = mean(predictions_noinfest[predictions_noinfest$noinfest_pred>32.9 & predictions_noinfest$noinfest_pred<33,]$time_pred, na.rm=TRUE)
recovery_noinfest
# Forest "recovers" 29.3 years after fire for areas without previous infestation

2.2.3 Visualize ‘recovery’

plot(VCF_post_fire_noinfest_mean, ylim=range(0:50), xlab="Time since fire (years)", ylab="Percent woody vegetation (VCF)", pch=20, col="green")
lines(VCF_post_fire_noinfest_plus_sd, lty=3, col="green")
lines(VCF_post_fire_noinfest_minus_sd, lty=3, col="green")
abline(h=VCF_post_fire_noinfest_mean[1,2], lty=4, col="green")

points(VCF_post_fire_infest_mean, ylim=range(0:50), xlab="Time since fire (years)", ylab="Percent tree cover (VCF)", pch=20, col="red")
lines(VCF_post_fire_infest_plus_sd, lty=3, col="red")
lines(VCF_post_fire_infest_minus_sd, lty=3, col="red")
abline(h=VCF_post_fire_infest_mean[1,2], lty=4, col="red")
legend(2, 50, legend=c("Beetle infestation prior to fire", "No infestation prior to fire"), lty=c(1,1), lwd=c(2.5, 2.5), col=c("red", "green"))

lines(time_pred, infest_pred, lty=1, col="red")
lines(time_pred, noinfest_pred, lty=1, col="lightgreen")

text(recovery_infest-5, VCF_post_fire_infest_mean[1,2]+3, labels=c("Recovery=31.8 years"), col="red")
text(recovery_noinfest-5, VCF_post_fire_noinfest_mean[1,2]+3, labels=c("Recovery=28.8 years"), col="green")

Report

For your lab report this week, please include
1. Each of the two figures we generated along with a brief, simple description. Pretend you are explaining the figures to someone with no context and include anything they would need to know to understand what they are depicting. At a minimum, please include the following for each of the figures:
a. What data went into making this figure?
b. What are the ‘dots’ on the figures?
b. What are the dashed horizontal lines?
c. What are the solid, curved lines?
d. What are the dashed lines on either side of the curved lines?
e. What does it mean when the solid, curved line crosses the dashed horizontal line?

Please also answer the following questions. There is no one right answer for any of these, and you’ll receive full points for thinking critically

2. What assumptions or biases are there in the data we used, and how do you think that might effect the results? You can speak to any of the data we used: the Southern Rockies ecoregion AOI, MTBS fire data, ADS beetle damage data, and / or MODIS VCF percent tree cover data.
   
3. How are we defining ‘recovery’ with this analysis? Does this seem like a reasonable approach?
   
4. For figure 2, why does recovery occur approximately 3 years sooner for areas with no prior beetle infestation relative to areas with prior beetle infestation even though the recovery curves appear similar?
   
5. What alternative methods might you use to define ‘recovery’?
   

Bonus: Generate any other figure or statistic that piques your interest. For example, you may want to redefine ‘recovery’ for Q2, or try out different forms of the model for Qs 1 or 2