---
title: "Hessami et al., 2024_code"
author: "Mateen Hessami"
date: '2024-10-11'
output: github_document
---


```{r, engine = "Figure", fig.cap = "Figure 1"}
#Code and data (RV_Cutblock_Data) for Figure 1 to calculate forested area harvested in the Lake Revelstoke Valley annually. Map was created in QGIS. 

####Load Packages####
library(tidyverse)
library(rgeos)
library(rgdal)
library(ggplot2)

#### Load in Revelstoke Valley cutblock dataframe ####
#Call in RV_Cutblock_Data (ALL_cutblocks_LRV_25Sept_2022)
clip <- st_read ("ALL_cutblocks_LRV_25Sept_2022.shp")

clip <-unite(clip,"Harvest_Ye",c(DN,Harvest_Ye,HARVESTYR),remove=TRUE,na.rm=TRUE) 
%>%
mutate(Harvest_Ye=as.numeric(Harvest_Ye)) 
 

######### Calculate annual area cut #######
clip$areacalc <- as.numeric(st_area(clip))

#drop the unesecary columns - fine tune focus (hectares) 
clip %>% st_drop_geometry() %>% group_by(Harvest_Ye) %>% summarize(sum(areacalc)) 

#summarize data
clip2_df <- clip %>% st_drop_geometry() %>% dplyr::group_by(Harvest_Ye) %>% dplyr::summarize(sumarea=sum(areacalc)) %>%
arrange(Harvest_Ye) 

#reduce to just single year, organize sequentially, and arrange as one observation per year##
clip3_df <- clip2_df %>% dplyr::group_by(Harvest_Ye) %>% dplyr::summarize(sumarea=sum(sumarea/10000)) %>% arrange(Harvest_Ye)


###adding cumulative cut area per year, “sumarea”###
clip3_df <- clip3_df %>% mutate(cumarea=cumsum(sumarea))

####### Create Plot ######
#Plot prep
clip3_df_long <- clip3_df %>% tidyr::gather(key=variable,value=value,c(2,)) %>% 
  arrange(Harvest_Ye,desc(variable))

Trans <- ggplot(clip3_df_long,
           aes(x = Harvest_Ye, y=value)) +
  geom_bar(stat="identity",fill="dark green", width = .5)+
  labs(x="year", y= "Forested area harvested (ha)") +
 theme(text=element_text(size = 18),
       legend.position="right",
         panel.background = element_rect(fill = "transparent",
                                  colour = NA_character_),
       plot.background = element_rect(fill = "transparent", color = NA),
         panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
         panel.border=element_rect(color="black", fill=NA,size=.6))+

   scale_x_continuous(name="Year",labels = waiver(), breaks= c(1961,1970,1980,1990,2000,2010,2020)) +
   scale_y_continuous(expand = c(0,0))

```


```{r, engine = "Figure", fig.cap = "Figure 2, Panel A"}
#Code and data (Used_Avail_Data.csv) for moose habitat selection analysis, derived from Resource Selection Functions and a Generalized Additive Mixed Model (Figure 2, Panel A).  

####Load Packages####
library(ggplot2)
library(getViz)
library(gamm4)
library(qgam)

# Calling in cleaned used and available data"
Used_Avail_Data <- read.csv("Used_Avail_Data.csv")

###### Run GAMM and plot ##### 
MooseGamm<- gamm(OBSERVED~s(cutage3,k=1), random=list(burst2=~1), family = "binomial",data = Used_Avail_Data)
plot(MooseGAMM$gam,shade=TRUE,lwd=2)
abline(h=0)
abline(v=14)
abline(v=36)


```



```{Figure 2, Panel B}
#Code and data (Annual_Hectares_Cutblocks.csv) Plot 2, Panel B hectares of optimal (e.g., carrying capacity) moose habitat from 1992-2020.  

####Load Packages####
library(roll)
library(RcppRoll)
library(Rcpp)
library (ggplot2)

#Read in dataframe that includes annual hectares cut in the study area 
MooseKTrend <- read.csv("Annual_Hectares_Cutblocks.csv")

#Calculate the total area of habitat in 3 categories (recent, optimal, old)
MooseKTrend $recent <- roll_sum(MooseKTrend$sumarea, n = 14,fill=NA, align = "right")
MooseKTrend$"Landscape K" <- roll_sum(MooseKTrend$sumarea, n = 22, fill=NA, align = "right")
MooseKTrend$old <- roll_sum(MooseKTrend$sumarea, n= 36, fill=NA, align = "right")
MooseKTrend$cumarea <- cumsum(MooseKTrend$sumarea)
#View(MooseKTrend)

#defining windows for recent, optimal, and old moose habitat 
MooseKTrend$recent[1:14] <- MooseKTrend$cumarea[1:14] 
MooseKTrend$"Landscape K"[14:36] <- MooseKTrend$cumarea[14:36] 
MooseKTrend$old[36:61] <- MooseKTrend$cumarea[36:60] 
MooseKTrend$cumarea <- NULL

#rename year to Harvest_Ye
library(dplyr)
MooseKTrend2 <- MooseKTrend %>% dplyr::rename(year = Harvest_Ye)

#clean up MooseKTrend2 for plot
library(dplyr)
library(tidyr)
MooseKTrend2 <- MooseKTrend2 %>% dplyr::select(-sumarea) %>% replace_na(list("Landscape K"=0)) %>% replace_na(list(old=0)) %>%  pivot_longer(-year, names_to = "variable", values_to = "value")

#Plot
filter(MooseKTrend2,year>=2003 & year<=2020 & variable =="Landscape K") %>% 
ggplot(aes(x=year, y=value,width=0.5)) +
  geom_line()+ 
  labs(x="year", y= "Moose K (hectares)") +
  theme(text=element_text(size = 12),
        panel.background = element_blank(),
        panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
        panel.border=element_rect(color="black", fill=NA,size=.6))+
  scale_x_continuous(breaks = seq(1993,2020,by=2))+
  scale_y_continuous(breaks = seq(0,60000,by=2000))

```


