inst/ModelBuild/svrDAEstimate.R

In waternumbers/FloodForT: Flood Forecasting Toolbox

## ####################################
## server code for Estimate Data Assimilation tab
## exception is the reactive function to
## estimate the models which is in
## svrDAEstimateModel.R
## ####################################

## ######################################
## table of results


## plots of hydrographs
## plots of hydrographs
output$plotEstDA <- renderDygraph({
    ## check if any names
    nms <- input$selDAMdl
    if( is.null(nms) ){
        return(NULL)
    }

    ## check if any horizons
    hrz <- input$selDAHorizon
    if( is.null(hrz) ){
        return(NULL)
    }


    ## make observed data
    if(input$ckValDataDA){
        idx <- paste(
            format(analysisRecord$mdlData$Periods['Valid','start'],tz='GMT'),
            format(analysisRecord$mdlData$Periods['Valid','finish'],tz='GMT'),
            sep="::")
    }else{
        idx <- paste(
            format(analysisRecord$mdlData$Periods['Calib','start'],tz='GMT'),
            format(analysisRecord$mdlData$Periods['Calib','finish'],tz='GMT'),
            sep="::")
    }
    tmp <- analysisRecord$mdlData$Variables['output']
    X <- analysisRecord$baseData[idx,tmp]

    ## handle the case of no observations
    if(all(!is.finite(X))){
        return(NULL)
    }

    for(ii in nms){
        tmp <- analysisRecord$mdl[[ii]]
        if(input$ckValDataDA){
            tmp <- tmp$val
        }else{
             tmp <- tmp$cal
         }
        for(jj in hrz){
            ttmp <- tmp[,jj]
            names(ttmp) <- paste(ii,jj)
            if(substr(jj,1,1)=="f"){
                index(ttmp) <- index(ttmp)+ as.numeric(substr(jj,3,nchar(jj))) * analysisRecord$mdl[[ii]]$param$mdl[1,'dt']
            }
            X <- cbind(X,ttmp)
        }
    }

    return( dygraph(X) %>% 
               dyOptions(useDataTimezone = TRUE) )# %>%
                #    dySeries(mtmp) )
})


