DEPONS2R: Read, Plot and Analyse Output from the DEPONS Model

Documented in read.DeponsBlockdyn

# Author: Jacob Nabe-Nielsen
# Date: 2 September 2020
# Licence GPL v3
# Description: Methods and classes for reading and summarizing DEPONS output
#    stored in 'PorpoisePerBlock'-files, i.e. population dynamics for different
#    parts of the landscape.


#' @title  DeponsBlockdyn-class
#' @description Stores objects containing population size for different parts of
#' the landscape (i.e. different 'blocks')
#' @slot title Character. Name of the object or simulation
#' @slot landscape Character. Identifier for the landscape used in the DEPONS
#' simulations. The landscapes 'DanTysk', 'Gemini', 'Kattegat', 'North Sea',
#' 'Homogeneous', and 'User defined' are distributed with the DEPONS model.
#' @slot simtime \code{\link{POSIXlt}} object with the date and time when the simulation was
#' finished. This is read from the name of the imput file.
#' @slot startday POSIXlt object with the first day of the simulation, i.e.
#' the first day in the period that the simulations are intended to represent in
#' the real world.
#' @slot dyn Data frame with simulation output.
#' @details The \code{dyn} slot contains a data frame with the columns 'tick',
#' which indicates the number of half-hourly time steps since the start of the
#' simulation; a column 'block' indicating the region of the landscape where
#' animals were counted, a 'count' column with the number of animals in that block
#' and tick. The 'real.time' column shows the real-world equivalent to 'tick, i.e.
#' the time that has passed since 'startday'.
#' @exportClass DeponsBlockdyn
#' @examples a.DeponsBlockdyn <- new("DeponsBlockdyn")
#' a.DeponsBlockdyn
#' @note DeponsBlockdyn-objects are usually read in from csv files produced during
#' DEPONS simulations. These files are named 'PorpoisePerBlock.XXX.csv', where XXX
#' indicates the date and time when the simulation was finished.
#' @seealso \code{\link[DEPONS2R]{plot.DeponsBlockdyn}} and
#' \code{\link[DEPONS2R]{read.DeponsBlockdyn}}.
setClass(Class="DeponsBlockdyn",
         slots=list(title="character", landscape="character", simtime="POSIXlt",
                    startday="POSIXlt", dyn="data.frame")
)


setMethod("initialize", "DeponsBlockdyn",
          function(.Object) {
            .Object@title <- "NA"
            .Object@landscape <- "NA"
            .Object@simtime <- as.POSIXlt(NA)
            .Object@startday <- as.POSIXlt(NA)
            .Object@dyn <- data.frame("tick"=NA, "count.X"=NA, "real.time"=NA)
            return((.Object))
          }
)

#' @name Summary-methods
#' @aliases summary
#' @aliases summary,DeponsBlockdyn-method
#' @rdname summary
#' @title Summary
#' @description Summarizes different kinds of objects created based on output
#'   from the DEPONS model
#' @param object Depons* object
#' @details The summary method is available for \code{\link{DeponsTrack-class}},
#' \code{\link{DeponsDyn-class}}, \code{\link{DeponsRaster-class}},
#' and \code{\link{DeponsBlockdyn-class}}-objects.
#' @return list summarizing the DeponsBlockdyn object
#' @exportMethod summary
setMethod("summary", "DeponsBlockdyn",
          function(object) {
            cat("class:    \t", "DeponsBlockdyn \n")
            cat("title:    \t", object@title, "\n")
            cat("landscape:\t", object@landscape, "\n")
            cat("simtime:  \t", as.character(format(object@simtime)), "\n")
            cat("startday:  \t", as.character(object@startday), "\n")
            cat("n.ticks:  \t", as.character(max(object@dyn$tick)), "\n")
            cat("n.days:  \t", as.character(round(max(object@dyn$tick)/48, 2)), "\n\n")
            b.nos <- sort(unique(object@dyn$block))
            rnd <- function(n) sprintf(n, fmt='%#.2f')
            cat("Count per block:\n")
            cat("     block \t min   \t mean \t\t max \n" )
            for (b in b.nos) {
              cnt <- object@dyn$count[object@dyn$block==b]
              if(max(cnt)<1000) cat("   \t", b, "\t", min(cnt), "   \t", rnd(mean(cnt)), "   \t", max(cnt), "\n")
              else cat("  \t", b, "\t", min(cnt), "\t", rnd(mean(cnt)), "\t", max(cnt), "\n")
            }
            nticks <- max(object@dyn$tick)
            ndays <- max(object@dyn$tick)/48
            out <- list(
              "title" <- object@title,
              "landscape" <- object@landscape,
              "simtime" <- object@simtime,
              "startday" <- object@startday,
              "nticks" <- nticks,
              "ndays" <- object@ext
            )
            return(invisible(out))
          }
)



