In RETURN-project/GEE.aux: What the Package Does (One Line, Title Case)

knitr::opts_chunk$set(echo = TRUE, eval = F)

library(zoo)
library(GEE.aux)
library(tidyverse)
library(bfast)
library(strucchange)
library(sf)
library(raster)
library(ggplot2)
library(lubridate)

User inputs

# path to folder with data
ifolder <- '/home/wanda/Documents/data/upscaleRecovery/test/sample_n10000_GEEfire'
# path to folder where outputs should be stored
ofolder <- '/home/wanda/Documents/data/upscaleRecovery/test/sample_n10000_GEEfire'
# name of data files
ifile <- 'samples_forest_fire_nat_n10000_' # base name of the input file
startPix <- 0 # first time series of the first input file to be processed
endPix <- 348000# first time series of the last input file to be processed
intPix <- 2000# number of time series per input file

confLevFire <- 95# mimimum confidence level to recognize a fire event

obspyr <- 12 # number of observations per year
nPre <- 2 #number of years used to quantify the pre-disturbance state
nDist <- 1 # number of years used to quantify the disturbance state
nPost <- 2 # number of years used to quantify the post-disturbance state
nPostStart <- 4# number of years between the disturbance and start of the post-disturbance period
nDelta <- 2# number of years used to quantify the post-disturbance state for the YrYr metric
nDeltaStart <- 4# number of years between the disturbance and start of the post-disturbance period for the YrYr metric
selBreak <- 'first'# which break should be used to calculate the recovery metrics?
chkBrk <- T# check if no additional breaks occur in the period that recovery is assessed
timeThres <- 1.5# the maximum allowed time difference between the reported disturbance and detected break

Read data into tidy data frame

