R/frcstPlot.R

In waternumbers/FloodForT: Flood Forecasting Toolbox

#' FloodForT REAL TIME
#' Generate the hydrograph and exceedence probability plots for a given set of forecast data
#'
#' @param frcst forecast data structure of quantiles
#' @param frcstTime time frcst issued - POSIXct object must be in index of frcst
#' @param lvl numeric vector containing named values ''alert'' and ''warning'' 
#' @param stt Station identifier - used both to label the plots and in the filename
#' @param outFile String specifying the file to generate the plot in
#' @param errMsg String containing any error message to add to the plot
#' 
#' @return NULL
#' @export
#
frcstPlot <- function(frcst,frcstTime,lvl,stt,outFile=NA,errMsg=NULL){

    ## extract the times from the xts object
    ts <- index(frcst) # POSIXct
    dts <- as.numeric(ts[2])-as.numeric(ts[1])  # in seconds

    ## determine timestep forecast is issued
    idx <- which(ts %in% frcstTime)
    npast <- idx
    ens <- names(frcst)
    ens <- ens[!(ens %in% 'obs')]
    
    ## make probability of exceedances
    nf <- nrow(frcst)-idx
    prob <- matrix(NA,nf,2)
    rownames(prob) <- paste(1:nf)
    
    for(ii in 1:nf){
        ## if( is.finite(frcst[idx+ii,'obs']) ){
        ##     prob[ii,1] <- as.numeric( frcst[idx+ii,'obs']>= lvl['danger'] )
        ##     prob[ii,2] <- as.numeric( frcst[idx+ii,'obs']>= lvl['warning'] )
        ## }else{
        prob[ii,1] <- 1 - (sum( lvl['danger'] >= frcst[idx+ii,ens] )/(length(ens)+1))
        prob[ii,2] <- 1 - (sum( lvl['warning'] >= frcst[idx+ii,ens])/(length(ens)+1))
        ## }
    }


    ## specify the quantiles and colour scheme
    nquanti <- floor(length(ens)/2)
    colseq <- rgb(seq(0.97,0.1,length=nquanti),seq(1,0.6,length=nquanti),1)

    ## set the plot limits
    ylm <- range(frcst,na.rm=TRUE)
    ylm[1] <- 0.5 * max(0 , min( c(ylm[1],lvl) , na.rm=TRUE ) , na.rm=TRUE)
    ylm[2] <- 1.1 * max( c(ylm[2],lvl) , na.rm=TRUE )

    xlm <- as.numeric(range(ts))
    xlm <- c(xlm[1]-dts/2,xlm[2]+3*dts/2)

    ## initialise the plot and layout
    if(is.na(outFile) | is.null(outFile)){
#        x11()
    }else{
         svg(filename = outFile,
             width = 8, height = 5.6,
             pointsize = 12, bg = "white", 
             antialias = "subpixel")
     }
    
    nf <- layout(matrix(c(1,1,1,2,1,1,1,3),2,4,byrow=TRUE))
    
    ## main plot
    plot.new()
    par(mar=c(5,5.5,4,2)+0.1)
    plot.window(xlm,ylm)
    
    ## plot the data
    for(ii in 1:nrow(frcst)){
    
        for( jj in 1:nquanti ){
            xx <- c(ts[ii]-dts,ts[ii])
            yy <- as.numeric(frcst[ii,c(ens[jj],ens[length(ens)+1-jj])])
            xx <- rep(xx,each=2)
            yy <- c(yy,rev(yy))
            if(jj == 1){
                polygon(xx, yy, col = colseq[jj],
                       ,lty=1,lwd = 0.3)
            }else{
                 polygon(xx, yy, col = colseq[jj],
                         border = NA)
             }
        }
        
        if( is.finite(frcst[ii,'obs']) ){
            xx <- c(ts[ii]-dts,ts[ii])
            yy <- c(frcst[ii,'obs'],frcst[ii,'obs'])
            lines(xx,yy,col=1,lwd=2)
        }
    
    }


    ## plot the levels
    abline(h=lvl['danger'],col=2,lwd=2)
    text(xlm[1]+2*dts,lvl['danger'],"Danger",pos=3,cex=2)
    
    abline(h=lvl['warning'],col=8,lty=3,lwd=2)
    text(xlm[1]+2*dts,lvl['warning'],"Warning",pos=3,cex=2)

    ## plot the forecast time
    abline(v = frcstTime, lty=2)
    
    ## add the legend
    legend(xlm[1]+floor(npast/4)*dts,ylm[2]*1.05,
           legend=c("Observation",
               "Predictive Uncertainty"),
           col=c(1,colseq[4]),lty=c(1,-1),bty="n",
           lwd=c(2,-1),pch=c(-1,-1),
           fill = c(0,colseq[4]),
           border = c(0,colseq[4]),
           cex=2)

    ## add axes
    axis(1,at=ts-dts/2,
         labels= (difftime(ts,frcstTime,'secs')/dts),
         cex.axis=2)
    axis(2,cex.axis=2)
    axis(3,at=as.numeric(trunc(ts,'hours')),
         labels=format(trunc(ts,'hours'),"%Y-%m-%d %H:%M"),
         cex.axis=2)
    
    ## add title and axes labels
    title(xlab="Leadtime",
          ylab = 'Water Level [m]',
          cex.lab=2)

    ## add error messages
    if(!is.null(errMsg)){
        text(xlm[2]+4*dts,0.5*(ylm[1] + ylm[2]),
             errMsg)
    }

    ## ##################
    ## first subplot
    barplot(prob[,2],
            space=0,
            xlab="Leadtime",
            ylab="Probability",
            main="Probability of \n Warning",
            ylim=c(0,1),
            col = colseq[length(colseq)],
            cex.main=2,
            cex.axis=2,
            cex.names=2,
            cex.lab=2)
    abline(h=0)


    ## second subplot
    barplot(prob[,1],
            space=0,
            xlab="Leadtime",
            ylab="Probability",
            main="Probability of \n Danger",
            ylim=c(0,1),
            col = colseq[length(colseq)],
            cex.main=2,
            cex.axis=2,
            cex.names=2,
            cex.lab=2)
    abline(h=0)

    ## shut the plot
    if(!is.na(outFile) & !is.null(outFile)){
        dev.off()
    }

    return(NULL)
}