#' @title Reading simulated population count for blocks
#' @description Function for reading DEPONS simulation output with number of
#' animals per block for each time step.
#'
#' @param fname Name of the file (character) that contains movement data
#' generated by DEPONS. The name includes the path to the directory if this is
#' not the current working directory.
#' @param title Optional character string giving name of simulation
#' @param landscape The landscape used in the simulation
#' @param simtime Optional text string with date of simulation (format:
#' yyyy-mm-dd). If not provided this is obtained from name of input file
#' @param startday The start of the period that the  simulation represents, i.e.
#' the real-world equivalent of 'tick 1' (POSIXlt)
#' @seealso See \code{\link{DeponsBlockdyn-class}} for details on what is stored in
#' the output object and \code{\link{read.DeponsParam}} for reading the parameters
#' used in the simulation.
#' @export read.DeponsBlockdyn
#' @return \code{DeponsBlockdyn} object
#' @examples
#' \dontrun{
#' # File loaded from default location
#' the.file <- "/Applications/DEPONS 2.1/DEPONS/PorpoisePerBlock.2020.Sep.02.20_24_17.csv"
#' file.exists(the.file)
#' porpoise.blockdyn <- read.DeponsBlockdyn(fname=the.file,
#'   title="Test simulation with two blocks", landscape="North Sea")
#' porpoise.blockdyn
#'
#' # Get the latest simulation
#' the.file <- get.latest.sim(type="blockdyn", dir="/Applications/DEPONS 2.1/DEPONS")
#' owd <- getwd()
#' setwd("/Applications/DEPONS 2.1/DEPONS")
#' porpoise.blockdyn <- read.DeponsBlockdyn(fname=the.file)
#' setwd(owd)
#' }
read.DeponsBlockdyn <- function(fname, title="NA", landscape="NA", simtime="NA",
                           startday="NA") {
  raw.data <- utils::read.csv(fname, sep=";")
  # Get sim date and time from file name
  if (simtime=="NA")  simtime <- get.simtime(fname)
  if (startday=="NA")  startday <- NA
  all.data <- new("DeponsBlockdyn")
  all.data@title <- title
  all.data@landscape <- landscape
  all.data@simtime <- as.POSIXlt(simtime)
  all.data@startday <- as.POSIXlt(startday)
  the.data <- utils::read.csv(fname, sep=",")
  names(the.data) <- c("tick", "block", "count")
  # Add "real time" if startday is not NA
  tick.1.secs <- as.numeric(tick.to.time(1))
  secs.since.start <- the.data$tick-tick.1.secs
  the.data$real.time <- as.POSIXct(startday)+secs.since.start
  all.data@dyn <- the.data
  return(all.data)
}



