R/swemod.R

In swemod: Snow Water Equivalent Modeling with Parametric and Statistical Models

swemod <- function(data, model = c("deltasnow","jonas","pistocchi","sturm"), alt, region, snowclass, verbose = FALSE) {
  # wrapper for computation of SWE with different methods
  # 
  #-----------------------------------------------------------------------
  # some checks
  # data=data.frame(date=c("2010-01-01","2010-01-02","2010-01-03","2010-01-04","2010-10-01","2010-11-15","2011-02-01","2011-06-30","2012-07-01"),hs=c(100,110,115,111,65,70,20,5,2),stringsAsFactors=FALSE)
  # # check class of datum: must be character or POSIXlt with format YYYY-MM-DD
  # if (inherits(datum,"character")){
  #   datum <- as.POSIXlt(datum,format="%Y-%m-%d")
  # } else if (!inherits(datum,"POSIXlt")){
  #   stop("deltasnow: datum must be either of class character or POSIXlt")
  # } 
  
  #data <- data.frame(date=as.character(seq(as.Date("2010-01-01"), as.Date("2010-01-22"),"day")),hs=c(0,0,0,1,2,3,4,0.9,0.7,0.5,0.2,0,0,0.9,0.5,1,3,2,0,0,0,0))
  
  #globalVariables(c("x.values", "y.values"))
  
  # if(!inherits(data,"data.frame"))
  #   stop("swemod: data must be given as data.frame")
  # 
  # if(!any((is.element(colnames(data), c("hs","date")))))
  #   stop("swemod: data must contain at least two columns named 'hs' and 'date'")
  
  # z <- zoo(data$hs,as.Date(data$date))
  # #if(!is.regular(z, strict = TRUE))
  # #  stop("swemod: date column must be strictly regular")
  # if(!all(z[,1] >= 0))
  #   stop("swemod: snow depth data must not be negative")
  # if(!all(!is.na(z)))
  #   stop("swemod: snow depth data must not be NA")
 
  # if(any(is.na(data)))
  #   stop("swemod: now depth data must not be NA")
  # if(!all(data >= 0))
  #   stop("swemod: snow depth data must not be negative")
  # if(!is.numeric(data))
  #   stop("swemod: snow depth data must be numeric")
  #   
  model <- match.arg(model, several.ok = TRUE)
  if(length(model) == 0) 
    model <- "deltasnow"

  
  #-----------------------------------------------------------------------
  # split into different models
  swe <- list()
  for(m in model){
   if(m == "deltasnow"){
     swe[["swe"]][[m]] <- swe.deltasnow(data, rho.max=440, rho.null=107, c.ov=0.004838623, k.ov=0.1996423, k.exp=0.02634508, tau=0.02, timestep = 24, verbose)
   } else if (m == "jonas"){
     swe[["swe"]][[m]] <- swe.jonas(data,alt,region)
   } else if (m == "pistocchi"){
     swe[["swe"]][[m]] <- swe.pistocchi(data)
   } else if (m == "sturm"){
     swe[["swe"]][[m]] <- swe.sturm(data,snowclass)
   }
  }
  
  swe[["date"]] <- data$date
  swe[["hs"]] <- data$hs
  
  class(swe) <- "swemod"
  return(swe)
  
}


# S3 function summary
summary.swemod <- function(object, ...){
  
  if(class(object) != "swemod")
    stop("swemod: Object must be of class 'swemod'.")
  
  models <- names(object$swe)
  if(length(models) == 0)
    stop("swemod: Cannot plot. No model was computed.")
  
  res <- c()
  for(m in models){
    res <- rbind(res,object$swe[[m]])  
  }
  rownames(res) <- models
  print(t(apply(res,1,summary)))
  
}




# S3 function plot
plot.swemod <- function(x, title = NULL, ...){
  
  if(class(x) != "swemod")
    stop("swemod: Object must be of class 'swemod'.")
  
  models <- names(x$swe)
  if(length(models) == 0)
    stop("swemod: Cannot plot. No model was computed.")
  
  colors <- c("red","green","blue","orange")
  
  # define maximum swe for plot outline
  ymax <- c()
  for(m in models){
    ymax <- c(ymax,max(na.omit(c(x$swe[[m]]))))  
  }
  ymax <- max(ymax)
  plot(as.Date(x$date,),x$swe[[1]],type="n",xlab="",ylab="HS [cm] / SWE [kg/m2]",ylim=c(0,ymax*1.2))
  axis.Date(1, at = seq(as.Date(x$date[1]), as.Date(x$date[length(x$date)]), by = "month"), format = "%m")
  n <- 1
  for(m in models){
    lines(as.Date(x$date),x$swe[[m]],type="l",col=colors[n])
    n <- n + 1
  }
  lines(as.Date(x$date),x$hs*100,type="l",lty=2,col="black")
  t <- ifelse(is.null(title),"SWE",title) #paste0("Chartreuse (",alts[s],"m)")
  legend(title=t,"topleft", 
         legend=c(models,"HS"),
         col=c(colors[1:length(models)],"black"), 
         lty=c(rep(1,length(models)),2), cex=0.8, bg="transparent", bty = "n")
  
  invisible(x)
}