```{Figure 3, Panel A}
#Code and data (Hunter_Harvest.csv) for plot 3, panel A, where observed and theoretical (2.5% and 6%) cumulative moose harvest statistics are plotted. 

####Load Packages####
library(readr)
library(dplyr)
library(lubridate)
library(readxl)
library(tidyverse)
library(ggplot2)

###Read in hunter survey dataframe### 
Hunter_Harvest <- read.csv("Hunter_Harvest.csv")

###Summarize total kills###
Hunter_Harvest <- group_by(Hunter_Harvest,Year) %>% 
  dplyr::summarize(TotalKillsLRV=sum(Kills)) %>% 
  dplyr::filter(Year>=2003 & Year<=2020) 


###Now weave in hunter “kill statistics” into Hunter_Harvest - small table gets joined second###
Hunter_Harvest2 <- left_join(HunterHarvest,LRVkillStats,by=c('year'='Year')) %>% 
    dplyr::filter(year>=2003 & year<=2020) 

####Create new column. 0.065=moose per hectare at carrying capacity###
Hunter_Harvest2$mooseNproj <- Hunter_Harvest2$LRV_K*0.065
                                                  
###Create the 3 scenarios - (1) observed - cumulative total kills. (2) 2.5%. (3) 6%. For 2.5% and 6% scenario apply calculation to theoretical harvest###

Hunter_Harvest2$CUMkillsobs <- ave(Hunter_Harvest2$TotalKillsLRV, FUN = cumsum)
Hunter_Harvest2$CUMkillsTheoTwo <- ave(Hunter_Harvest2$mooseNproj*0.025, FUN = cumsum)                        
Hunter_Harvest2$CUMkillsSix <- ave(Hunter_Harvest2$mooseNproj*0.06, FUN = cumsum)   

###changing Y axis to read zero #####
Year2002 <- Hunter_Harvest2[1,]
Year2002[,3:19] <- 0
Year2002[1,2] <- 2002
#names
# CUMkillsobs
# CUMkillsTheoTwo               
# CUMkillsSix

####Plotting####
library(ggplot2)
bind_rows(Year2002,Hunter_Harvest2) %>% 
#Hunter_Harvest2 %>% 
  mutate(CUMkillsTheoTwo=CUMkillsTheoTwo,CUMkillsobs=CUMkillsobs,CUMkillsSix=CUMkillsSix)%>% 
  dplyr::select(year,CUMkillsobs,CUMkillsTheoTwo,CUMkillsSix) %>% 
  pivot_longer(c(CUMkillsTheoTwo,CUMkillsobs,CUMkillsSix),names_to = "Cumulative Kills", values_to = "value") %>% 
ggplot(aes(x=year,y=value,col=`Cumulative Kills`,linetype=`Cumulative Kills`))+
geom_line(se=FALSE)+
  xlab("Year")+
  ylab('Cumalative Moose Kills')+
  scale_x_continuous(name="Year",breaks= c(2001,2003,2005,2007,2009,2011,2013,2015,2017,2019,2021),expand=c(0,0))+
  scale_y_continuous(name="Cumaltive Moose Harvested",breaks= c(0,100,200,300,400,500,600,700,800,900,1000,1100,1200,1300,1400,1500),expand=c(0,0))+
# scale_fill_viridis_d(option = "viridis", end = 0.7) +
#   scale_color_viridis_d(option = "viridis", end = 0.7)+
  scale_colour_manual(name="Legend",values=c("CUMkillsobs"="blue","CUMkillsTheoTwo"="red","CUMkillsSix"="green"))+
   labs(color="legend")+
  theme_classic())
  #scale_linetype_manual(values=c("solid","dashed"))
  # theme(axis.text.x=element_blank(),
  #  axis.ticks.x=element_blank())




```



