In bcgov/forestDroughtTool: Predict Forest Drought Hazard In British Columbia

knitr::opts_chunk$set(echo = FALSE,message=FALSE)

# Load libraries
library(tidyverse)
library(forestDroughtTool)
library(lubridate)
library(magrittr)

# Load some functions
# source("C:/Users/hgriesba/Documents/git/forestDroughtTool/R/asmrCalc.R")
# source("C:/Users/hgriesba/Documents/git/forestDroughtTool/R/yearSelect.R")


# Load some functions
source("C:/Users/hgriesba/Documents/git/bcgov/forestDroughtTool/R/asmrCalc.R")
source("C:/Users/hgriesba/Documents/git/bcgov/forestDroughtTool/R/yearSelect.R")

# imputeTS has now changed so have to recode function

cleanData<-function(climData) {

  # Omit years with 10 or more consecutive missing values in any climate variable

  climData %>%
  # filter(stn==stnID) %>% 
  group_by(year) %>%
  summarise(pptNA=imputeTS::statsNA(ppt,print_only=FALSE)$longest_na_gap,
            tmxNA=imputeTS::statsNA(tmx,print_only=FALSE)$longest_na_gap,
            tmnNA=imputeTS::statsNA(tmn,print_only=FALSE)$longest_na_gap)%>%
  mutate(pptNA, pptNA = ifelse(is.na(pptNA), 0, pptNA)) %>%
  mutate(tmxNA, tmxNA = ifelse(is.na(tmxNA), 0, tmxNA)) %>%
  mutate(tmnNA, tmnNA = ifelse(is.na(tmnNA), 0, tmnNA)) %>%
  filter(pptNA<10&tmxNA<10&tmnNA<10) %>%
  select(year) %>%

  # Impute NA values using adjacent values
  left_join(climData,by="year") %>%
  mutate(ppt2=imputeTS::na_locf(ppt),ppt3=imputeTS::na_locf(ppt,option="nocb")) %>%
  mutate(tmx2=imputeTS::na_locf(tmx),tmx3=imputeTS::na_locf(tmx,option="nocb")) %>%
  mutate(tmn2=imputeTS::na_locf(tmn),tmn3=imputeTS::na_locf(tmn,option="nocb")) %>%
  rowwise() %>%
  mutate(ppt_filled=mean(ppt2,ppt3)) %>%
  mutate(tmx_filled=mean(tmx2,tmx3)) %>%
  mutate(tmn_filled=mean(tmn2,tmn3)) %>%
  dplyr::select(year,month,day,ppt_filled,tmx_filled,tmn_filled) %>% 

     return()
}

# This won't run, takes a long time