for (i in seq(startPix,endPix,by =intPix)){
  if(i == endPix){
    endpnt = 'None'
    endStr <- 'None'
  } else{
    endpnt = i+intPix
    endStr <- sprintf('%i',endpnt)}#'None'
    # read the input file
  x <- readLines(file.path(ifolder,paste0(ifile, sprintf('%i',i), '_', endStr, '.csv')))
  # column names of the data 
  col_names <- c("id","ecoReg","coords","img_names","series","img_doy_","doy_fire_","img_conf_","conf_fire_")
  # generate a dataframe from the input data 
  dat_tidy <- tidy_data(x,col_names)
  # write the dataframe to a csv file
  write_csv(dat_tidy, file.path(ofolder, paste0('tidy_',ifile, sprintf('%i',i), '_', endStr, '.csv')))
  save(dat_tidy, file = file.path(ofolder, paste0('tidy_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))

  # column names of the dataframe
  cols <- names(dat_tidy)# all column names
  dtcols <- cols[startsWith(cols,"p__")]# column names of the columns containing the observation dates
  obscols <- cols[startsWith(cols,"series__")]# column names of the columns containing the observations
  auxcols <- cols[!startsWith(cols,"p__") & !startsWith(cols,"series__") & !startsWith(cols,"img_doy_") & !startsWith(cols,"doy_fire_") & !startsWith(cols,"img_conf_") & !startsWith(cols,"conf_fire_")]# column names of the columns containing the auxiliary information
  firecols <- cols[!startsWith(cols,"p__") & !startsWith(cols,"series__")]# column names of the columns containing fire-related info 
  col_names <- list(cols = cols, dtcols = dtcols, obscols = obscols, auxcols = auxcols, firecols = firecols)
  save(col_names, file = file.path(ofolder, paste0('col_names_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))
  rm(x, col_names, dat_tidy, cols, dtcols, obscols, auxcols, firecols)
}

Extract fire dates

for (i in seq(startPix,endPix,by =intPix)){
  if(i == endPix){
    endpnt = 'None'
    endStr <- 'None'
    }else{
      endpnt = i+intPix
      endStr <- sprintf('%i',endpnt)}
  # load data
  load(file.path(ofolder, paste0('tidy_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))
  load(file.path(ofolder, paste0('col_names_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))
  firecols <- col_names$firecols
  # get the columns with fire related information
  dat_fire <- dat_tidy[,firecols]
  # extract time series of fire doy and fire confidence
  tidy_fire <- tidy_fire_date(dat_fire)
  fire_doy <- tidy_fire$fire_doy
  fire_conf <- tidy_fire$fire_conf
  firedts <- as.Date(names(fire_conf)[startsWith(names(fire_conf),"1") | startsWith(names(fire_conf),"2")])

  # get fire events per time series
  fireDts <- lapply(1:dim(fire_doy)[1], function(i){getFireDate(as.numeric(fire_conf[i,as.character(firedts)]), as.numeric(fire_doy[i,as.character(firedts)]),firedts, confLevFire)})
  # save results
  save(fireDts, file = file.path(ofolder, paste0('fire_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))# save the result
  rm(firecols, dat_tidy, col_names, tidy_fire, fire_doy, fire_conf, firedts, fireDts)
}

convert the irregular time series to regular ones

library(parallel)
numCores <- parallel::detectCores()-1 # Requires library(parallel)
print(numCores)

mclapply(seq(startPix,endPix,by =intPix), function(i){
# for (i in seq(startPix,endPix,by =intPix)){
  if(i == endPix){
    endpnt = 'None'
    endStr <- 'None'
  }else{
      endpnt = i+intPix
      endStr <-  sprintf('%i',endpnt)}
  load(file.path(ofolder, paste0('tidy_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))
  load(file.path(ofolder, paste0('col_names_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))
  dtcols <- col_names$dtcols
  obscols <- col_names$obscols
  auxcols <- col_names$auxcols
  # set desired start and end dates of the regular time series
  dat_tidy <- dat_tidy%>%mutate_at(dtcols,function(x){as.numeric(x)})# make dates numeric
  mnyr <- format(as.Date(as.character(min(dat_tidy[,dtcols], na.rm = T)),'%Y%m%d'),'%Y')# start year of the irregular time series
  mxyr <- format(as.Date(as.character(max(dat_tidy[,dtcols], na.rm = T)),'%Y%m%d'),'%Y')# end year of the irregular time series
  startdt <- as.Date(paste0(mnyr,'-01-01'))# desired start date of regular time series
  enddt <- as.Date(paste0(mxyr,'-12-31'))# desired end date of regular time series
  dts <- seq(startdt,enddt,by='month')# dates of regular series

  # extract dates and associated observations of the irregular time series
  # add a NA observation at the desired start and end date of the regular series
  # the latter is required to ensure that the time span of the regular series is fixed between these two dates
  dtsi <- dat_tidy[,dtcols]# dates of the irregular time series
  dtsi <- cbind(rep(as.numeric(paste0(mnyr,'0101')),dim(dtsi)[1]),rep(as.numeric(paste0(mxyr,'1231'),dim(dtsi)[1])),dtsi)# dates, extended with desired start and end date
  dat_tidy <- dat_tidy%>%mutate_at(obscols,function(x){as.numeric(x)})# make dates numeric
  obsi <- dat_tidy[,obscols]# observations of the irregular time series
  obsi <- cbind(rep(NA,dim(obsi)[1]),rep(NA,dim(obsi)[1]),obsi)# observations, extended with two NA values

  # Generate a regular time series for each pixel
  out <- sapply(1:dim(obsi)[1], function(i) getRegularTS(dtsi[i,], obsi[i,]))# generate a regular time series
  df_out <- as.data.frame(t(out))# convert to a dataframe
  names(df_out) <- dts# add dates of observations to column names
  df_out <- cbind(dat_tidy[,auxcols],df_out)# add auxiliary data to dataframe

  save(df_out, file = file.path(ofolder, paste0('reg_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))# save the result
  rm(col_names, dat_tidy, dtcols, obscols, auxcols, mnyr, mxyr, startdt, enddt, dts, dtsi, obsi, out, df_out)
}, mc.cores = numCores)

Piecewise regression

library(parallel)
numCores <- parallel::detectCores()-1 # Requires library(parallel)
print(numCores)

mclapply(seq(startPix,endPix,by =intPix), function(i){
  if(i == endPix){endpnt = 'None'}else{endpnt = i+intPix}
  load(file.path(ofolder, paste0('reg_', ifile, sprintf('%i',i), '_', sprintf('%i',endpnt), '.Rdata')))
  # number of pixels to be processed
  nts <- dim(df_out)[1]
  # fit piecewise regression for each pixel
  seg <- lapply(1:nts, function(i){getSegments(as.numeric(df_out[i,-c(1:4)]),12,0.15,T)})
  # save result
  save(seg, file = file.path(ofolder, paste0('seg_', ifile, sprintf('%i',i), '_', sprintf('%i',endpnt), '.Rdata')))
  rm(nts, df_out, seg)
}, mc.cores = numCores)

# for (i in seq(startPix,endPix,by =intPix)){
#   if(i == endPix){endpnt = 'None'}else{endpnt = i+intPix}
#   load(file.path(ofolder, paste0('reg_', ifile, sprintf('%i',i), '_', sprintf('%i',endpnt), '.Rdata')))
#   # number of pixels to be processed
#   nts <- dim(df_out)[1]
#   # fit piecewise regression for each pixel
#   seg <- lapply(1:nts, function(i){getSegments(as.numeric(df_out[i,-c(1:4)]),12,0.15,T)})
#   # seg <- mclapply(1:nts, function(i){getSegments(as.numeric(df_out[i,-c(1:4)]),12,0.15,T)}, mc.cores = numCores)
#   # save result
#   save(seg, file = file.path(ofolder, paste0('seg_', ifile, sprintf('%i',i), '_', sprintf('%i',endpnt), '.Rdata')))
#   rm(nts, df_out, seg)
# }

Calculate recovery metrics

for (i in seq(startPix,endPix,by =intPix)){
  if(i == endPix){
    endpnt = 'None'
    endStr <- 'None'
  }else{
      endpnt = i+intPix
      endStr <-  sprintf('%i',endpnt)}
  # load(file.path(ofolder,paste0('rec_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))
  load(file.path(ofolder,paste0('seg_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))
  load(file.path(ofolder,paste0('reg_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))
  load(file.path(ofolder,paste0('fire_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))

  nts <- dim(df_out)[1]
  dts <- as.Date(names(df_out)[-c(1:4)])

  # calculate recovery metrics
  df_lst <-  lapply(1:nts,function(ii){
    # get observation number of reported disturbance
    fdts <- fireDts[[ii]][1]# only take the first fire
    fdts[!is.na(fdts)] <- rollback(fdts[!is.na(fdts)], roll_to_first = T)
    tdist <- which(dts %in% fdts)# convert fire date to observation number
    # trend component of piecewise regression
    tsi <- seg[[ii]]$trend
    # breakpoints of piecewise regression
    tbp <- seg[[ii]]$breakpoints
    tbp <- tbp[tbp!=0 & tbp!=length(tsi)]
    # calculate recovery metric for each fire event
    rec <- sapply(1:length(fdts), function(iii){
      if(!is.na(tdist[iii])){
        calcRecMetrics(tsi, tdist[iii], obspyr, nPre, nDist, nPost, nPostStart, nDelta, nDeltaStart, tbp, selBreak, chkBrk, timeThres)
      }else{
        list('RRI' = NA, 'R80P' = NA, 'YrYr' = NA , 'impact' = NA, 'Vpre' = NA, 'Break' = NA,'timeChck' = NA,'postChck' = NA, 'distChck' = NA, 'brkChck' = NA, 'brkChckPre' = NA)
      }})
  })
  df_fr <- as.data.frame(t(data.frame(do.call(cbind.data.frame,df_lst))))
  df_fr$RRI <- unlist(df_fr$RRI)
  df_fr$R80P <- unlist(df_fr$R80P)
  df_fr$YrYr <- unlist(df_fr$YrYr)
  df_fr$impact <- unlist(df_fr$impact)
  df_fr$Vpre <- unlist(df_fr$Vpre)
  df_fr$Break <- unlist(df_fr$Break)
  df_fr$timeChck <- unlist(df_fr$timeChck)
  df_fr$postChck <- unlist(df_fr$postChck)
  df_fr$distChck <- unlist(df_fr$distChck)
  df_fr$brkChck <- unlist(df_fr$brkChck)
  df_fr$brkChckPre <- unlist(df_fr$brkChckPre
                             )

  # add auxiliary data to the dataframe
  df_rec <- cbind(df_out[,c(1:4)], df_fr)#dat_tidy[,auxcols]
  # save result
  save(df_rec, file = file.path(ofolder,  paste0('rec_', ifile, sprintf('%i',i), '_', endStr, '.Rdata')))
  rm(df_out, nts, dts, df_lst, fireDts, seg, df_fr, df_rec)#fdts,tdist, tsi,tbp,  rec,
}

Merge data

for (i in seq(startPix,endPix,by =intPix)){
  if(i == endPix){
    endpnt = 'None'
    endStr <- 'None'
  }else{
      endpnt = i+intPix
      endStr <-  sprintf('%i',endpnt)}
  # if(i == endPix){endpnt = 'None'}else{endpnt = i+intPix}

  load(file.path(ofolder,paste0('rec_', ifile, sprintf('%i',i),'_',endStr,'.Rdata')))
  if(i==0){df_rec_tot <- df_rec}else{df_rec_tot <- rbind(df_rec_tot, df_rec)}
  rm(df_rec)

  load(file.path(ofolder,paste0('seg_', ifile,sprintf('%i',i),'_',endStr,'.Rdata')))
  if(i==0){seg_tot <- seg}else{seg_tot <- c(seg_tot, seg)}
  rm(seg)

  load(file.path(ofolder,paste0('reg_', ifile,sprintf('%i',i),'_',endStr,'.Rdata')))
  if(i==0){df_out_tot <- df_out}else{df_out_tot <- rbind(df_out_tot, df_out)}
  rm(df_out)

  load(file.path(ofolder,paste0('fire_', ifile,sprintf('%i',i),'_',endStr,'.Rdata')))
  if(i==0){fireDts_tot <- fireDts}else{fireDts_tot <- c(fireDts_tot, fireDts)}
  rm(fireDts)
}

save(df_rec_tot, file = file.path(ofolder,paste0('df_rec_tot')))
save(df_out_tot, file = file.path(ofolder,paste0('df_out_tot')))
save(fireDts_tot, file = file.path(ofolder,paste0('fireDts_tot')))
save(seg_tot, file = file.path(ofolder,paste0('seg_tot')))

Classes of pixels that were not processed

df_rec_tot <- df_rec_tot[complete.cases(df_rec_tot),]
# fraction of pixels with break
sum(df_rec_tot$Break,na.rm = T)/ length(df_rec_tot$Break)
# fraction of pixels with correct time span between break and fire
sum(df_rec_tot$timeChck == F,na.rm = T)/ length(df_rec_tot$Break)
# positive post-disturbance slope
sum(df_rec_tot$postChck == F,na.rm = T)/ length(df_rec_tot$Break)
# negative break
sum(df_rec_tot$distChck == F,na.rm = T)/ length(df_rec_tot$Break)
# no negative break in recovery period
sum(df_rec_tot$brkChck == F,na.rm = T)/ length(df_rec_tot$Break)
# no negative break in pre-disturbance period
sum(df_rec_tot$brkChckPre == F,na.rm = T)/ length(df_rec_tot$Break)
# total processed
sum(is.na(df_rec_tot$R80P)==F)/ length(df_rec_tot$Break)
# pixels 
sum((df_rec_tot$Break & df_rec_tot$timeChck & df_rec_tot$postChck & df_rec_tot$distChck & df_rec_tot$brkChck & df_rec_tot$brkChckPre),na.rm = T)/ length(df_rec_tot$Break)

totChck <- (df_rec_tot$Break & df_rec_tot$timeChck & df_rec_tot$postChck & df_rec_tot$distChck & df_rec_tot$brkChck & df_rec_tot$brkChckPre)

which (is.na(df_rec_tot$R80P) & totChck == T)

# not processed since not enough data before or after break
sum(is.na(df_rec_tot$R80P) & totChck == T, na.rm = T)/ length(df_rec_tot$Break)

which(is.na(df_rec_tot$Break))

Generate shapefile of all valid recovery indicators

df_rec_val <- df_rec_tot[complete.cases(df_rec_tot),]
df_rec_val <- df_rec_val %>%
      separate(coords,into=c('lon','lat'),sep=", ",convert=TRUE)
coordinates(df_rec_val)=~lon+lat
proj4string(df_rec_val)<- CRS("+proj=longlat +datum=WGS84")
raster::shapefile(df_rec_val, file.path(ofolder, paste0(ifile, "VAL.shp")))
# this shapefile should be uploaded to the GEE in order to dll covariates