R/evaluate_TLstatistics.R

In pik-piam/validation: MAgPIE validation library

evaluate_TLstatistics<-function(TLstatistics,aggregate=TRUE){
  ###################################################
  #Construct output objects
  ###################################################
  #construct output object for the aggregated traffic light
  TrafficLight<-new.magpie(getCells(TLstatistics),NULL,fulldim(TLstatistics)[[2]][[5]])
  #construct intermediate traffic light output for individual tests
  TLtmp<-new.magpie(getCells(TLstatistics),NULL,fulldim(TLstatistics)[[2]][[4]])
  
  #####################################################
  #Define thresholds. 3 thresholds for each parameter 
  #   for the two classes green and yellow. 
  #   If none of the thresholds is met, the TL (TrafficLight) will be red.
  #1. Upper threshold (up)
  #2. lower threshold (low)
  #3. define a value that if it is included in the 
  #   range of statistic tests among different datasets,
  #   sets the TL to yellow or green (includes)
  ########################################################
  
  #Define object
  names<-as.vector(outer(fulldim(TLstatistics)[[2]][[4]],c("up","low","includes"),paste,sep="."))
  names<-as.vector((outer(names,c("yellow","green"),paste,sep=".")))
  thresholds<-new.magpie("GLO",NULL,names)
  
  #Thresholds for GRI:
  thresholds[,,"L_GRI"][,,"yellow"][,,"up"]<-1.2
  thresholds[,,"L_GRI"][,,"yellow"][,,"low"]<-1
  thresholds[,,"L_GRI"][,,"yellow"][,,"includes"]<-NA
  
  thresholds[,,"L_GRI"][,,"green"][,,"up"]<-1.1
  thresholds[,,"L_GRI"][,,"green"][,,"low"]<-1
  thresholds[,,"L_GRI"][,,"green"][,,"includes"]<-NA
  
  #Thresholds for MRSR:
  thresholds[,,"L_MRSR"][,,"yellow"][,,"up"]<-3
  thresholds[,,"L_MRSR"][,,"yellow"][,,"low"]<- 0
  thresholds[,,"L_MRSR"][,,"yellow"][,,"includes"]<-NA
  
  thresholds[,,"L_MRSR"][,,"green"][,,"up"]<-1
  thresholds[,,"L_MRSR"][,,"green"][,,"low"]<- 0
  thresholds[,,"L_MRSR"][,,"green"][,,"includes"]<-NA
  
  #Thresholds for RD:
  thresholds[,,"L_RD"][,,"yellow"][,,"up"]<-0.3
  thresholds[,,"L_RD"][,,"yellow"][,,"low"]<- -0.3
  thresholds[,,"L_RD"][,,"yellow"][,,"includes"]<-0
  
  thresholds[,,"L_RD"][,,"green"][,,"up"]<- 0.1
  thresholds[,,"L_RD"][,,"green"][,,"low"]<- -0.1
  thresholds[,,"L_RD"][,,"green"][,,"includes"]<-0
  
  #Thresholds for 40 year trend p value of Mann Kendall test:
  thresholds[,,"T_MK_p"][,,"yellow"][,,"up"]<-1
  thresholds[,,"T_MK_p"][,,"yellow"][,,"low"]<- 0.01
  thresholds[,,"T_MK_p"][,,"yellow"][,,"includes"]<-NA
  
  thresholds[,,"T_MK_p"][,,"green"][,,"up"]<- 1
  thresholds[,,"T_MK_p"][,,"green"][,,"low"]<- 0.1
  thresholds[,,"T_MK_p"][,,"green"][,,"includes"]<-NA
  
  #Thresholds for 40 year trend m_annual
  thresholds[,,"T_MK_m_annual"][,,"yellow"][,,"up"]<-0.01
  thresholds[,,"T_MK_m_annual"][,,"yellow"][,,"low"]<- -0.01
  thresholds[,,"T_MK_m_annual"][,,"yellow"][,,"includes"]<-0
  
  thresholds[,,"T_MK_m_annual"][,,"green"][,,"up"]<- 0.005
  thresholds[,,"T_MK_m_annual"][,,"green"][,,"low"]<- -0.005
  thresholds[,,"T_MK_m_annual"][,,"green"][,,"includes"]<-0
  
  #Thresholds for overlap p value of Mann Kendall test:
  thresholds[,,"O_MK_p"][,,"yellow"][,,"up"]<-1
  thresholds[,,"O_MK_p"][,,"yellow"][,,"low"]<- 0.01
  thresholds[,,"O_MK_p"][,,"yellow"][,,"includes"]<-NA
  
  thresholds[,,"O_MK_p"][,,"green"][,,"up"]<- 1
  thresholds[,,"O_MK_p"][,,"green"][,,"low"]<- 0.1
  thresholds[,,"O_MK_p"][,,"green"][,,"includes"]<-NA
  
  #Thresholds for overlap m_annual
  thresholds[,,"O_MK_m_annual"][,,"yellow"][,,"up"]<-0.01
  thresholds[,,"O_MK_m_annual"][,,"yellow"][,,"low"]<- -0.01
  thresholds[,,"O_MK_m_annual"][,,"yellow"][,,"includes"]<-0
  
  thresholds[,,"O_MK_m_annual"][,,"green"][,,"up"]<- 0.005
  thresholds[,,"O_MK_m_annual"][,,"green"][,,"low"]<- -0.005
  thresholds[,,"O_MK_m_annual"][,,"green"][,,"includes"]<-0  
  
  
  
  ####################################################
  #calculate traffic light for individual tests 
  #0 = red
  #1 = yellow
  #2 = green
  #####################################################
  TLtmp[,,]<-0
  for(stat in getNames(TLtmp)){
    #test for the two classes yellow and green
    for(class in c("yellow","green")){
      #test range(x TLstatistics between low and up thresholds
      if(all(!is.na(thresholds[,,stat][,,class][,,c("low","up")])) & !all(is.na(TLstatistics[,,stat]))){
        low_test <- as.numeric(thresholds[,,stat][,,class][,,"low"])<=TLstatistics[,,stat]
        up_test  <- as.numeric(thresholds[,,stat][,,class][,,"up"])>=TLstatistics[,,stat]
        range_test<-low_test&up_test
      } else {
        range_test<-TLstatistics[,,stat]
        range_test[,,]<-NA
      }
      #test if TLstatistics for different datasets include "includes" value
      if(!is.na(thresholds[,,stat][,,class][,,"includes"]) & !all(is.na(TLstatistics[,,stat]))){
        includes_test_up<-  as.numeric(thresholds[,,stat][,,class][,,"includes"])<=max(TLstatistics[,,stat],na.rm=TRUE)
        includes_test_low<-  as.numeric(thresholds[,,stat][,,class][,,"includes"])>=min(TLstatistics[,,stat],na.rm=TRUE)
        includes_test<-includes_test_up & includes_test_low
      } else {
        includes_test<-NA
      }
      #if there is one TRUE, the test has been passed
      passed<-suppressWarnings(any(c(range_test,includes_test)))
      #Dirty piece to set to FALSE if NAs and no TRUES are present
      if(is.na(passed)){
        if(sum(any(!c(range_test,includes_test)),na.rm=TRUE)>0) passed<-FALSE
      }
      if(is.na(passed)){
        TLtmp[,,stat]<-NA
      } else {
        if(passed) TLtmp[,,stat]<-TLtmp[,,stat]+1 
      }
    }      
  }
  
  if(aggregate){
    ##############################3
    #Aggregate the p value and m_annual test
    #for the trend by taking the better one of both
    ##############################
    TLtmp[,,"T_MK_p"]<-suppressWarnings(max(TLtmp[,,c("T_MK_p","T_MK_m_annual")],na.rm=TRUE))
    if(is.infinite(TLtmp[,,"T_MK_p"])) TLtmp[,,"T_MK_p"]<-NA
    TLtmp[,,"T_MK_m_annual"]<-NA

    TLtmp[,,"O_MK_p"]<-suppressWarnings(max(TLtmp[,,c("O_MK_p","O_MK_m_annual")],na.rm=TRUE))
    if(is.infinite(TLtmp[,,"O_MK_p"])) TLtmp[,,"O_MK_p"]<-NA
    TLtmp[,,"O_MK_m_annual"]<-NA
    
    
    ####################################
    #Take the mean of all tests 
    #and round to get the final traffic light
    ####################################
    TrafficLight[,,]<-mean(TLtmp,na.rm=TRUE)
    TrafficLight<-round(TrafficLight)
    if(is.infinite(TrafficLight)) TrafficLight[,,]<-NA
  } else {
    TrafficLight<-TLtmp
  }
  
  return(TrafficLight)
}