# Compile climate data into a list
X<-list(

    apex<-
      readxl::read_excel("ApexRoadside_daily.xlsx") %>% 
      mutate(year=year(Day),
             month=month(Day),
             day=day(Day)) %>% 
      rename(tmn="Tmin",
             tmx="Tmax",
             ppt="Precip") %>% 
     mutate(stn="apex") %>% 
      dplyr::select(stn,year,month,day,tmn,tmx,ppt) %>% 
      mutate_at(vars(year:ppt),as.numeric), 


    beaverdell<-
      readxl::read_excel("Beaverdell_daily.xlsx") %>% 
      mutate(year=year(Day),
             month=month(Day),
             day=day(Day)) %>% 
      rename(tmn="Tmin",
             tmx="Tmax",
             ppt="Precip") %>% 
    mutate(stn="beaverdell") %>% 
      dplyr::select(stn,year,month,day,tmn,tmx,ppt)%>% 
      mutate_at(vars(year:ppt),as.numeric),

    bigCreek<-
      readxl::read_excel("BigCreek_1893-1998.xlsx",range="B2:J38686",skip=1,col_types = "numeric") %>% 
      transmute(stn="bigCreek",
                year=Year,
             month=Month,
             day=Day,
             tmn=Tmin,
             tmx=Tmax,
             ppt=Total)%>% 
      mutate_at(vars(year:ppt),as.numeric),

    bmn<-
      readxl::read_excel("BMN_1993-2020.xlsx",skip=1) %>% 
      setNames(c("time","ppt","tmx","tmn")) %>% 
      mutate(stn="bmn",
             year=year(time),
             month=month(time),
             day=day(time)) %>% 
      dplyr::select(stn,year,month,day,tmn,tmx,ppt)%>% 
      mutate_at(vars(year:ppt),as.numeric),


    fording<-
      readxl::read_excel("FordingRiverCominco_1970-2017.xlsx",skip=0) %>% 
        transmute(stn="fording",
                year=Year,
             month=Month,
             day=Day,
             tmn=Tmin,
             tmx=Tmax,
             ppt=Total)%>% 
      mutate_at(vars(year:ppt),as.numeric),

    greenstone<-
      readxl::read_excel("GreenstoneHub_daily.xlsx") %>%   
      mutate(year=year(Day),
             month=month(Day),
             day=day(Day)) %>% 
      rename(tmn="Tmin",
             tmx="Tmax",
             ppt="Precip") %>% 
    mutate(stn="greenstone") %>% 
      dplyr::select(stn,year,month,day,tmn,tmx,ppt)%>% 
      mutate_at(vars(year:ppt),as.numeric),

    idaBell<-
      readxl::read_excel("IdaBell3_daily.xlsx") %>%   
      mutate(year=year(Day),
             month=month(Day),
             day=day(Day)) %>% 
      rename(tmn="Tmin",
             tmx="Tmax",
             ppt="Precip") %>% 
    mutate(stn="idaBell") %>% 
      dplyr::select(stn,year,month,day,tmn,tmx,ppt)%>% 
      mutate_at(vars(year:ppt),as.numeric),

    jellicoe<-
      readxl::read_excel("JellicoeDaily_1995-2020.xlsx") %>% 
        transmute(stn="jellicoe",
                year=Year,
             month=Month,
             day=Day,
             tmn=Tmin,
             tmx=Tmax,
             ppt=Total)%>% 
      mutate_at(vars(year:ppt),as.numeric),

    kimberley<-
      readxl::read_excel("KimberleyPCC_1975-2020.xlsx") %>% 
        transmute(stn="kimberley",
                year=Year,
             month=Month,
             day=Day,
             tmn=Tmin,
             tmx=Tmax,
             ppt=Total)%>% 
      mutate_at(vars(year:ppt),as.numeric),

    pennask <- 
      readxl::read_excel("PennaskSummit_daily.xlsx") %>%   
      mutate(year=year(Day),
             month=month(Day),
             day=day(Day)) %>% 
      rename(tmn="Tmin",
             tmx="Tmax",
             ppt="Precip") %>% 
    mutate(stn="pennask") %>% 
      dplyr::select(stn,year,month,day,tmn,tmx,ppt)%>% 
      mutate_at(vars(year:ppt),as.numeric),

    redLake<-
        readxl::read_excel("RedLakeDaily_1974-2020.xlsx") %>% 
        transmute(stn="redLake",
                year=Year,
             month=Month,
             day=Day,
             tmn=Tmin,
             tmx=Tmax,
             ppt=Total)%>% 
      mutate_at(vars(year:ppt),as.numeric),


marysville<-
  readxl::read_excel("Marysville_IDFdm2BEC7_cleaned.xlsx") %>% 
  pivot_longer(cols=contains("Day"),values_to="value",names_to="day") %>% 
  pivot_wider(names_from="Element",values_from="value") %>% 
  rename(ppt="12",
         tmx="1",
         tmn="2") %>% 
  mutate(day=str_remove(day,pattern="Day ")) %>% 
  mutate(day=as.integer(day)) %>%
  mutate(stn="marysville") %>% 
  dplyr::select(stn,year=Year,month=Month,day,tmn,tmx,ppt),

    sparwood<-
        readxl::read_excel("Sparwood_1980-2020.xlsx") %>% 
        transmute(stn="sparwood",
                year=Year,
             month=Month,
             day=Day,
             tmn=Tmin,
             tmx=Tmax,
             ppt=Total)%>% 
      mutate_at(vars(year:ppt),as.numeric),

    wasa<-
       readxl::read_excel("Wasa_1924-2017.xlsx",skip=1,col_types = "numeric") %>% 
      transmute(stn="wasa",
                year=Year,
             month=Month,
             day=Day,
             tmn=Tmin...6,
             tmx=Tmax...5,
             ppt=Total...10)%>% 
      mutate_at(vars(year:ppt),as.numeric)

) # close list