#' @title Plot a DeponsBlockdyn object
#' @description Plot population dynamics simulated with DEPONS
#' @aliases plot.DeponsBlockdyn
#' @param x DeponsBlockdyn object
#' @param y Not used
#' @param dilute Integer. Plot only one in every 'dilute' values. Defaults to
#' 5, which yields a plot of the first simulated value and one in every five of
#' the following values.
#' @param ... Optional plotting parameters
#' @examples data("porpoisebdyn")
#' my.col <- c("red", "darkgreen", "orange")
#' plot(porpoisebdyn, col=my.col)
#' legend("bottomright", bty="n", fill=my.col, legend=paste("Block", 0:2))
#'
#' # Show all data points for small range of x-values
#' plot(porpoisebdyn, xlim=c(1950, 2050), ylim=c(4850, 5050), type="p", dilute=1, col=my.col)
##' @importFrom graphics points
#' @return \code{data.frame} listing blocks where no animals were counted
#' (returned invisibly)
#' @note The function returns a data frame with numbers of blocks with no agents.
setMethod("plot", signature("DeponsBlockdyn", "missing"),
          function(x, y, dilute=5, ...)  {
            oldpar <- graphics::par(no.readonly = TRUE)
            on.exit(graphics::par(oldpar))
            all.messages <- data.frame()
            if (!(dilute %% 1 == 0)) stop("'dilute' must be an integer")
            use.row <- x@dyn$tick %% dilute == 0
            z <- x@dyn[use.row ,] # use only every X rows
            use.row[1] <- TRUE
            graphics::par(mar=c(4.2, 4.2, 4, 4.2))
            if(!hasArg(pch)) {
              pch <- 16
            } else {
              pch <- list(...)[["pch"]]
            }
            if(!hasArg(cex)) {
              cex <- 0.8
            } else {
              cex <- list(...)[["cex"]]
            }
            if(!hasArg(ylab)) {
              ylab <- "count"
            } else {
              ylab <- as.character(list(...)[["ylab"]])
            }
            if(!hasArg(lwd)) {
              lwd <- 2
            } else {
              lwd <- as.character(list(...)[["lwd"]])
            }
            if(!hasArg(lty)) {
              lty <- 1
            } else {
              lty <- as.character(list(...)[["lty"]])
            }
            if(!hasArg(col)) {
              col <- "blue"
              col <- grDevices::rainbow(length(unique(x@dyn$block)))
            } else {
              col <- as.character(list(...)[["col"]])
              if(length(col) != length(unique(x@dyn$block))) stop("The nunber of colours must match the number of blocks in the landscape")
            }
            if(!hasArg(type)) {
              type <- "l"
            } else {
              type <- as.character(list(...)[["type"]])
            }
            if(!hasArg(axes)) {
              axes <- TRUE
            } else {
              axes <- list(...)[["axes"]]
            }
            if(!hasArg("ylim")) {
              ylim <- c(min(x@dyn$count), max(x@dyn$count))
            } else {
              ylim <- list(...)[["ylim"]]
            }
            if(!hasArg("main")) {
              main <- "DEPONS simulation output"
            } else {
              main <- x@title
            }
            # Make plot with either date or tick on x-axis
            the.blocknos <- sort(unique(x@dyn$block))
            if (!is.na(x@startday)) {
              if(!hasArg("xlim")) {
                xlim <- NULL
              } else {
                xlim <- list(...)[["xlim"]]
              }
              if(!hasArg("xlab")) {
                xlab <- "time"
              } else {
                xlab <- list(...)[["xlab"]]
              }
              plot(z$real.time, z$count,
                   xlab=xlab, ylab=ylab, main=main, col=col, type="n",
                   xlim=xlim, ylim=ylim, axes=axes, lwd=lwd, lty=lty)

              for (i in 1:length(the.blocknos)) {
                b <- the.blocknos[i]
                if(sum(z$count[z$block==b])==0) {
                  all.messages <- rbind(all.messages, paste("No obs in block", b))
                  names(all.messages) <- "messages"
                  next
                }
                if(type=="l") lines(z$real.time[z$block==b], z$count[z$block==b],
                                    col=col[i], lty=lty, lwd=lwd)
                if(type=="p") graphics::points(z$real.time[z$block==b], z$count[z$block==b],
                                     col=col[i], pch=pch, cex=cex)
              }
            } else {
              if(!hasArg("xlim")) {
                xlim <- c(min(x@dyn$tick), max(x@dyn$tick))
              } else {
                xlim <- list(...)[["xlim"]]
              }
              if(!hasArg("xlab")) {
                xlab <- "tick"
              } else {
                xlab <- list(...)[["xlab"]]
              }
              plot(z$tick, z$count,
                   xlab=xlab, ylab=ylab, main=main, col=col, type="n",
                   xlim=xlim, ylim=ylim, axes=axes, lwd=lwd, lty=lty)
              for (i in 1:length(the.blocknos)) {
                b <- the.blocknos[i]
                if(sum(z$count[z$block==b])==0) {
                  all.messages <- rbind(all.messages, paste("No obs in block", b))
                  names(all.messages) <- "messages"
                  next
                }
                if(type=="l") lines(z$tick[z$block==b], z$count[z$block==b],
                                    col=col[i], lty=lty, lwd=lwd)
                if(type=="p") points(z$tick[z$block==b], z$count[z$block==b],
                                     col=col[i], pch=pch, cex=cex)
              }
            }
            invisible(all.messages)
          }
)