R/histoPlot.R

In pik-piam/validation: MAgPIE validation library

histoPlot<-function(magpie,historical=NULL,projection=NULL,years=NULL,xlim=NULL,ylim=NULL,index=FALSE,same_yscale=TRUE,text_size=15,title="",ylab="",legend_names=c("Model output","Historical Data","Other projections"),raw_plot=FALSE,raw_legend=FALSE){
  require(ggplot2, quietly = TRUE)
  require(RColorBrewer, quietly = TRUE)
  require(plyr, quietly = TRUE)
  require(gridExtra, quietly = TRUE)

  if(raw_legend&!raw_plot)stop("Conflicting arguments: raw_legend=TRUE only works with raw_plot=TRUE")
  ####################################################
  #Get the magpie data
  ####################################################
  if(!is.list(magpie)){
    tmp<-list()
    tmp[["MAgPIE"]]<-magpie
    magpie<-tmp
    rm(tmp)
  }
  if(length(magpie)==1 & is.null(names(magpie))) names(magpie)<-"MAgPIE"
  #Get the required years if the years argument is not NULL
  reduce_years<-function(x,years){
    x_years<-getYears(x,as.integer=TRUE)
    if(is.na(years[1])) years[1]<-min(x_years)
    if(is.na(years[2])) years[2]<-max(x_years)
    x<-x[,which(x_years>=years[1]&x_years<=years[2]),]
    return(x)
  }
  if(!is.null(years)) magpie<-lapply(X=magpie,FUN=reduce_years,years)
  ####################################################
  #Determine the aggregation level of the magpie data
  #Could be replaced by .aggrLevel 
  ####################################################
  level<-.aggrLevel(magpie[[1]])
  if(level=="unknown") stop("Unknown aggregation level specified.")
  if(level=="country") stop("country aggregation level not supported.")
  if(level=="cell") stop("Cellular aggregation level not supported.")
  ####################################################
  #Transform the historical and projection datasets
  ####################################################
  if(!is.null(historical)){
    if(!is.list(historical)){
      tmp<-list()
      tmp[["historical"]]<-historical
      historical<-tmp
      rm(tmp)
    }
    if(length(historical)==1 & is.null(names(historical))) names(historical)<-"Historical"
  }
  #Get the required years if the years argument is not NULL
  if(!is.null(years)) historical<-lapply(X=historical,FUN=reduce_years,years)
  
  if(!is.null(projection)){
    if(!is.list(projection)){
      tmp<-list()
      tmp[["projection"]]<-projection
      projection<-tmp
      rm(tmp)
    }
    if(length(projection)==1 & is.null(names(projection))) names(projection)<-"Projection"
  }
  #Get the required years if the years argument is not NULL
  if(!is.null(years)) projection<-lapply(X=projection,FUN=reduce_years,years)
  
  data<-list()
  data[["historical"]]<-historical
  data[["projection"]]<-projection
  ####################################################
  #Calculate share of legend columns
  ####################################################  
  col1 <- length(magpie)
  col2 <- length(data[["historical"]])
  col3 <- length(data[["projection"]])
  col1 <- round_any(col1/5, 1, f = ceiling)
  col2 <- round_any(col2/5, 1, f = ceiling)
  col3 <- round_any(col3/5, 1, f = ceiling)
  allcol <- col1 + col2 + col3
  col1 <- col1/allcol
  col2 <- col2/allcol
  col3 <- col3/allcol
  if (col2 < 0.2 & col2 > 0) col2 <- 0.2
  if (col3 < 0.2 & col3 > 0) col3 <- 0.2
  col1 <- 1-col2-col3
  ####################################################
  #Convert to index if index=T
  ####################################################  
  if(index){
    index_year<-getYears(magpie[[1]])[1]
    reference<-setNames(setYears(magpie[[1]][,index_year,],NULL),NULL)
    for(i in 1:length(magpie)){
      magpie[[i]]<-magpie[[i]]/reference
    }
    index_transform<-function(x,reference){return(x/reference)}
    for(i in names(data)){
      data[[i]]<-lapply(data[[i]],index_transform,reference=reference)
    }
  }
  ####################################################
  #Prepare the data for plotting with ggplot
  ####################################################
  
  #Get the correct structure: 
  #-One Magpie object with all data points
  #-One magpie object with upper and lower errorbounds for datasets with nyears(data)==1
  #-One magpie object with upper and lower errorbounds for datasets with nyears(data)>1
  #All of them should have a row indicating historical/projection/magpie result
  #All of them should have a row indicating the scenario
  plotdata<-list()
  ploterrbar<-list()
  ploterrbar$up<-list()
  ploterrbar$lo<-list()
  ploterrbar$data<-list()
  ploterrshade<-list()
  ploterrshade$up<-list()
  ploterrshade$lo<-list()
  ploterrshade$data<-list()
  
  for(i in 1:length(magpie)){
    tmp<-as.data.frame(setNames(magpie[[i]],paste("MAgPIE",names(magpie)[i],sep=".")))
    names(tmp)[which(names(tmp)=="Data1")]<-"Historical"
    names(tmp)[which(names(tmp)=="Data2")]<-"Scenario"
    plotdata[[i]]<-tmp
  }
  for(hist in names(data)){
    for(sourc in names(data[[hist]])){
      #Get the datapoints
      tmp<-as.data.frame(setNames(data[[hist]][[sourc]][,,"data"],paste(hist,sourc,sep=".")))
      names(tmp)[which(names(tmp)=="Data1")]<-"Historical"
      names(tmp)[which(names(tmp)=="Data2")]<-"Scenario"
      plotdata[[length(plotdata)+1]]<-tmp
      #get the errors if available
      if(all(c("up","lo") %in% fulldim(data[[hist]][[sourc]])[[2]][[3]])){
        for(bound in c("up","lo","data")){
          tmp<-as.data.frame(setNames(data[[hist]][[sourc]][,,bound],paste(hist,sourc,sep=".")))
          names(tmp)[which(names(tmp)=="Data1")]<-"Historical"
          names(tmp)[which(names(tmp)=="Data2")]<-"Scenario"
          if(nyears(data[[hist]][[sourc]])==1){
            ploterrbar[[bound]][[length(ploterrbar[[bound]])+1]]<-tmp
          } else {
            ploterrshade[[bound]][[length(ploterrshade[[bound]])+1]]<-tmp
          }
        }
      }
    }
  }

  #Check what information is available
  errshade<-FALSE
  errbar<-FALSE
  if(length(ploterrshade$data)>0) errshade<-TRUE
  if(length(ploterrbar$data)>0) errbar<-TRUE
  
  #Combine all the datapoints into one dataframe
  plotdata<-suppressMessages(melt(plotdata,measure.vars="Value",value.name="Value"))
  plotdata$variable<-NULL
  plotdata$L1<-NULL
  plotdata$Cell<-factor(plotdata$Cell)
  plotdata$Region<-factor(plotdata$Region)
  plotdata$Value<-as.numeric(plotdata$Value)
  plotdata$Year <- as.Date(paste(plotdata$Year, "-01-01",sep = ""), format = "%Y-%m-%d")
# return(plotdata)
# }
  #Combine all the errorinformations for shaded error regions into one dataframe
  if(errshade){
    tmp<-suppressMessages(melt(ploterrshade[["up"]],measure.vars="Value",value.name="Value"))
    tmp$variable<-NULL
    tmp$L1<-NULL
    tmp$Cell<-factor(tmp$Cell)
    tmp$Region<-factor(tmp$Region)
    tmp$Value<-as.numeric(tmp$Value)
    names(tmp)[which(names(tmp)=="Value")]<-"up"
    tmp2<-suppressMessages(melt(ploterrshade[["lo"]],measure.vars="Value",value.name="Value"))
    tmp$lo<-as.numeric(tmp2$Value)
    tmp2<-suppressMessages(melt(ploterrshade[["data"]],measure.vars="Value",value.name="Value"))
    tmp$Value<-as.numeric(tmp2$Value)
    ploterrshade<-tmp
    ploterrshade$Year <- as.Date(paste(ploterrshade$Year, "-01-01",sep = ""), format = "%Y-%m-%d")
  }
  
  
  #Combine all the errorinformations for errorbars into one dataframe
  if(errbar){
    tmp<-suppressMessages(melt(ploterrbar[["up"]],measure.vars="Value",value.name="Value"))
    tmp$variable<-NULL
    tmp$L1<-NULL
    tmp$Cell<-factor(tmp$Cell)
    tmp$Region<-factor(tmp$Region)
    tmp$Value<-as.numeric(tmp$Value)
    names(tmp)[which(names(tmp)=="Value")]<-"up"
    tmp2<-suppressMessages(melt(ploterrbar[["lo"]],measure.vars="Value",value.name="Value"))
    tmp$lo<-as.numeric(tmp2$Value)
    tmp2<-suppressMessages(melt(ploterrbar[["data"]],measure.vars="Value",value.name="Value"))
    tmp$Value<-as.numeric(tmp2$Value)
    ploterrbar<-tmp
    ploterrbar$Year <- as.Date(paste(ploterrbar$Year, "-01-01",sep = ""), format = "%Y-%m-%d")
  }
# return(plotdata)
# }
  ######################################################################
  #Now finally do the plotting
  ######################################################################
  #Define some plotting styles
  linetypes<-list(MAgPIE=1,historical=1,projection=4)
  shapes<-list(MAgPIE=16,historical=1,projection=16)
  alpha_point<-list(MAgPIE=1,historical=1,projection=0.2)
  alpha_line<-list(MAgPIE=1,historical=0.3,projection=0.2)
  #Get the mapping between Scenario and Historical
  mapp<-vector()
  for(i in levels(plotdata$Scenario)){
    mapp<-c(mapp,unique(as.vector(plotdata$Historical[which(plotdata$Scenario==i)])))
    names(mapp)[length(mapp)]<-i
  }
  #Determine line and point sizes
  if(level %in% c("reg","SSP")){
    linewidth<-1
    pointwidth<-3
  } else{
    linewidth<-1.5
    pointwidth<-4
  }
  #First construct an empty plot
  p <- ggplot(plotdata, aes(x=Year, y=Value,fill=Scenario, colour=Scenario,alpha=Scenario,linetype=Scenario,shape=Scenario,group=Scenario))
  #Get the colors right
  #MAgPIE_colours<-c("#FF0000","#000000","#006666")
#  MAgPIE_colours <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33","#A65628", "#F781BF")
  set1 <- colorRampPalette(brewer.pal(9,"Set1"))  
  MAgPIE_colours <- set1(length(magpie))
  data_colours<-c("#0000FF","#006600","#FF00FF","#993300","#FF9900","#00CCCC")
  colours<-vector()
  historical_colours<-vector()
  projection_colours<-vector()
  magpies<-0
  datasets<-0
  for(i in mapp){
    if(i=="MAgPIE"){
      colours<-c(colours,MAgPIE_colours[magpies+1])
      magpies<-magpies+1
    } else {
      colours<-c(colours,data_colours[datasets+1])
      if(i=="historical"){
        historical_colours<-c(historical_colours,data_colours[datasets+1])
      } else{
        projection_colours<-c(projection_colours,data_colours[datasets+1])
      }
      datasets<-datasets+1
    }
  }
  drop<-TRUE
  p<-p+scale_fill_manual(values=colours,breaks=unique(plotdata$Scenario),drop=drop,name="Data type")
  p<-p+scale_colour_manual(values=colours,breaks=unique(plotdata$Scenario),drop=drop,name="Data type")
  #Get the correct linetypes
  p<-p+scale_linetype_manual(values=as.numeric(gsub("MAgPIE",linetypes[["MAgPIE"]],gsub("projection",linetypes[["projection"]],gsub("historical",linetypes[["historical"]],mapp)))),name="Data type",breaks=unique(plotdata$Scenario),drop=drop)
  #Get the correct pointshape
  p<-p+scale_shape_manual(values=as.numeric(gsub("MAgPIE",shapes[["MAgPIE"]],gsub("projection",shapes[["projection"]],gsub("historical",shapes[["historical"]],mapp)))),name="Data type",breaks=unique(plotdata$Scenario),drop=drop)
  #Get the correct alpha
  p<-p+scale_alpha_manual(values=as.numeric(gsub("MAgPIE",alpha_line[["MAgPIE"]],gsub("projection",alpha_line[["projection"]],gsub("historical",alpha_line[["historical"]],mapp)))),name="Data type",breaks=unique(plotdata$Scenario),drop=drop)
  #Add the lines
  p<-p+geom_line(size=linewidth)  
  #Add the points for historical values
  if(length(subset(plotdata,plotdata$Historical=="historical")$Value)!=0){
    p<-p+geom_point(data=subset(plotdata,plotdata$Historical=="historical"),size=pointwidth,alpha=1,show.legend=FALSE)
  }
  #Add the points for projections including MAgPIE projections
  if(length(subset(plotdata,plotdata$Historical!="historical")$Value)!=0){
    p<-p+geom_point(data=subset(plotdata,plotdata$Historical!="historical"),size=pointwidth)
  }  
  #Add Errorbars
  if(errbar){
    p<-p+geom_errorbar(data=ploterrbar,aes(ymax=up,ymin=lo,colour=Scenario),size=2*linewidth,alpha=0.15,show.legend=FALSE)
  }
  #Now add the shaded error regions
  if(errshade) {
    p<-p+geom_smooth(data=ploterrshade,aes(ymax=up,ymin=lo,fill=Scenario),stat="identity",colour=NA,alpha=0.15,size=0.1,show.legend=FALSE)
  }
  #Add a vertical line where the MAgPIE simulation starts
  p<-p+geom_vline(xintercept=as.numeric(min(subset(plotdata,plotdata$Historical=="MAgPIE")$Year)),linetype=2)
  #Do the splitting if regional outpout is tested
  if(level %in% c("reg","SSP")){
    if(same_yscale){
      scales<-"fixed"
    } else {
      scales<-"free_y"
    }
    p<-p+facet_wrap( ~ Region, ncol=5,scales=scales)
  }  
  #scale the axes
  if(!is.null(xlim)){
    xlim<-as.Date(paste(xlim,"-01-01",sep=""))
    p<-p+geom_vline(xintercept=as.numeric(xlim))
    if(!is.null(ylim)){
      p<-p+coord_cartesian(ylim=ylim,xlim=xlim)  
      p<-p+geom_hline(yintercept=ylim)
    } else {
      p<-p+coord_cartesian(xlim=xlim)  
    }
  } else if(!is.null(ylim)){
    p<-p+geom_hline(yintercept=ylim)
    p<-p+coord_cartesian(ylim=ylim)  
  }
  #Format the plot in a nice way
  p  <- p + labs(y=ylab) + ggtitle(title) + theme(panel.background = element_rect(fill='white', colour='black'),plot.title=element_text(size=text_size+4,face="bold",vjust=1.5), strip.text.x=element_text(size=text_size), axis.title.y=element_text(angle=90,size=text_size,face="bold",vjust=0.3), axis.text.y=element_text(size=text_size,colour="black"), axis.title.x=element_text(size=text_size,face="bold",vjust=-0.3), axis.text.x=element_text(size=text_size, angle=90, hjust=1.2,colour="black"),legend.position="none")
  
  ###################################################################
  #Construct the legend (unfortunately a bit complicated)
  ###################################################################
  leg<-list()
  legend_colours<-list(MAgPIE=MAgPIE_colours,historical=historical_colours,projection=projection_colours)
  lnames<-list(MAgPIE=legend_names[1],historical=legend_names[2],projection=legend_names[3])
  #Get one legend for each data type
  for(type in as.vector(unique(plotdata$Historical))){
    tmp<-ggplot(subset(plotdata,Historical==type), aes(x=Year, y=Value,fill=Scenario, colour=Scenario,group=Scenario))
    tmp<-tmp+scale_fill_manual(values=legend_colours[[type]][1:length(unique(subset(plotdata,Historical==type)$Scenario))],breaks=unique(subset(plotdata,Historical==type)$Scenario),drop=drop,name=lnames[[type]])
    tmp<-tmp+scale_colour_manual(values=legend_colours[[type]][1:length(unique(subset(plotdata,Historical==type)$Scenario))],breaks=unique(subset(plotdata,Historical==type)$Scenario),drop=drop,name=lnames[[type]])
    tmp<-tmp+geom_line(size=linewidth,linetype=linetypes[[type]],alpha=alpha_line[[type]])  
    tmp<-tmp+geom_point(size=pointwidth,shape=shapes[[type]],alpha=alpha_point[[type]]) 
    tmp<-tmp+theme(legend.title=element_text(size=text_size,face="bold"), legend.text=element_text(size=text_size),legend.background = element_rect(fill="white"),legend.key=element_blank())
    tmp<-tmp+guides(fill=guide_legend(nrow=5))
    tmp<-suppressWarnings(ggplot_gtable(ggplot_build(tmp)))
    leg[[type]]<-tmp$grobs[[8]]
  }
  rm(tmp)
  args<-leg
  args[["ncol"]]<-length(args)
  width<-vector()
  if("MAgPIE" %in% as.vector(unique(plotdata$Historical)))width<-c(width,col1)
  if("historical" %in% as.vector(unique(plotdata$Historical)))width<-c(width,col2)
  if("projection" %in% as.vector(unique(plotdata$Historical)))width<-c(width,col3)
  args[["widths"]]<-width
  if(!raw_legend){
    leg <- do.call(arrangeGrob,args=args)
  } else {
    combine_legend<-function(x){return(do.call(arrangeGrob,args=x))}
    leg <- args
    leg[["combine"]] <- combine_legend
  }
  
  ###########################################################################
  #Construct the final plot
  ###########################################################################
  combine<-function(graph, legend){
    if("combine" %in% names(legend)) legend<-legend$combine(legend[-which(names(legend)=="combine")])
    out<-arrangeGrob(graph,legend,ncol=1,heights=c(0.76,0.24))
    return(out)
  }
    
  if(!raw_plot){
    out<-combine(p,leg)
  } else {
    out<-list()
    out$figure<-p
    out$legend<-leg      
    out$combine<-function(x){out<-combine(x[[1]],x[[2]]);return(out)}
  }

  return(out)
}