# now rbind the data frames
rawData<-
  bind_rows(X) %>% 
  mutate(stn=factor(stn)) %>% 
  ungroup()

# try cleanData a second way
stnList=unique(rawData$stn)

for (i in 1:length(stnList)) {

  if (i==1) {cleanedData=data.frame()}

  cleanedData<-
    rawData %>% 
    filter(stn==stnList[i]) %>% 
    cleanData() %>% 
    mutate(stn=stnList[i]) %>% 
    dplyr::select(stn,everything()) %>% 
    rbind(cleanedData)
  }


climData<-
  list(rawData=rawData,cleanedData=cleanedData)

save(climData,file="climData.RData")

Here is table summarizing the number of years for each station that have sufficient climate data for AET/PET calculation, as per DeLong et al. (2019):

# Load climate data (compiled from previous step)
load("climData.RData")

# Summarize data
climData %>% 
  pluck("cleanedData") %>% 
  rename(Station="stn") %>% 
  mutate(Station=str_to_title(Station)) %>% 
  group_by(Station) %>% 
  summarise(NumYrs=length(unique(year)),
            From=min(year),
            To=max(year)) %>% 
  knitr::kable()

load(here("data","bgcClim.rda"))

clim.X<-
  bgcClim %>% 
  pivot_longer(cols=Tmax01:last_col()) %>% 
  mutate(Period=paste(period,Scenario,sep="_")) %>% 
  dplyr::select(-period,-Scenario) %>% 
  pivot_wider(names_from="Period",values_from="value")



future.ppt<-
  clim.X %>% 
  dplyr::filter(str_detect(name,"PPT")) %>% 
  mutate_at(vars(matches("_rcp")), list(dif = (~ . / `1961 - 1990_NA`)))

future.tmx<-
  clim.X %>% 
  dplyr::filter(str_detect(name,"Tmax")) %>% 
  mutate_at(vars(matches("_rcp")), list(dif = (~ . - `1961 - 1990_NA`)))


future.tmn<-
  clim.X %>% 
  dplyr::filter(str_detect(name,"Tmin")) %>% 
  mutate_at(vars(matches("_rcp")), list(dif = (~ . - `1961 - 1990_NA`)))


# Future climate data
  future.Data<-
    rbind(future.ppt,future.tmx,future.tmn) %>% 
    mutate(month=str_sub(name,start=-2,end=-1)) %>% 
    mutate(climVar=str_sub(name,start=1,end=-3)) %>% 
    dplyr::select(BGC,climVar,month,contains("dif")) %>% 
    pivot_longer(`2025_rcp45_dif`:last_col(),names_to="scenario",values_to="values") %>% 
    mutate(climVar=paste(climVar,"future",sep="_")) %>% 
    pivot_wider(names_from="climVar",values_from="values") %>% 
    mutate(month=as.integer(month))

save(future.Data,file=here::here("data-raw","WxData_newBECvariants","futureData.RData"))

load(here::here("data-raw","WxData_newBECvariants","futureData.RData"))

BGC unit ASMR

