R/constructplant.R

In YplantQMC: Plant Architectural Analysis with Yplant and QuasiMC

#'Construct a 3D plant
#'
#'@description Read legacy-style Yplant input files into a special object, to be used in
#'any analysis in YplantQMC.
#'
#'The function constructs an object of class \code{plant3d}, based on Yplant
#'input files (.p and .l/.lf). Various methods exist for \code{plant3d}
#'objects, in particular \code{\link{plot.plant3d}} and
#'\code{\link{summary.plant3d}}.
#'
#'For batch analyses, the function \code{readplantlist} reads a number of files
#'at a time, and stores the results in a special list (of class
#'\code{plant3dlist}).
#'
#'Three plants are provided with \code{YplantQMC} (and automatically loaded)
#':\code{toona}, \code{pilularis} and \code{sugarmaple}. See
#'\code{\link{plantexamples}}.
#'
#'To learn about the format of P and L files, read the detailed account on the
#'Prometheus wiki (Pearcy, Falster & Duursma 2011): \url{http://goo.gl/Hmyv6}.
#'
#'For legacy Yplant users (Pearcy and Yang 1996, see
#'\url{http://goo.gl/Hmyv6}), you will find that \code{constructplant} is much
#'more robust with respect to malformed input files. It will also attempt to
#'write error messages when things go wrong.
#'
#'The \bold{Q file} format is an alternative to .P files, and is much easier to
#'use if the virtual plant does not have stem sections. There are seven
#'columns:
#'
#'\describe{ \item{X,Y,Z}{ Coordinates of the leaf base} \item{ang,az}{ Angle
#'and azimuth of the normal to the leaf surface} \item{or}{Orientation (azimuth
#'angle) of the midrib} \item{L}{Leaf length} }
#'
#'The file is space-delimited (such as the output of \code{write.table}), and
#'includes column headers (exactly named as above).
#'
#'The \bold{leafplantkey} file is a convenient way to organize a large number
#'of plant files. This is a simple comma-separated text file without headers.
#'The order is pfile,lfile (without quotes). For example, a "leafplantkey.txt"
#'file may look like this: \preformatted{ acaflo1.p,acaflo.l acaflo2.p,acaflo.l
#'acaflo3.p,acaflo.l acamyr1.p,acamyr.l acamyr2.p,acamyr.l acamyr3.p,acamyr.l
#'acasua1.p,acasua.l acasua2.p,acasua.l acasua4.p,acasua.l acasuaR02.p,acasua.l
#'acasuaR05.p,acasua.l acasuaR09.p,acasua.l }
#'
#'@param input One of several possible inputs that contain the plant structure.
#'@param lfile A leaf file. If not provided, a built-in triangle leaf is used.
#'@param multiplier Multiplies length dimensions, e.g. to change units.
#'@param X0,Y0,Z0 New x,y,z coordinate of the stem base.
#'@param warn If TRUE, writes warnings of minor issues with P file format.
#'@param quiet If TRUE, no messages are ever shown.
#'@param pfiles,lfiles Vectors of .p files and .l files.
#'@param lpk Optionally, a 'leafplantkey' file. See Details.
#'
#'@return In the case of \code{constructplant}, an object of class
#'\code{plant3d}. For \code{readplantlist}, an object of class
#'\code{plant3dlist} (which is simply a list of objects as generated by
#'\code{constructplant} to ease batch analyses).
#'@author Remko Duursma
#'@seealso \code{\link{plot.plant3d},\link{readp}}
#'@keywords misc
#'@examples
#'
#'
#'\dontrun{
#'# Read one plant:
#'myplant <- constructplant("sompfile.p","somelfile.l")
#'
#'# Pfile was in cm - should be in mm. Multiply all length dimensions by 10.
#'myplant <- constructplant("sompfile.p","somelfile.l", multiplier=10)
#'
#'# Read a couple of plants. 
#'myplants <- constructplant(pfiles=c("plant1.p","plant2.p"), lfiles=rep("leaf.l",2))
#'
#'}
#'
#'
#'@aliases constructplant readplantlist
#'@rdname constructplant
#'@export
constructplant <- function(input=NULL,  # either pfile,pdata, qfile,or qdata.
                           lfile=NULL,  
                           multiplier=1.0,
                           X0=0, Y0=0, Z0=0, 
                           warn=FALSE, quiet=FALSE){
  
  # Constants
  k <- pi/180
  
  # Local functions
  getxyz <- function(azim = 45, angle=30, len=3, origin=list(x=0,y=0,z=0)){
    
    azim <- azim * k
    angle <- angle * k
    
    z <- origin$z + len*sin(angle)  # z coordinates of line starting and endpoint
    d <- len*cos(angle)  			# distance of line from origin, as seen from above.
    x <- origin$x + d*sin(azim)     # x coordinates
    y <- origin$y + d*cos(azim)     # y coordinates
    
    return(list(x=x,y=y,z=z))
  }
  
  #--- Read plant input
  # Decide the inputformat (Q or P), and whether a filename or 
  # dataframe is given.
  if(is.character(input)){
    
    if(grepl("[.][qQ]$", input)){
      inputformat <- "Q"
      qfile <- input
      pfile <- NA
      pdata <- NA
      
      qdata <- read.table(qfile, header=TRUE)
      names(qdata) <- c("X","Y","Z","An.3","Az.3","Or","L.3")
      
      # For now, only one leaf type allowed...
      qdata$Lt <- 1
    }
    
    if(grepl("[.][pP]$", input)){
      inputformat <- "P"
      pfile <- input
      qfile <- NA
      pdata <- readp(pfile)
      qdata <- NULL
    }
  }
  if(is.data.frame(input)){
    
    qvars <- c("ang","az","len","or","x","y","z")  
    pvars <- c("mn","n","ste","lt","az")
    
    if(all(qvars %in% tolower(names(input)))){
      inputformat <- "Q"
      qdata <- input
      names(qdata) <- c("X","Y","Z","An.3","Az.3","Or","L.3")
      qdata$Lt <- 1
      pfile <- qfile <- "from dataframe"
    }
    if(all(pvars %in% tolower(names(input)))){
      inputformat <- "P"
      pdata <- input
      qdata <- NA
      pfile <- qfile <- "from dataframe"
    }
  }
  
  #--- Read leaf file
  # If file name given , read it:
  if(is.character(lfile))
    ldata <- readl(lfile)
  
  # If a leaf object is given, use that.
  if(inherits(lfile, "leaffile"))
    ldata <- lfile
  
  # Leaf data
  if(is.null(lfile)){
    if(.Platform$OS.type != "windows" || !interactive())
      stop("Please specify a leaf file.")
    
    # .... needs to be updated again (readl has changed).
    if(!quiet)message("Using built-in triangle leaf.")
    ldata <- structure(list(structure(list(XYZ = structure(c(0, 5, 0, -5, 
                                                             0, 0, 10, 0, 0, 0, 0, 0), 
                                                           .Dim = c(4L, 3L), 
                                                           .Dimnames = list(NULL, 
                                                                            c("X", "Y", "Z"))),
                                           midribpoints = c(1, 3), 
                                           leafpars = structure(c(-999, -999, -999, 
                                                                  -999, -999, NA), 
                                           .Names = c("Amax", "Rd","QY", "shape",
                                                     "absorp", "")), leaftype = 1L, 
                                           nleaftypes = 1L, midriblen = 10, 
                                           leafshape = 0.5), 
                                           .Names = c("XYZ", "midribpoints", "leafpars", 
                                                 "leaftype", "nleaftypes", "midriblen", 
                                                 "leafshape"), class = "leaffile")), 
                                            class = "leaffile")
    lfile <- "triangle"
  }
  
  # If still no leaf read, stop here:
  if(is.null(lfile)){
    stop("Need a leaf file to continue.\n")
  }
  
  # Normalize ldata, so that 'length' (length of the midrib!) = 1
  # ldata$XYZ <- ldata$XYZ / max(ldata$XYZ[,"Y"])
  ldata_scaled <- ldata
  for(i in 1:length(ldata)){
    ldata_scaled[[i]]$XYZ <- ldata[[i]]$XYZ / ldata[[i]]$midriblen
  }
  
  #--- P Format
  if(inputformat == "P"){
    
    pdata$Lt[pdata$L.3 == 0] <- 0
    
    # Convert units (if needed).
    if(multiplier != 1.0){
      pdata$L <- multiplier * pdata$L
      pdata$L.1 <- multiplier * pdata$L.1
      pdata$L.2 <- multiplier * pdata$L.2
      pdata$L.3 <- multiplier * pdata$L.3
      pdata$D <- multiplier * pdata$D
      pdata$D.1 <- multiplier * pdata$D.1
      pdata$D.2 <- multiplier * pdata$D.2
    }
    
    # Any L == 0? Cannot be - set to 0.1.
    pdata$L[pdata$L == 0] <- 0.1
    
    endcoordinates_stem <- vector("list", nrow(pdata))
    endcoordinates_branch <- vector("list", nrow(pdata))
    Petioles <- vector("list", nrow(pdata))
    Stems <- vector("list", nrow(pdata))
    Branches <- vector("list", nrow(pdata))
    
    for(i in 1:nrow(pdata)){
      
      # Note the +1 : in the p format, 0 is the first node.
      node <- pdata$N[i] + 1
      mothernode <- pdata$MN[i] + 1
      thismothernode <- which(pdata$N == mothernode-1) # actual storage nr. in list
      
      # If mothernode does not exist, find the closest one in rank that does exist:
      if(length(thismothernode) == 0){
        thismothernode <- which.min(abs(pdata$N - mothernode))
      }
      
      if(length(thismothernode) > 1){
        if(warn)warning("Duplicated node number. First one used.\n")
        thismothernode <- thismothernode[1]
      }
      
      ste <- pdata$ste[i]
      
      # Find coordinates of starting point for current segment.
      if(i == 1){
        origin <- list(x=X0,y=Y0,z=Z0)   # Use zero for first segment.
      } else {
        if(ste==1)origin <- endcoordinates_stem[[thismothernode]]
        if(ste==2)origin <- endcoordinates_branch[[thismothernode]]
        
        # For robustness: if a node is missing from the file, connect segment to
        # previous node. 
        if(is.null(origin)){
          if(ste==1)origin <- endcoordinates_stem[[thismothernode-1]]
          if(ste==2)origin <- endcoordinates_branch[[thismothernode-1]]
          if(warn)warning("Origin NULL at node", node, 
                          "- current node connected to previous node")
        }	
      }
      
      # Get end coordinates of a stem segment ('origin' is the start coordinate)
      endcoordinates_stem[[i]] <- getxyz(azim = pdata$Az[i], angle=pdata$An[i],
                                         len=pdata$L[i], origin)
      Stems[[i]] <- list()
      Stems[[i]]$xyz <- list(from=unname(unlist(origin)), 
                             to=unname(unlist(endcoordinates_stem[[i]])))
      Stems[[i]]$mothernode <- mothernode
      Stems[[i]]$node <- node
      Stems[[i]]$diam <- pdata$D[i]
      
      # Branch segment.
      endcoordinates_branch[[i]] <- getxyz(azim = pdata$Az.1[i], angle=pdata$An.1[i], 
                                           len=pdata$L.1[i], origin)
      Branches[[i]] <- list()
      Branches[[i]]$xyz <- list(from=unname(unlist(origin)), 
                                to=unname(unlist(endcoordinates_branch[[i]])))
      Branches[[i]]$mothernode <- mothernode
      Branches[[i]]$node <- node
      Branches[[i]]$diam <- pdata$D.1[i]
      
      # Get petiole coordinates, if leaf length of leaf greater than zero.
      # Note that angles of petiole and leaf are multiplied by -1, by definition.
      if(pdata$L.3[i] > 0){
        petiole_endpoint <- getxyz(azim = pdata$Az.2[i], angle=pdata$An.2[i],
                                   len=pdata$L.2[i], origin)
        Petioles[[i]] <- list()
        Petioles[[i]]$xyz <- list(from=unname(unlist(origin)), 
                                  to=unname(unlist(petiole_endpoint)))
        Petioles[[i]]$mothernode <- mothernode
        Petioles[[i]]$node <- node
      }
    }
  }
  
  
  
  # This list will store xyz coordinates of all leaf edges.
  newleaves <- list()
  
  # Data with only the leaf information.
  if(inputformat == "P"){
    whichleaves <- which(pdata$Lt >= 1)
    pdatL <- pdata[whichleaves,]
  } 
  if(inputformat == "Q"){
    pdatL <- qdata
    whichleaves <- 1:nrow(pdatL)
  }
  
  # Shorthand
  OR <- pdatL$Or      # orientation
  AZ <- pdatL$Az.3    # azimuth
  AN <- pdatL$An.3	  # angle
  AZ[AZ < 0] <- AZ[AZ < 0] + 360
  AN[AN==90] <- 89.999
  LEN <- pdatL$L.3    # length
  nleaves <- nrow(pdatL)
  
  # Get leaf base (=petiole end) and leaf tip coordinates	
  leafbasecoor <- list()
  leaftipcoor <- list()
  leafnodenumber <- list()
  ROLL <- INC <- c()
  
  #--- Calculate leaf edge coordinates for all leaves
  if(nleaves > 0){
    if(inputformat == "Q"){
      
      # Shift the plant base.
      qdata$X <- qdata$X + X0
      qdata$Y <- qdata$Y + Y0
      qdata$Z <- qdata$Z + Z0
      
      for(i in 1:nrow(qdata)){
        leafbasecoor[[i]] <- c(qdata$X[i],qdata$Y[i],qdata$Z[i])
      }
    }
    
    # Save the coordinates of the leaf base, and the node numbers.
    if(inputformat == "P"){
      for(i in 1:nleaves){
        leafbasecoor[[i]] <- Petioles[[whichleaves[i]]]$xyz$to
        leafnodenumber[[i]] <- Petioles[[whichleaves[i]]]$node
      }
    }
    
    # Calculate leaf edge coordinates, and tip coordinate. 
    # (based on 'madeleafdirection' and 'leaf.init' in 'Dpunit.pas' in YPLANT.)
    mds <- list()
    for(i in 1:nleaves){
  
      mds[[i]] <- madeleafdirection(k*OR[i], k*AZ[i], k*AN[i])
      leaftipcoor[[i]] <- mds[[i]]$ld * LEN[i] + leafbasecoor[[i]]
      
      ld <- mds[[i]]$ld
      wd <- mds[[i]]$wd
      
      # find leaftype (Only supported for P files at the moment).
      if(inputformat == "P")
        leaft <- pdatL$Lt[i]
      else
        leaft <- 1
      
      ldatcur <- ldata_scaled[[leaft]]
      
      leafxy <- ldatcur$XYZ[,1:2]
      npoints <- nrow(leafxy)	
      X <- Y <- Z <- c()
      
      for(j in 1:npoints){
        
        tw <- leafxy[,"X"][j]
        tl <- leafxy[,"Y"][j]
        
        X[j] <- (tw*wd[1] + tl*ld[1]) * LEN[i] + leafbasecoor[[i]][1]
        Y[j] <- (tw*wd[2] + tl*ld[2]) * LEN[i] + leafbasecoor[[i]][2]
        Z[j] <- (tw*wd[3] + tl*ld[3]) * LEN[i] + leafbasecoor[[i]][3]
        
      }
      newleaves[[i]] <- ldatcur
      newleaves[[i]]$XYZ <- cbind(X,Y,Z)
      newleaves[[i]]$leafnodenumber <- ifelse(inputformat == "P", leafnodenumber[[i]], NA)
      newleaves[[i]]$leafnormal <- getleafnormal(newleaves[[i]])
      
    }
  } #END if(nleaves > 0)
  
  
  #--- Output
  l <- list()
  l$leaves <- newleaves
  l$nleaves <- nleaves

  if(inputformat=="P")l$stems <- Stems else l$stems <- NA
  if(inputformat=="P")l$branches <- Branches else l$branches <- NA
  if(inputformat=="P")l$petioles <- Petioles else l$petioles <- NA
  if(inputformat=="P")l$pdata <- pdata else l$pdata <- NA
  l$ldata <- ldata
  
  l$inputformat <- inputformat
  l$pfile <- pfile
  l$lfile <- lfile
  l$qdata <- qdata
  l$qfile <- qfile
  
  l$leaftipcoor <- do.call("rbind",leaftipcoor)
  l$leafbasecoor <- do.call("rbind",leafbasecoor)
  
  l$leafdata <- leafdata(l)
  l$leafarea <- sum(l$leafdata$area)*10^-6  
  
  # Make QMC input files.
  m <- makeQMCinput(l, writefile=FALSE)
  l$qmcinput <- m$qmcinput
  l$qmcinputfile <- m$inputfile
  l$qmcoutputfile <- m$outputfile
  
  class(l) <- "plant3d"
  return(l)
}

  
#'@rdname constructplant
#'@return An object of class \code{plant3dlist}
#'@export
readplantlist <- function(pfiles=NA, lfiles=NA, lpk="leafplantkey.txt", multiplier=1.0){
	  
	  # Init for 'bad' plants that cannot be read
	  nbad <- 0
	  whichbad <- c()
	  
	  if(length(pfiles) != length(lfiles))
	    stop("Provide equal number of plant and leaf files.")
	  # Read from leafplantkey
	  if(all(is.na(pfiles)) && file.exists(lpk)){
	    lpk <- read.csv(lpk, header=FALSE, stringsAsFactors=FALSE, strip.white=TRUE)
	    pfiles <- lpk[,1]
	    lfiles <- lpk[,2]
	  } 
	  
	  # in case pfiles,lfiles were provided as factors.
	  pfiles <- as.character(pfiles)
	  lfiles <- as.character(lfiles)
	  
	  # Multiplier (to change units, or fix strange files)
	  if(length(multiplier) != length(pfiles))
	    multipliers <- rep(multiplier, length(pfiles))
	  else
	    multipliers <- multiplier
	  
	  plants <- list()
	  for(i in 1:length(pfiles)){
	    tm <- system.time(plants[[i]] <- try(constructplant(pfiles[i], lfiles[i], 
	                                                        multiplier=multipliers[i]), silent=TRUE))
	    if(inherits(plants[[i]], "try-error")){
	      warning("Error constructing from ",pfiles[i],";",lfiles[i]," - plant skipped.")
	      nbad <- nbad + 1
	      whichbad <- c(whichbad, i)
	      next
	    } else {
	      cat("Plant",i,"of",length(pfiles),"(", pfiles[i], ") constructed in", unname(tm[3]),"sec.\n")
	      flush.console()
	    }
	  }
	  
	  plants[whichbad] <- NULL
	  pfilesnotread <- pfiles[whichbad]
	  lfilesnotread <- lfiles[whichbad]
	  
	  if(length(whichbad) > 0)
	    message("Some plants could not be read. See attributes(MyPlants)$notread")
	  
	  # pfiles, lfiles that were actually read in will be stored as attributes
	  pfiles <- sapply(plants, "[[", "pfile")
	  lfiles <- sapply(plants, "[[", "lfile")
	  
	  attributes(plants) <- list(pfiles=pfiles, lfiles=lfiles, nplants=length(plants),
	                             notread=data.frame(pfile=pfilesnotread,lfile=lfilesnotread))
	  class(plants) <- "plant3dlist"
	  
	  return(plants)
	}