In RETURN-project/BenchmarkRecovery: Benchmark recovery indicators derived from remote sensing time series

knitr::opts_chunk$set(echo = TRUE, 
                      eval = FALSE,
                      fig.width=18,
                      fig.height=10)
library(BenchmarkRecovery)
library(reshape2)
library(plyr)

# inputs
# input folder
ifolder <- '../data/'#'/home/wanda/Documents/data/benchmarkRecovery/Run_20201016_expRecPeriod/'#ifolder <- '/home/wanda/Documents/data/benchmarkRecovery/Run_20201005_4steps/'#'../data/'#ifolder <- '/home/wanda/Documents/data/benchmarkRecovery/Run_20201005_4steps/'
# Folder where outputs will be written
ofolder <- '../data/'#'/home/wanda/Documents/data/benchmarkRecovery/Run_20201016_expRecPeriod/Figures/'#ofolder <-'/home/wanda/Documents/data/benchmarkRecovery/Run_20201005_4steps/Figures'#'../data/' #ofolder <- '/home/wanda/Documents/data/benchmarkRecovery/Run_20201005_4steps/Figures'
# Name of the input dataset
basename <- 'LSTS_RndmSample_NoFire_5_Tree_80_scl_30_npnt_20000_VI' 
caseList <- c('seasAmp','remSd','distT','distRec','missVal','distMag')# 'distMag',evaluated time series characteristics for 

simFullName  <- list('Disturbance magnitude',
                     'Disturbance timing',
                     'Recovery period',
                     'Seasonal amplitude',
                     'Noise level',
                     'Missing values')
names(simFullName) <- c('distMag',
                        'distT',
                        'distRec',
                        'seasAmp',
                        'remSd',
                        'missVal')
metric_list <- c('MAPE', 'R2', 'nTS', 'RMSE')
metric_names <- c('MAPE', 'R²', 'nTS', 'RMSE'); names(metric_names) <- c('MAPE', 'R2', 'nTS', 'RMSE')
tempRes_list <- c('quarterly', 'annual', 'dense', 'all')

dat_list <- list()
for(mm in 1:length(metric_list)){
  metric <- metric_list[mm]
  for(vr in 1:length(caseList)){
  evr <- caseList[vr]# name of parameter 

  RRI_dat <- loadRData(file.path(ifolder, paste0(basename, '_RRI_', metric, '_' , evr, '.rda')))
  R80p_dat <- loadRData(file.path(ifolder, paste0(basename, '_R80p_', metric, '_' , evr, '.rda')))
  YrYr_dat <- loadRData(file.path(ifolder, paste0(basename, '_YrYr_', metric, '_' , evr, '.rda')))

  tot_dat <- melt(rbind(RRI_dat, R80p_dat, YrYr_dat))

  tot_dat$Period <- revalue(factor(tot_dat$nPostMin), c("1"="Short", "4"="Long"))

   if((evr == 'remSd') || (evr == 'seasAmp') || (evr == 'missVal')) {
    tot_dat$variable <-mapvalues(tot_dat$variable, from = levels(tot_dat$variable), to = c("very low", "low", 'medium', 'high'))
    }  else{
    tot_dat$variable <-mapvalues(tot_dat$variable, levels(tot_dat$variable), to = c("very low", "low", 'medium', 'high'))
    }

  tot_dat$param = simFullName[[caseList[[vr]]]]
  if(vr == 1){totp_dat <- tot_dat}else{totp_dat <- rbind(totp_dat,tot_dat)}
  }
  totp_dat$param[totp_dat$param == 'SD remainder'] <- 'Noise level'
  totp_dat$paramType <- 'Environmental parameter'
  totp_dat[(totp_dat$param == 'Disturbance magnitude' | totp_dat$param == 'Recovery period' | totp_dat$param == 'Disturbance timing'), ]$paramType <- 'Disturbance parameter'
  totp_dat$param <- factor(totp_dat$param, levels = simFullName)

  dat_list[[metric]] <- totp_dat
}

Compare the performance of each recovery indicator

for(mm in 1:length(metric_list)){
  metric <- metric_list[mm]
  # plot 
  plt <- plotMet(dat_list[[metric]],  'Metric', metric_names[[metric]])
  png(file.path(ofolder, paste0(basename, '_', metric, 'Met.png')),width = 1311,height =628 )
  print(plt)
  dev.off()
  print(plt)
}

Which characteristics influence the performance the most?

for(tr in 1:length(tempRes_list)){
  tempRes <- tempRes_list[tr]
  for(mm in 1:length(metric_list)){
    metric <- metric_list[mm]
    data <- dat_list[[metric]]
    if(tempRes != 'all'){data <- data[data$Dense == tempRes, ]}
    xlbl <- 'Parameter value'
    ylbl <- metric_names[[metric]]
    scales = 'free_y'
    # plot 
    plt <- plotSensBar(data, xlbl, ylbl, scales)
    png(file.path(ofolder, paste0(basename, '_', tempRes, '_',metric,'_Env.png')),width = 1911,height =1828 )
    print(plt)
    dev.off()
  }

}

for(mm in 1: length(metric_list)){
  metric <- metric_list[mm]
  for(pp in 1:length(simFullName)){
    data <- dat_list[[metric]]
    data <- data[data$param == simFullName[pp],]
    data$Dense <- revalue(factor(data$Dense, levels = c('dense', 'quarterly', 'annual')), c("dense" = "no", "annual"="annual", "quarterly" = "quarterly"))
  data$Smooth <- revalue(factor(data$Smooth, levels = c('raw', 'smoothed', 'segmented')), c("raw"="no", "smoothed"="rolling mean", "segmented" = "segmentation"))
  xlbl <- simFullName[pp]
  ylbl <- metric_names[[metric]]
  scales = 'free_y'
  plt <- plotSens(data, xlbl, ylbl, scales)
  png(file.path(ofolder, paste0(basename, '_',metric,'_', names(simFullName[pp]),'_Sens.png')),width = 1911,height =1828 )
    print(plt)
    dev.off()
  }
}

How can we improve the performance?

for(mm in 1:length(metric_list)){
  metric <- metric_list[mm]
  data <- dat_list[[metric]]
  data$Dense <- revalue(factor(data$Dense, levels = c('dense', 'quarterly', 'annual')), c("dense" = "no", "annual"="annual", "quarterly" = "quarterly"))
  data$Smooth <- revalue(factor(data$Smooth, levels = c('raw', 'smoothed', 'segmented')), c("raw"="no", "smoothed"="rolling mean", "segmented" = "segmentation"))

  xlbl <- 'Temporal aggregation'
  ylbl <- metric_names[[metric]]
  scales = 'free_y'
  # plot 
  plt <-pltPrepBox(data, xlbl, ylbl, scales) 
  png(file.path(ofolder, paste0(basename, '_', metric,'_Prep.png')),width = 1911,height =1828 )
  print(plt)
  dev.off()
}

RETURN-project/BenchmarkRecovery documentation built on July 13, 2021, 5:47 p.m.