asmrRun1<-function(stnID,bgc) {

  # set seed to control random year 
  set.seed(5)  

    # Generate random years
  asmrData<-
    climData %>% 
    pluck("cleanedData") %>% 
    .[.$stn%in%stnID,] %>% # for some reason dplyr::filter isn't working
    filter(stn %in% stnID) %>% 
    rename(ppt="ppt_filled",
           tmx="tmx_filled",
           tmn="tmn_filled") %>% 
    mutate(date=paste(year,month,day,sep="-")) %>% 
    asmrCalc()


  yearData<-
    asmrData %>% 
    yearSelect(excl=c(1970,2020),win=c(1980,2005),yrs=10)



  asmrSummary<- 
   asmrData %>% 
    filter(year %in% yearData) %>% # filter for years
    group_by(month) %>% 
    summarise_at(vars(ends_with(".ASMR")),mean) %>%
        ungroup() %>%

        rbind(c(13,colMeans(.)[-1])) %>%
        mutate(month=cut(month,breaks=13,labels=c(month.abb,"Annual"))) %>% 
    mutate(Station=str_to_title(stnID),BGC=bgc) %>% 
    mutate_at(vars(ends_with(".ASMR")),round,3) %>% 
    dplyr::select(Station,BGC,everything()) 



    # generate future AET/PET ratios
 future.X<-
    climData %>% 
    pluck("cleanedData") %>% 
    .[.$stn%in%stnID,] %>% # for some reason dplyr::filter isn't working
    filter(stn %in% stnID) %>% 
    rename(ppt="ppt_filled",
           tmx="tmx_filled",
           tmn="tmn_filled") %>% 
    right_join(filter(future.Data,BGC==bgc),by="month") %>% 
    mutate(ppt=PPT_future*ppt,
           tmx=Tmax_future+tmx,
           tmn=Tmin_future+tmn) %>% 
   mutate(date=paste(year,month,day,sep="-")) %>% 
   dplyr::select(stn,BGC,year,month,day,date,ppt,tmx,tmn,scenario)

   futureASMR<-rbind(data.frame(Scenario="rcp45_2025",asmrCalc(filter(future.X,scenario=="2025_rcp45_dif"))),
    data.frame(Scenario="rcp45_2055",asmrCalc(filter(future.X,scenario=="2055_rcp45_dif"))),
    data.frame(Scenario="rcp85_2025",asmrCalc(filter(future.X,scenario=="2025_rcp85_dif"))),
    data.frame(Scenario="rcp85_2055",asmrCalc(filter(future.X,scenario=="2055_rcp85_dif")))) %>% 

   # create summaries
       filter(year %in% yearData) %>% # filter for years
    group_by(Scenario,month) %>% 
    summarise_at(vars(ends_with(".ASMR")),mean) %>%
        ungroup() %>%
        mutate(Station=str_to_title(stnID),BGC=bgc) %>% 
    mutate_at(vars(ends_with(".ASMR")),round,3) %>% 
     mutate(month=as.integer(month)) %>% 
    dplyr::select(Station,BGC,everything()) 


  futureSummary<- 
     futureASMR %>% 
       mutate(month=as.integer(month)) %>% 
      ungroup() %>% 
     group_by(Scenario) %>% 
     summarise_at(vars(ends_with(".ASMR")),mean) %>% 
        mutate(Station=str_to_title(stnID),BGC=bgc,month=13) %>%
           dplyr::select(Station,BGC,month,everything()) %>% 
       rbind(futureASMR) %>% 
          arrange(Scenario,month) %>% 
         mutate(month=cut(month,breaks=13,labels=c(month.abb,"Annual"))) %>% 
    mutate_at(vars(ends_with(".ASMR")),round,3) %>% 
   mutate(Period=str_split_fixed(Scenario,patter="_",n=2)[,2]) %>% 
     mutate(Scenario=str_split_fixed(Scenario,patter="_",n=2)[,1]) %>% 
    dplyr::select(-Station) %>% 
    dplyr::select(BGC,Period,Scenario,Month=month,everything())



  return(list(years=yearData,asmr=asmrSummary,future=futureSummary))




}

For each BGC unit, I chose the station with the best climate data coverage. we can run this again if there's a better station to use.

IDFdm1

stnID="marysville"
bgc="IDFdm1"

x<-asmrRun1(stnID,bgc)

For this BGC unit, I used r stnID station, and randomly selected the following years: r x$years. Here is the monthly and annual AET/PET ratio for those 10 years:

x$asmr %>% 
  knitr::kable()