## plots of thresholded CRPS
output$plotTcrpsDA <- renderPlot(
    {
    ## check if any names
    nms <- input$selDAMdl
    if( is.null(nms) ){
        return(NULL)
    }

    ## check if any horizons
    hrz <- input$selDAHorizon
    if( is.null(hrz) ){
        return(NULL)
    }

    ## make array of values
    y <- NULL
    nmsy <- NULL
    cnt <- 0
    for(ii in nms){
        tmp <- analysisRecord$mdl[[ii]]
        if(input$ckValDataDA){
            crps <- tmp$val.perf$crps
            vcrps <- tmp$val.perf$vcrps
        }else{
             crps <- tmp$cal.perf$crps
             vcrps <- tmp$cal.perf$vcrps
         }

        ## if values don't exist
        if(is.null(crps) |is.null(vcrps)){
            next
        }


        for(jj in hrz){
            if(jj %in% names(crps)){
                
                if(is.null(y)){
                    y <- array(NA,c(nrow(crps),3,length(nms)*length(hrz)))
                }
                nmsy <- c(nmsy,paste(ii,jj))
                cnt <- cnt + 1
                y[,1,cnt] <- sqrt( crps[,jj] - 2*sqrt(vcrps[,jj]) )
                y[,2,cnt] <- sqrt( crps[,jj] )
                y[,3,cnt] <- sqrt( crps[,jj] + 2*sqrt(vcrps[,jj]) )
            }
        }
    }
    
    ## if nothing to plot
    if(is.null(y)){return(NULL)}

    
    y <- y[,,1:cnt,drop=FALSE]
    nmsy <- nmsy[1:cnt]
    dimnames(y) <- list(NULL,NULL,nmsy)

    yrng <- c(0.99*min(y,na.rm=TRUE),1.01*max(y,na.rm=TRUE))
    lty <- rep(1,length(nmsy))
    names(lty) <- nmsy
    col <- 1:length(nmsy)
    names(col) <- nmsy

    matplot(crps[,'thold'],y[,2,],
            type="l",
            lty=lty,
            col=col,
            ylim=yrng,
            xlab="Threshold",
            ylab="CRPS"            
    )

    for(ii in nmsy){
        for(jj in 1:(dim(y)[1])){
            segments(crps[jj,'thold'], y[jj,1,ii], y1 = y[jj,3,ii],
                     col = col[ii],lty=lty[ii])
        }
    }

    legendStr <- nmsy
    legend("topleft",legendStr,
           lty = lty, lwd = 1,
           pch = NULL,
           col = col)
})
## plots of thresholded RMSE
output$plotTrmseDA <- renderPlot(
    {
    ## check if any names
    nms <- input$selDAMdl
    if( is.null(nms) ){
        return(NULL)
    }

    ## check if any horizons
    hrz <- input$selDAHorizon
    if( is.null(hrz) ){
        return(NULL)
    }

    ## make array of values
    y <- NULL
    nmsy <- NULL
    cnt <- 0
    for(ii in nms){
        tmp <- analysisRecord$mdl[[ii]]
        if(input$ckValDataDA){
            mse <- tmp$val.perf$mse
            vmse <- tmp$val.perf$vmse
        }else{
             mse <- tmp$cal.perf$mse
             vmse <- tmp$cal.perf$vmse
         }

        ## if values don't exist
        if(is.null(mse) |is.null(vmse)){
            next
        }

        for(jj in hrz){
            if(jj %in% names(mse)){
                
                if(is.null(y)){
                    y <- array(NA,c(nrow(mse),3,length(nms)*length(hrz)))
                }
                nmsy <- c(nmsy,paste(ii,jj))
                cnt <- cnt + 1
                y[,1,cnt] <- sqrt( mse[,jj] - 2*sqrt(vmse[,jj]) )
                y[,2,cnt] <- sqrt( mse[,jj] )
                y[,3,cnt] <- sqrt( mse[,jj] + 2*sqrt(vmse[,jj]) )
            }
        }
    }

    ## if nothing to plot
    if(is.null(y)){ return(NULL)}

    y <- y[,,1:cnt,drop=FALSE]
    nmsy <- nmsy[1:cnt]
    dimnames(y) <- list(NULL,NULL,nmsy)

    yrng <- c(0.99*min(y,na.rm=TRUE),1.01*max(y,na.rm=TRUE))
    lty <- rep(1,length(nmsy))
    names(lty) <- nmsy
    col <- 1:length(nmsy)
    names(col) <- nmsy

    matplot(mse[,'thold'],y[,2,],
            type="l",
            lty=lty,
            col=col,
            ylim=yrng,
            xlab="Threshold",
            ylab="RMSE"            
    )

    for(ii in nmsy){
        for(jj in 1:(dim(y)[1])){
            segments(mse[jj,'thold'], y[jj,1,ii], y1 = y[jj,3,ii],
                     col = col[ii],lty=lty[ii])
        }
    }

    legendStr <- nmsy
    legend("topleft",legendStr,
           lty = lty, lwd = 1,
           pch = NULL,
           col = col)
})



