R/plot.phenologydata.R

In phenology: Tools to Manage a Parametric Function that Describes Phenology and More

#' plot.phenologydata plots the observation data of phenology.
#' @title Plot the phenology from a add_phenology().
#' @author Marc Girondot \email{marc.girondot@@gmail.com}
#' @return Nothing.
#' @param x A result file generated by add_phenology
#' @param ... Parameters used by plot
#' @param series Name or number of series to be plotted or 'all'
#' @param season Which season to plot
#' @param moon If TRUE, the moon phase is ploted. Default is FALSE
#' @param col.observations Color of the points
#' @param col.grouped.observations Color of the lines indicating grouped observations
#' @param col.minimum.observations Color of the points indicating minimum counts
#' @description The function plot.phenologydata plots the observations.\cr
#' @family Phenology model
#' @examples
#' \dontrun{
#' library(phenology)
#' # Read a file with data
#' data(Gratiot)
#' # Generate a formatted list nammed data_Gratiot 
#' data_Gratiot <- add_phenology(Gratiot, name = "Complete", 
#' 		reference = as.Date("2001-01-01"), format="%d/%m/%Y")
#' 
#' plot(data_Gratiot)
#' }
#' @method plot phenologydata
#' @export

#plot.phenology <- function(x, ...) {

plot.phenologydata <- 
  function(x, ..., 
           series="all"                                                                                , 
           moon=FALSE                                                                                  , 
           season = NULL                                                                               , 
           col.observations = "black"                                                                  , 
           col.minimum.observations = "blue"                                                           ,
           col.grouped.observations = "green"                                                          ) {
    
    # x=NULL; series="all"; moon=FALSE; level=0.95; replicate.CI=1000; progressbar=TRUE; growlnotify=TRUE; show.plot=TRUE; resultmcmc = NULL; chain = 1; replicate.CI.mcmc = "all"; plot.objects = c("observations", "ML", "ML.SD", "MCMC.SD", "ML.quantiles", "MCMC.quantiles"); col.ML="black"; col.SD="red"; col.MCMC.quantiles="purple"; col.ML.quantiles="black"; col.observations = "black"; col.grouped.observations = "green"; col.observations = "black"; col.minimum.observations = "blue"; col.SD.polygon=rgb(red = 1, green = 0, blue = 0, alpha = 0.2); col.MCMC.quantiles.polygon=rgb(red = 160/255, green = 32/255, blue = 240/255, alpha = 0.2); col.ML.quantiles.polygon=rgb(red = 0, green = 0, blue = 0, alpha = 0.2) 
    

    p3p <- list(...)
    
    # result <- x
    if (series == "all") {
    data <- x
    } else {
    data <- x[series]
    }
        
    # kseries <- 1
    for(nmser in names(x)) {
      
      reference <- attributes(data[[nmser]])$reference
      if (is.null(reference)) {
        reference <- data[[nmser]][1, "Date"]
        referen_end <- data[[nmser]][nrow(data[[nmser]]), "Date"]
        premiermois <- as.POSIXlt(reference)$mon+1
        derniermois <- as.POSIXlt(referen_end)$mon+1
        refencours <- as.character(reference)
        substr(refencours, 9, 10) <- "01"
        
        # 1 8 <- 1
        # 8 1 <- 7
        # 1 1 <- 1
        # 8 9 <- 7
        
        if (premiermois < 7) {
          substr(refencours, 6, 7) <- "01"
          reference <- as.Date(refencours)
        } else {
          substr(refencours, 6, 7) <- "07"
          reference <- as.Date(refencours)
        }
      }
      
      nday <- 366 * (max(unlist(lapply(data, FUN = function(x) max(c(x[, "ordinal"], x[, "ordinal2"]), na.rm = TRUE))) %/% 366) + 1)
      
      
      # nday <- ifelse(as.POSIXlt(reference+365)$mday==as.POSIXlt(reference)$mday, 365, 366)
      
      vmaxx <- c(reference, reference+nday)
      
      vmaxy <- c(0, 0.1)
      # if (any((is.na(data[[nmser]]$ordinal2)))) {
      #   vmaxy[2] <- max(data[[nmser]]$nombre[(is.na(data[[nmser]]$ordinal2)) & 
      #                                          (!is.na(data[[nmser]]$nombre))])
      # }
      
      vmaxy[2] <- max(vmaxy[2], data[[nmser]]$nombre/ifelse(is.na(data[[nmser]]$ordinal2), 1, 
                                                            data[[nmser]]$ordinal2-data[[nmser]]$ordinal+1), na.rm = TRUE)
      # Si j'ai ordinal2, je dois prend n/(ordinal2-ordinal+1)
      
      if (vmaxy[2] ==0) vmaxy[2] <- 0.01
      
      x <- seq(from=reference, to=reference+nday-1, by="1 day")
      
      
      # je prépare une base
      # par(new=FALSE)

      pnp <- modifyList(list(xlab="Months", ylab="Counts", main=nmser, 
                             pch=16, cex=0.5, xlim=vmaxx, ylim=vmaxy, type="n", bty="n"), p3p)
      do.call("plot", modifyList(pnp, list(x=x, y=rep(0, length(x)))))
      
      if (moon) {
        par(xpd=TRUE)
        moony <- ScalePreviousPlot()$ylim["begin"] + (ScalePreviousPlot()$ylim["end"] - ScalePreviousPlot()$ylim["begin"]) * 1.06
        mp<-moon.info(x, phase=TRUE)
        mpT1<-ifelse((mp!="FM") | (is.na(mp)), FALSE, TRUE)
        mpT2<-ifelse((mp!="NM") | (is.na(mp)), FALSE, TRUE)
        
        xnewmoon <- ifelse(x[mpT1]>=ScalePreviousPlot()$xlim["begin"] & x[mpT1]<=ScalePreviousPlot()$xlim["end"], TRUE, FALSE)
        xfullmoon <- ifelse(x[mpT2]>=ScalePreviousPlot()$xlim["begin"] & x[mpT2]<=ScalePreviousPlot()$xlim["end"], TRUE, FALSE)
        points(x[mpT1][xnewmoon], rep(moony, length(x[mpT1]))[xnewmoon], cex=1, bg="black", col="black", pch=21, xpd=TRUE)
        points(x[mpT2][xfullmoon], rep(moony, length(x[mpT2]))[xfullmoon], cex=1, bg="white", col="black", pch=21, xpd=TRUE)
        par(xpd=FALSE)
      }
              
        col_ec <- ifelse(data[[nmser]]$CountTypes[is.na(data[[nmser]]$Date2)] == "exact", col.observations, col.minimum.observations)
        points(x = data[[nmser]]$Date[is.na(data[[nmser]]$Date2)], 
               y=data[[nmser]]$nombre[is.na(data[[nmser]]$Date2)], 
               pch=16, 
               col=col_ec, cex=0.5)
        
        for(i in 1:dim(data[[nmser]])[1]) {
          
          if (!is.na(data[[nmser]]$ordinal2[i])) {
            x0<-data[[nmser]]$Date[i]
            x1<-data[[nmser]]$Date2[i]
            lgt01<-as.numeric(data[[nmser]]$Date2[i]-data[[nmser]]$Date[i]+1)
            y0<-data[[nmser]]$nombre[i]/lgt01
            y1<-y0
            segments(x0, y0, x1=x1, y1=y1, col=col.grouped.observations, lwd=2)
          }
        }
      }
    
    return(invisible())
    
  }