R/dynatop.R

In waternumbers/dynatop: An Implementation of Dynamic TOPMODEL Hydrological Model in R

#' R6 Class for Dynamic TOPMODEL
#' @examples
#' ## the vignettes contains further details of the method calls.
#' 
#' data("Swindale") ## example data
#' mdl <- Swindale$model
#' mdl$map <- system.file("extdata","Swindale.tif",package="dynatop",mustWork=TRUE)
#' ctch_mdl <- dynatop$new(mdl$hru,map=mdl$map) ## create with model
#' ctch_mdl$add_data(Swindale$obs) ## add observations
#' ctch_mdl$initialise() ## initialise model
#' ctch_mdl$sim(Swindale$model$output_flux) ## simulate model
#' @export
dynatop <- R6Class(
            "dynatop",
    public = list(
        #' @description Creates a dynatop class object from the a list based model description as generated by dynatopGIS.
        #'
        #' @param model a dynamic TOPMODEL list object
        #' @param map file name of the map layers for the model
        #' @param use_states logical if states should be imported
        #' @param drop_map logical if the map should be dropped
        #' @param delta error term in checking redistribution sums
        #'
        #' @return invisible(self) suitable for chaining
        #'
        #' @details This function makes some basic consistency checks on a list representing a dynamic TOPMODEL model. The checks performed and basic 'sanity' checks. They do not check for the logic of the parameter values nor the consistency of states and parameters. Sums of the redistribution matrices are checked to be in the range 1 +/- delta.
        initialize = function(model, map=NULL, use_states=FALSE, delta = 1e-13){
            ## digest model with checks - will fail if errors identified
            private$digest_model(model,use_states,delta)
            if( !is.null( map ) ){
                if( !file.exists(map) ){ stop( "Map file does not exist") }
                private$digest_map(map)
            }
            invisible(self)
        },
        #' @description Adds observed data to a dynatop object
        #'
        #' @param obs_data an xts object of observed data
        #'
        #' @return invisible(self) suitable for chaining
        #'
        #' @details This function makes some basic consistency checks on the observations to ensure they have uniform timestep and all required series are present.
        add_data = function(obs_data){
            self$clear_data()
            ## check input and get model timestep
            private$digest_obs(obs_data)
            invisible(self)
        },
        #' @description Clears all forcing and simulation data except current states
        #'
        #' @return invisible(self) suitable for chaining
        clear_data = function(){
            private$time_series <- list()
        },
        #' @description Initialises a dynatop object in the simplest way possible.
        #'
        #' @param vtol tolerance for the solution for the saturated zone storage (as volume)
        #' @param ftol tolerance for the solution of the saturated zone storage (as difference of function from 0)
        #' @param max_it maximum number of iterations to use in the solution of the saturated zone
        #'
        #' @return invisible(self) suitable for chaining
        initialise = function(vtol = sqrt(.Machine$double.eps), ftol = sqrt(.Machine$double.eps), max_it = 1000){
            ## check the solver options
            vtol <- as.double(vtol); ftol <- as.double(ftol)
            if(any( c(vtol,ftol) < .Machine$double.eps)){
                stop("A solution tolerance is set lower then machine eps")
            }
            max_it <- as.integer(max_it)
            if(max_it < 10){stop("Please use at least 10 iterations")}

            private$init(vtol,ftol,max_it)
            invisible(self)
        },
        #' @description Simulate the hillslope and channel components of a dynatop object
        #' @param output_defn a description of the output series
        #' @param keep_states a vector of POSIXct objects (e.g. from xts) giving the time stamp at which the states should be kept
        #' @param sub_step simulation timestep in seconds, default value of NULL results in data time step
        #' @param vtol tolerance on width of bounds in the numeric search for surface and saturated zone solutions (as volume)
        #' @param ftol - not currently used
        #' @param max_it maximum number of iterations to use in the solution of the saturated zone
        #'
        #' @details Saving the states at every timestep and keeping the mass balance can generate very large data sets!!
        #'
        #' @return invisible(self) for chaining
        sim = function(output_defn,keep_states=NULL,sub_step=NULL,
                       vtol=0.001,ftol=sqrt(.Machine$double.eps), max_it=1000){
            
            ## check the solver options
            vtol <- as.double(vtol); ftol <- as.double(ftol)
            if(any( c(vtol,ftol) < .Machine$double.eps)){
                stop("A solution tolerance is set lower then machine eps")
            }
            max_it <- as.integer(max_it)
            if(max_it < 10){stop("Please use at least 10 iterations")}

            ## check digest output defn
            private$digest_output_defn(output_defn)
            
            ## check presence of finite states
            tmp <- sapply(private$model, function(x){ x$properties["area"]==0 | all(is.finite(x$states)) })
            if( !all(tmp) ){
                stop("Model states have non-finite values")
            }

            if( !is.null(sub_step) && !is.finite(sub_step[1]) ){
                stop("sub_step should be a single finite value")
            }
            
            ## check presense of obs
            if( length(private$time_series$index) < 2 ){
                stop("Insufficent data to perform a simulation")
            }
            
            ## check keep_states is valid
            if( length(keep_states)>0 ){
                if( !("POSIXct" %in% class(keep_states)) ){
                    stop("Times for returning states should be POSIXct object")
                }
            }
            keep_states <- keep_states[keep_states %in% private$time_series$index]

            #browser()
            private$sim_dyna(keep_states,sub_step,vtol,ftol,max_it)
            invisible(self)
        },
        ## ############
        ## Functions for extracting and plotting data
        #' @description Return channel inflow as an xts series or list of xts series
        #' @param name one or more output series to return
        get_output = function(name=colnames(private$time_series$output)){
            if( length(private$time_series$output) == 0 ){
                stop("No output data available")
            }
            
            name <- match.arg(name, colnames(private$time_series$output), several.ok=TRUE)
            x <- xts::xts(private$time_series$output[,name,drop=FALSE],
                          order.by=private$time_series$index)
            return(x)
        },
        #' @description Plot the channel inflow
        #' @param name of series to plot
        plot_output = function(name=colnames(private$time_series$output)){
            x <- self$get_output(name)
            plot(x)
            
            
            ## if(seperate){
            ##     oldpar <- par(no.readonly = TRUE)
            ##     on.exit(par(oldpar))
            ##     nc <- floor(sqrt(length(name)))
            ##     nr <- ceiling( length(name)/nc )
            ##     par(mfrow=c(nr,nc))
            ##     for(ii in name){
            ##         plot(x[,ii])
            ##     }
            ## }else{
            ##     plot(x)
            ## }
        },
        #' @description Get the observed data
        get_obs_data = function(){
            xts::xts(private$time_series$obs,
                     order.by=private$time_series$index)
        },
        #' @description Return the model
        get_model = function(){
            private$reform_model()
        },
        #' @description Return the model
        get_mass_errors = function(){
            ## if( !("mass_balance" %in% names(private$time_series)) ){
            ##     stop("Mass errors are not available")
            ## }
            xts::xts(private$time_series$mass_balance,
                     order.by=private$time_series$index)
        },
        #' @description Return states
        #' @param record logical TRUE if the record should be returned. Otherwise the current states returned
        get_states = function(record=FALSE){
            
            if( record ){
                return( setNames(private$time_series$state_record,
                                 private$time_series$index) )
            }else{
                tmp <- do.call(rbind, lapply(private$model, function(x){ c(id=x$id, x$states) }))
                return( as.data.frame(tmp) )
            }
       },
       #' @description Plot a current state of the system
       #' @param state the name of the state to be plotted
       # #' @param add_channel Logical indicating if the channel should be added to the plot
       plot_state = function(state=c("s_sf","s_rz","s_uz","s_sz")){ #,add_channel=TRUE){
           state <- match.arg(state)
           
           if( is.null(private$map) | ( length(private$map)==0) ){
               stop("The model contains no map of HRU locations")
           }

           x <- self$get_states()
           rst <- terra::subst(private$map[["hru"]], x$id, x[[state]])
           
           terra::plot(rst)
           ## if( add_channel & file.exists(private$model$map$channel) ){
           ##   chn <- terra::vect(private$model$map$channel)
           ##    terra::plot(chn,add=TRUE)
           ## }
           
       }
       
    ),
    private = list(
        ## stores of data
        version = "0.3.0",
        model = list(), # storage for model object
        map  = NULL, # storage for map object
        output_defn = list(), ## definition of output
        time_series = list(), ## storage for time series data
        info = list(sf = setNames(as.integer(1:3),c("cnst","kin","comp")),
                    rz = setNames(as.integer(1),c("orig")),
                    uz = setNames(as.integer(1),c("orig")),
                    sz = setNames(as.integer(1:4),c("exp","bexp","dexp","cnst")),
                    
                    output = setNames(1:14, c("precip","pet","aet",
                                              "q_sf","q_sf_in","q_sz","q_sz_in",
                                              "s_sf","s_rz","s_uz","s_sz",
                                              "v_sf_rz","v_rz_uz","v_uz_sz"))
                    ),
        digest_hru = function(h, use_states, delta){ ## check HRU returns a text string of errors
            etxt = character(0)
            
            ## check id
            if("id" %in% names(h)){
                if( length( h$id ) > 1 ){ etxt <- c(etxt, paste0(h$id[1], ": ID should be of length 1")) }
                if( !is.integer( h$id ) ){ etxt <- c(etxt, paste0( h$id, ": ID should be an integer")) }
            }else{
                etxt <- c(etxt, "No ID is specified")
                h$id <- NA ## for error messages
            }
            
            ## check properties
            if("properties" %in% names(h)){
                prpnm <- c("area", "width", "Dx", "gradient")
                if( !is.numeric(h$properties) ){ etxt <- c( etxt, paste0(h$id, ": properties should be a numeric vector") ) }
                if( all(prpnm %in% names(h$properties)) ){
                    if( !all( h$properties[c("gradient","width")] >0 ) ){
                        etxt <- c( etxt, paste0(h$id, ": gradient and width but be greater then 0") )
                    }
                    if( h$properties["area"] < 0 ){
                        etxt <- c( etxt, paste0(h$id, ": area must not be negative"))
                    }
                    h$properties <- h$properties[ c( prpnm, setdiff(names(h$properties),prpnm)) ]
                }else{
                    etxt <- c(etxt, paste0(h$id, ": properties is missing named values") )
                }
            }else{
                etxt <- c(etxt,paste0(h$id, ": properties is missing") )
            }
                
            ## check states
            snm <- c("s_sf","s_rz","s_uz","s_sz")
            if("states" %in% names(h)){
                if( !is.numeric(h$states) ){ etxt <- paste(etxt, paste0(h$id, ": states should be a numeric vector"), sep="\n") }
                if( !all(snm %in% names(h$states)) ){
                    etxt <- c(etxt, paste0(h$id, ": states is missing named values"))
                }
                h$states <- h$states[ c(snm, setdiff(names(h$states),snm)) ] ## make sure states are in correct order
            }else{
                etxt <- c(etxt,paste0(h$id, ": states is missing") )
            }
            
            ## check sf, rz, uz, sz
            for(ii in c("sf","rz","uz","sz")){
                if(!(ii %in% names(h))){
                    etxt <- c(etxt,paste0(h$id, ": ",ii," definition is not present"))
                    next
                }
                if( !all(c("type","parameters") %in% names(h[[ii]])) ){
                    etxt <- c(etxt,paste0(h$id, ": ",ii, " definition is missing type and.or parameters"))
                    next
                }
                if( length( h[[ii]]$type ) >1 ){
                    etxt <- c(etxt, paste0(h$id[1], ": ", ii, " type should be of length 1"))
                    next
                }
                
                if( !( h[[ii]]$type %in% names(private$info[[ii]])) ){
                    etxt <- c(etxt, paste0(h$id[1], ": ", ii, " type is not valid"))
                    next
                }
                pnm <- switch( paste0(ii, "_", h[[ii]]$type), ## make a unique code
                              "sf_cnst" = c("c_sf","d_sf","s_raf","t_raf"),
                              "sf_kin" = c("n","s_raf","t_raf"),
                              "sf_comp" = c("v_sf_1","d_sf_1","s_1","v_sf_2","d_sf_2"),
                              "rz_orig" = c("s_rzmax"),
                              "uz_orig" = c("t_d"),
                              "sz_exp" = c("t_0","m"),
                              "sz_bexp" = c("t_0","m","h_sz_max"),
                              "sz_cnst" = c("v_sz","h_sz_max"),
                              "sz_dexp" = c("t_0","m","m2","omega"),
                              stop("Invalid options for pname")
                              )
                if( !is.numeric( h[[ii]]$parameters )){
                    etxt <- c(etxt, paste0(h$id[1], ": ", ii, " parameters should be a numeric vector"))
                    next
                }
                if( !all( pnm %in% names(h[[ii]]$parameters)) ){
                    etxt <- c(etxt,paste0(h$id, ": ",ii, " is missing parameters"))
                    next
                }
                if( !all( h[[ii]]$parameters[ pnm ] >=0 ) ){
                    etxt <- c(etxt, paste0(h$id, ": some ", ii, " parameters are negatve"))
                    next
                }
                h[[ii]]$parameters <- h[[ii]]$parameters[ c(pnm,setdiff(names(h[[ii]]$parameters),pnm)) ] ## make sure parameters are in correct order
                ## print(h[[ii]]$parameters)
            }
            
            ## check precip and pet
            for(ii in c("precip","pet")){
                if( !all(c("name","fraction") %in% names(h[[ii]])) ){
                    etxt <- c(etxt, paste0(h$id, ": ", ii, " should contain names and fractions"))
                    next
                }
                if( !is.character(h[[ii]]$name) ){
                    etxt <- c(etxt, paste0(h$id, ": ", ii, " name should be a character vector"))
                    next
                }
                if( !is.numeric(h[[ii]]$fraction) ){
                    etxt <- c(etxt, paste0(h$id, ": ", ii, " fraction should be a numeric vector"))
                    next
                }
                if( length( h[[ii]]$name) != length(h[[ii]]$fraction) ){
                    etxt <- c(etxt, paste0(h$id, ": ", ii, " name and fraction should be the same length"))
                    next
                }
                if( any(h[[ii]]$fraction < 0) | ( abs( sum(h[[ii]]$fraction) -1) > delta) ){
                    etxt <- c(etxt, paste0(h$id, ": ", ii, " fractions should be positive and sum to 1"))
                    next
                }
            }

            ## check lateral flow
            for(ii in c("sf_flow_direction","sz_flow_direction")){
                if( !all(c("id","fraction") %in% names(h[[ii]])) ){
                    etxt <- c(etxt, paste0(h$id, ": ", ii, " should contain ids and fractions"))
                    next
                }
                if( !is.integer(h[[ii]]$id) ){
                    etxt <- c(etxt, paste0(h$id, ": ", ii, " id should be an integer vector"))
                    next
                }
                if( !is.numeric(h[[ii]]$fraction) ){
                    etxt <- c(etxt, paste0(h$id, ": ", ii, " fraction should be a numeric vector"))
                    next
                }
                if( length( h[[ii]]$id) != length(h[[ii]]$fraction) ){
                    etxt <- c(etxt, paste0(h$id, ": ", ii, " id and fraction should be the same length"))
                    next
                }
                if( any(h[[ii]]$id >= h$id) ){
                    etxt <- c(etxt, paste0(h$id, ": ", ii, " id value be less then current id"))
                    next
                }
                if( length(h[[ii]]$fraction)>0 ){
                    
                    if( any(h[[ii]]$fraction < 0) | ( abs( sum(h[[ii]]$fraction) -1) > delta) ){
                        etxt <- c(etxt, paste0(h$id, ": ", ii, " fractions should be positive and sum to 1"))
                        next
                    }
                }
            }

            ## fail if errors
            if( length( etxt ) >0 ){
                stop( paste(etxt,collapse = "\n") )
            }

            ## convert for C++
            for(ii in c("sf","rz","uz","sz")){ ## convert type to integer
                h[[ii]]$type <- private$info[[ii]][ h[[ii]]$type ]
            }

            if( !use_states ){ h$states[] <- NA }
            return(h)

        },
        regurge_hru = function(h){
            ## convert for C++
            for(ii in c("sf","rz","uz","sz")){ ## convert type to integer
                h[[ii]]$type <- private$info[[ii]][ h[[ii]]$type ]
            }
            return(h)
        },
        ## this code checks and digests the model
        digest_model = function(model, use_states, delta=1e-13){
            m <- lapply( model, private$digest_hru, use_states = use_states, delta = delta)
            ## check ids
            id <- sapply(m, function(x){x$id})
            idx <- order(id)
            id <- id[idx]
            if( !all( id == 0:(length(id)-1) ) ){ stop("ids are not in sequence") }
            
            private$model <- m[idx]
        },
        ## function to digest maps
        digest_map = function(mapFile){
            if( !requireNamespace("terra",quietly=TRUE) ){
                stop( "The terra package is required for plotting the maps - please install or add to libPath" )
            }
            private$map <- terra::rast(mapFile)
            if(!("hru" %in% names(private$map))){
                private$map <- NULL
                stop("Map should have a layer named hru")
            }
        },
        ## convert the form the internal storage to that input
        ## we presume the model has been checked!!
        reform_model = function(){            
            return( lapply( private$model, private$regurge_hru) )
        },        
        ## check and add observations
        digest_obs = function(obs){
            
            ## check types
            if(!is.xts(obs)){ stop("observations should be an xts object") }
            
            ## check constant time step
            tmp <- diff(as.numeric(index(obs)))
            if( !all( tmp == tmp[1] ) ){
                stop("Time steps in data are not unique")
            }

            ## get all the observed series names
            nm <- lapply(private$model,
                         function(h){unique(c(h$precip$name,h$pet$name))})
            nm <- unique(do.call(c,nm))

            ## check names
            idx <- nm %in% names(obs)
            if( !all(idx) ){
                stop(paste(c("Missing series:",nm[!idx]),collaspe=" "))
            }
            
            ## check all required values are finite
            if( !all(is.finite(obs[,nm])) ){
                stop("There are non finite values in the required time series")
            }

            nm = setNames(0:(ncol(obs)-1),colnames(obs))

            faddobs <- function(h,nm){
                h$precip$idx <- nm[ h$precip$name ]
                h$pet$idx <- nm[ h$pet$name ]
                h
            }
            
            private$model <- lapply( private$model, faddobs,nm = nm)
            private$time_series$obs_data <- as.matrix(obs)
            private$time_series$index <- index(obs)
            
        },
        ## digest the output definition
        digest_output_defn = function(defn){
            ## check table
            if( !is.data.frame(defn) ){ stop("Output definition should be a data frame") }
            if( !all(c("name","id","flux") %in% names(defn) ) ){ stop("Output definition must have variables name, id and flux") }
            if( !all(defn$id %in% (0:(length(private$model)-1))) ){
                stop(paste("id should be between 0 and",length(private$model)-1))
            }
            
            if( !("scale" %in% names(defn) ) ){
                warning("Output definition does not have scale - adding a vector of 1's")
                defn$scale <- 1
            }
            
            defn$name <- as.character(defn$name)
            defn$id <- as.integer(defn$id)
            defn$flux <- as.character(defn$flux)
            if( !all( defn$flux %in% names(private$info$output) ) ){
                stop("At least one flux type not recognised")
            }
            defn$scale <- as.numeric(defn$scale)
            unm <- unique(defn$name)
            defn$name_idx <- setNames(0:(length(unm)-1),unm)[ defn$name ]
            defn$flux_int <- private$info$output[ defn$flux ]
            private$output_defn <- defn
            private$time_series$output <- matrix(as.numeric(NA), length(private$time_series$index), length(unm))
            colnames( private$time_series$output ) <- unm
        },
        ## reform the output definition if required
        reform_output_defn = function(){
            defn <- private$output_defn
            defn$flux_int <- defn$name_idx <- NULL
            return( defn )
        },    
        ## compute the simulation timestep
        comp_ts = function(sub_step=NULL){
            ## work out time steps for use in simulation
            ts <- list()
            ts$step <- diff( as.numeric(private$time_series$index[1:2])) # seconds
            if(is.null(sub_step)){sub_step <- ts$step}
            ts$n_sub_step <- max(1,floor(ts$step/sub_step)) # dimensionless
            ts$sub_step <- ts$step / ts$n_sub_step
            private$info$ts <- ts
        },
        ## ###########################################
        ## Initialise the states
        init = function(vtol,etol,max_it){
            ##TODO we could initialise without obs - maybe change to allow this
            if( length(private$model[[1]]$precip$idx)==0 ){
                stop("Please add data before initialisation")
            }
            dt_init(private$model,
                    vtol,etol,max_it)
            
        },
        ## ###############################
        ## function to perform simulations
        sim_dyna= function(keep_states,sub_step,vtol,etol,max_it){
            
            ## compute time substep
            if(length(sub_step)>1){ sub_step <- sub_step[1] }
            ts <- private$comp_ts(sub_step)            

            ## Logical if states to be kept and store
            keep_states <- private$time_series$index %in% keep_states
            if(any(keep_states)){
                private$time_series$state_record <- rep(list(as.data.frame(NULL)),length(private$time_series$index))
            }else{
                private$time_series$state_record <- list()
            }
            
            ## Initialise the mass error store
            private$time_series$mass_balance <- matrix(as.numeric(NA),nrow(private$time_series$obs),6)
            colnames(private$time_series$mass_balance) <-
                c("initial_state","p","e_t","outflow","final_state","error")
            
            ## simulate
            dt_sim(private$model,
                   private$output_defn,
                   as.logical( keep_states ),
                   private$time_series$obs,
                   private$time_series$mass_balance,
                   private$time_series$output,
                   private$time_series$state_record,
                   as.double(ts$step),
                   ts$n_sub_step,
                   as.double(vtol),
                   as.double(etol),
                   as.integer(max_it))            
            
        }
        
    )
    
)