############################################
## plot of error models
output$plotErrDA <- renderPlot(
    {
    ## check if any names
    nms <- input$selDAMdl
    if( is.null(nms) ){
        return(NULL)
    }

    ## check if any horizons
    hrz <- input$selDAHorizon
    if( is.null(hrz) ){
        return(NULL)
    }


    ## make observed data
    if(input$ckValDataDA){
        idx <- paste(
            format(analysisRecord$mdlData$Periods['Valid','start'],tz='GMT'),
            format(analysisRecord$mdlData$Periods['Valid','finish'],tz='GMT'),
            sep="::")
    }else{
        idx <- paste(
            format(analysisRecord$mdlData$Periods['Calib','start'],tz='GMT'),
            format(analysisRecord$mdlData$Periods['Calib','finish'],tz='GMT'),
            sep="::")
    }
    tmp <- analysisRecord$mdlData$Variables['output']

    X <- analysisRecord$baseData[idx,tmp] ## will contain forecasts
    E <- analysisRecord$baseData[idx,tmp] ## will contain errors

    n <- length(nms)*length(hrz)
    legendStr <- NULL

    ## loop to get data for the scatter plot
    for(ii in nms){
        tmp <- analysisRecord$mdl[[ii]]
        if(input$ckValDataDA){
            yhat <- tmp$val
        }else{
             yhat <- tmp$cal
         }
        for(jj in hrz){
            if(!(jj %in% names(yhat))){
                next
            }
            
            tmp <- yhat[,jj]
            if(substr(jj,1,1)=="f"){
                index(tmp) <- index(tmp)+ as.numeric(substr(jj,3,nchar(jj))) * analysisRecord$mdl[[ii]]$param$mdl[1,'dt']
            }
            X <- cbind(X,tmp)
            E <- cbind(E,X[,1]-tmp)
            legendStr <- c(legendStr,paste(ii,jj))
        }
    }

    matplot(X[,-1,drop=FALSE],E[,-1,drop=FALSE],
            col=1:(ncol(X)-1),type="p",pch=1,
            xlab="Forecast",
            ylab="Error")

    ## add lines
    cnt <- 0
    for(ii in nms){
        tmp <- analysisRecord$mdl[[ii]]
        for(jj in hrz){
            if(!(jj %in% names(tmp$param))){
                next
            }

            ttmp <- tmp$param[[jj]]
            idx <- paste(ttmp$q) %in% paste((1:9)/10)

            cnt <- cnt + 1
            matlines(ttmp$x,ttmp$y[,idx],col=cnt,lty=1)
        }
    }
}
)
############################################
## plot of error models

output$plotRealTime <- renderPlot(
    {
    ## check if any names
    nms <- input$selDAMdl
    if( is.null(nms) ){
        return(NULL)
    }
    ## get the first model
    nms <- nms[1] # only use the first model
    mdl <- analysisRecord$mdl[[nms]]
    
    ## get the time the forcast is to issued
    time.index <- input$sldDAmovTime

    ## make observed data
    if(input$ckValDataDA){
        idx <- paste(
            format(analysisRecord$mdlData$Periods['Valid','start'],tz='GMT'),
            format(analysisRecord$mdlData$Periods['Valid','finish'],tz='GMT'),
            sep="::")
        Y <- mdl$val
    }else{
        idx <- paste(
            format(analysisRecord$mdlData$Periods['Calib','start'],tz='GMT'),
            format(analysisRecord$mdlData$Periods['Calib','finish'],tz='GMT'),
            sep="::")
        Y <- mdl$cal
    }
    tmp <- analysisRecord$mdlData$Variables['output']
    X <- analysisRecord$baseData[idx,tmp] # all obs data
    names(X) <- 'output'

    ## make vector of time stamps
    dt.issued <- index(Y)[time.index]
    dts <- seq(dt.issued - 12*mdl$param$mdl[1,'dt'],
               dt.issued + mdl$param$mdl[1,'lag']*mdl$param$mdl[1,'dt'],
               by = mdl$param$mdl[1,'dt'])
    
    ## make the matrix of data
    nq <- length(mdl$param[['f.0']]$q)
    frcst <- matrix(NA,length(dts),nq+1)
    colnames(frcst) <- c('obs',paste0("q",mdl$param[['f.0']]$q))

    for(ii in 1:length(dts)){
        jj <- which(index(X)==dts[ii])
        if(length(jj)>0){
            frcst[ii,'obs'] <- X[jj,'output']
        }
        
        ## find index of nearest timestep with an observation
        if(dts[ii] > dt.issued){
            jj <- which(index(Y)==dt.issued)
            ff <- (as.numeric(dts[ii]) - as.numeric(dt.issued))/mdl$param$mdl[1,'dt']
        }else{
             jj <- which(index(Y)==dts[ii])
             ff <- 0
         }

        if(length(jj)==0){next}

        str <- paste("f",ff,sep=".")
        es <- mdl$param[[str]]
        for(kk in 1:ncol(es$y)){
            frcst[ii,1+kk] <- Y[jj,str] + approx(es$x,es$y[,kk],xout=Y[jj,str],rule=2)$y
        }
    }

    frcst <- as.xts(frcst,order.by=dts)

    frcstPlot(frcst,dt.issued,mdl$param$lvl,'outBase',NA,FALSE)
})