```{Figure 3, Panel B}
#code and data (MasterFigure.csv) for plot for Figure 3 Panel B; changes in moose abundance (N) (dotted black line) based on aerial surveys. The red line is the ratio of calf to adult moose observed in aerial surveys. The light gray dotted trend line is Kill Per Unit Effort, an index of moose hunter harvest success. 

####Load Packages####
library(tidyverse)
library(ggplot2)

###Read in master database with KPUE, calf:cow ratios, and abundance estimates###
MasterFigure <- read.csv("MasterFigure.csv")

####Plotting#### 
MasterFigure %>% dplyr::select(Year,SRB,cc_ratio) %>% 
  mutate(cc_ratio=cc_ratio*20) %>% drop_na() %>% 
  pivot_longer(-Year,names_to = "method", values_to = "values") %>% 
  ggplot()+ 
  theme_classic()+
  geom_point(aes(x=Year,y=values,col=method, shape=method))+ 
   geom_line(aes(x=Year,y=values,col=method,shape=method,linetype=method))+
   geom_point(data=KPUE,mapping=aes(x=Year,y=KPUEperYear,color="KPUE",shape="KPUE"))+
  geom_line(data=KPUE,mapping=aes(x=Year,y=KPUEperYear,color="KPUE"), linetype="dotted")+
  scale_y_continuous(
    "Population estimate", 
    sec.axis = sec_axis(~ . /20, name = "Calves per 100 Adult Females"),
  expand = c(0,0))+
  scale_colour_manual(values=c("SRB"="black","KPUE"="#646464","cc_ratio"= "red"))+
  #scale_linetype_manual(values=c("SRB"="solid","KPUE"="dotted","cc_ratio"= "dashed"))+
   scale_x_continuous(breaks = seq(1994,2020,by=2))+
  expand_limits(y=0)


```


```{Figure 4, Panel B}
#Code and data (all within code below) for figure 4; recruitment (calves per 100 adult females) as a function of per capita resource availability. The positive effect of resource availability on cow:calf ratios support a density dependent response of moose  

####Load Packages####
library(tidyverse)
library(ggplot2)

###Read in master database with KPUE, calf:cow ratios, and abundance estimates###
MasterFigure <- read.csv("MasterFigure.csv")

####Plotting#### 
ggplot() +
  geom_point(data = dplyr::filter(MasterFigure, Year < 2018),
             mapping = aes(x = hpm, y = cc, label = Year),
             color = "blue", size = 2) + # Increase size for blue dots
  geom_point(data = dplyr::filter(MasterFigure, Year >= 2018),
             mapping = aes(x = hpm, y = cc, label = Year),
             color = "red", shape = 15, size = 2) + # Increase size for red squares
  geom_smooth(method = lm, se = FALSE, data = dplyr::filter(MasterFigure, Year < 2018),
              mapping = aes(x = hpm, y = cc),
              color = "black", size = 0.5) +
  theme_classic() +
  expand_limits(x = 0, y = 0) +
  coord_fixed() +
  labs(x = "Hectares Per Moose", y = "Recruitment (calves:100 cows)") +
  theme(legend.title = element_blank())

```



```{Figure 5, Panel B}
#Projected moose density based on moose carrying capacity under different forest harvest scenarios. The dashed line is the threshold density of moose in which wolf densities would be low enough to satisfy caribou recovery projections (Figure 5, both panels).   

####Load Packages###
library(roll)
library(RcppRoll)
library(Rcpp)
###Load in template dataframe (hectares cut per year from study area###
Template <-read.csv("RV_CutData.csv")
##Change sumarea to 640, the average area cut per year between 2010-2020###
Template[which(Template $Harvest_Ye==2020):nrow(Template),"sumarea"]=640* 0.25
Projections <- Template %>% mutate(cumarea=cumsum(sumarea))

####Calculuate rolling summaries of recent, landscape K (e.g., optimal), and old age classifications of forests###
Projections$recent <- roll_sum(Projections$sumarea, n = 14, fill=NA, align = "right")
Projections$"Landscape K" <- roll_sum(Projections$sumarea, n = 22, fill=NA, align = "right")
Projections$old <- roll_sum(Projections$sumarea, n = 36, fill=NA, align = "right")
Projections$cumarea <- cumsum(Projections$sumarea)
Projections$recent[1:14] <- Projections$cumarea[1:14] 
Projections$"Landscape K" [14:36] <- Projections$cumarea[14:36] 
Projections$old[36:81] <- Projections$cumarea[36:81] 
Projections$cumarea <- NULL

###Plot Prep####
Projections2 <- Projections %>% dplyr::rename(year = Harvest_Ye)
SF3 <- Projections2 %>% dplyr::select(-sumarea) %>% replace_na(list(optimal=0)) %>% replace_na(list(old=0)) %>% pivot_longer(-year, names_to = "scenario", values_to = "hectares")

Projections2$SF <- "75% reduction"

##Plotting####
dplyr::filter(SF3,year>=2020 & year<=2040 &Projections=='25% reduction'
 & scenario=="Landscape K") %>% 
ggplot(aes(x=year, y=`hectares`, fill=SF, width=0.5)) +
  geom_bar(stat='identity', position='stack')+ 
  labs(x="year", y= "Hectares") +
  theme(text=element_text(size = 12),
        
        
        panel.background = element_blank(),
        panel.grid.major = element_blank(),panel.grid.minor = element_blank(),
        panel.border=element_rect(color="black", fill=NA,size=.6))+
  scale_y_continuous(expand = c(0,0))

```

















