MAESPA

Documented in Plotstand

#' Plot the stand in 3D
#'
#' @description Reads the MAESTRA trees file, and plots the stand in 3D. Supports all
#' MAESTRA crown shapes except the box shape.  Looks for the 'trees.dat' file
#' in the current working directory, unless specified (see Examples).  The XY
#' coordinates *must be present* in the 'trees.dat' file. Users will typically
#' only use the 'Plotstand' function.
#'
#' Optionally reads the crown shape from the 'str.dat' file, and plots the
#' correct crown shape for each species in the stand by reading the
#' multi-species namelists in 'confile.dat' and 'trees.dat'.
#'
#' The target trees are colored red (unless specified otherwise, see Details),
#' if the 'itargets' is specified in the confile.
#'
#' Attempts to read indivradx, indivrady, indivhtcrown, indivdiam, and
#' indivhttrunk namelists from the 'trees.dat' file. If any of these fail, the
#' 'all' versions are tried ('allradx', etc.).  Although MAESTRA runs fine when
#' no XY coordinates are provided, this plot function crashes. A future
#' implementation will calculate XY coordinates in the same way as MAESTRA.
#'
#' If the 'strfiles' parameter is set in 'confile.dat' (one str.dat file for
#' each species in the stand), these files are opened and used to set the crown
#' shape by species. Alternatively, you may specify crownshape as a parameter,
#' and override reading of str.dat by setting readstrfiles=FALSE.
#'
#' The 'nz' and 'nalpha' arguments specify the 'smoothness' of the crowns:
#' higher values provide more detailed triangulation of the crowns, at the
#' expense of speed.
#'
#' @param treesfile By default, the 'trees.dat' file in the current dir.
#' @param strfile Not used, yet.
#' @param crownshape Character,
#' "cone","elipsoid","ellipsoid","halfellipsoid","paraboloid","cylinder", or
#' abbreviation.
#' @param addNarrow Logical. Add arrow pointing North?
#' @param xyaxes Logical. Add annotated X and Y axes?
#' @param labcex Relative size of X and Y axis labels.
#' @param axiscex Relative size of X and Y axis annotation.
#' @param verbose If TRUE, writes more info to the screen while plotting.
#' @param readstrfiles Read the 'str.dat' file(s) to find out crown shape?
#' @param targethighlight Plot the target trees in red?
#' @param idate If multiple dates are provided for tree size variables, which
#' one to display.
#' @param path The folder where the input files are stored.
#' @param \dots See Details for a list of additional arguments recognized by `Plotstand`.
#' @return An rgl device is opened.
#' @author Remko Duursma
#' 
#' @details 
#' 
#' For large stands, the plot takes quite a while to complete. This implementation is certainly not optimized for speed.
#' Also, minimize the rgl window to greatly speed up the plotting process.
#'
#' The `Plotstand` function accepts a number of additional arguments that are used by subsidiary functions, these are:
#' 
#' \enumerate{
  #' \item **CL**: Crown length (m).
  #' \item **CW**: Crown width (m).
  #' \item **crowncolor**: The color of the tree crowns. Default, obviously, `forestgreen`.
  #' \item **stemcolor**: The color of the tree stems. Default `brown`.
  #' \item **x0,y0,z0**: Coordinates of crown base when calculating 3D coordinates.
  #' \item **HCB**: Height of crown base (m).
  #' \item **X,Y**: X- and Y-coordinates of tree stem base (m).
  #' \item **dbh**: Stem diameter (m). Converted to m if appears to be in cm.
  #' \item **nz**: Number of z divisions (increase number to get smoother crowns).
  #' \item **nalpha**: Number of angular divisions (increase number to plot smoother crowns).
  #' \item **m**: 3xN matrix (x,y,z coordinates in rows). Optional.
#' }
#' @examples
#'
#'
#' \dontrun{
#'
#' # Plot the 'trees.dat' file in the current working directory:
#' Plotstand()
#'
#' # Open a dialog box to select a trees.dat file:
#' Plotstand(file.choose())
#'
#' # Save a snapshot to a .png file.
#' # Note: make sure to move the 3D plot into view
#' # (so that other windows are not blocking it!)
#' snapshot3d('myforest.png')
#'
#' # For publication-quality graphs:
#' Plotstand(nz=50, nalpha=50)
#'
#'
#' }
#'
#' @export
Plotstand <- function(treesfile="trees.dat",
          				    strfile="str.dat",
          					  crownshape=c("cone","ellipsoid","round","halfellipsoid","paraboloid","cylinder"),
          					  readstrfiles=TRUE,
          					  targethighlight=TRUE,
          					  addNarrow=TRUE,
                      xyaxes=TRUE,
                      labcex=1,
                      axiscex=1,
                      verbose=FALSE,
                      idate=1,
                      path="",
                      ...){

  o <- getwd()
  on.exit(setwd(o))
  if(path != "")setwd(path)

	notrees <- readPAR(treesfile, "notrees", "plot")

	crownshapes <- rep(NA,notrees)

	haveconfile <- file.exists("confile.dat")

	if(targethighlight & !haveconfile){
		warning("No confile.dat found - target trees not highlighted")
		targethighlight <- FALSE
		crowncolors <- rep("forestgreen",notrees)
	}

	if(targethighlight){
		crowncolors <- rep("forestgreen",notrees)
		itargets <- readPAR("confile.dat", "itargets", "treescon", fail=FALSE)
		if(all(is.na(itargets))){
			warning("itargets not read in confile.dat")
		} else crowncolors[itargets] <- "red"
	}

	if(!haveconfile){
		strfiles <- strfile
	} else {
		strfiles <- readPAR("confile.dat","strfiles","species",fail=FALSE)
		if(all(is.na(strfiles)))
			strfiles <- strfile
		else {
			strfiles <- gsub("'","",strfiles)
		}
		if(!all(file.exists(strfiles)))
			stop("Not all strfiles are in the current working directory.")
	}

	if(readstrfiles){
		species <- readPAR("trees.dat","ispecies","speclist",fail=FALSE)
		if(all(is.na(species)))
      species <- rep(1,notrees)

		for(i in 1:notrees){
			crownshapes[i] <- tolower(readPAR(strfiles[species[i]],"cshape","canopy"))
	  }
		crownshapes <- gsub("'","",crownshapes)

	} else crownshapes[] <- match.arg(crownshape)


  xycoor <- readPAR(treesfile, "xycoords", "xy", fail=FALSE)
  if(all(is.na(xycoor)))stop("XY coordinates must be present in trees.dat file")
  xycoor <- matrix(xycoor,ncol=2,byrow=TRUE)
  X <- xycoor[,1]
  Y <- xycoor[,2]

  Bearing <- readPAR(treesfile, "bearing", "plot")

  # varname w/o 'indiv' or 'all'!
  readVar <- function(varname, idate=1){

    indvar <- paste0("indiv",varname)
    allvar <- paste0("all",varname)

    vals <- readPAR(treesfile, "values", indvar,fail=FALSE)
    whichvar <- indvar

    if(all(is.na(vals))){
      vals <- rep(readPAR(treesfile, "values", allvar), notrees)
      whichvar <- allvar
    }

    # dates?
    ndates <- readPAR(treesfile, "nodates", whichvar)
    if(all(is.na(ndates)) || ndates==1)return(vals)

    # else...
    vals <- matrix(vals, ncol=ndates, byrow=TRUE)
    return(vals[,idate])
  }

  radx <- readVar("radx", idate)
  CW <- 2*radx

  CL <- readVar("htcrown", idate)
	DBH <- readVar("diam", idate)
	if(max(DBH) > 3)DBH <- 0.01 * DBH

  # Avoid small DBH, rgl does not like that.
  DBH <- pmin(0.05, DBH)

  HCB <- readVar("httrunk", idate)

  if(any(is.na(c(CL,CW,X,Y,HCB,DBH))))stop("Missing values somewhere!")

  Openstand(treesfile)

  for(i in 1:notrees){
	if(verbose)message("Plotting tree number : ", i)
      plottree(crownshape=crownshapes[i], CL=CL[i], CW=CW[i],
          HCB=HCB[i], X=X[i], Y=Y[i], dbh=DBH[i], crowncolor=crowncolors[i])
  }

	if(addNarrow){
		X0 <- readPAR(treesfile,"x0","plot", fail=FALSE)
		if(is.na(X0))X0 <- 0
		Y0 <- readPAR(treesfile,"y0","plot",fail=FALSE)
		if(is.na(Y0))Y0 <- 0

		Xmax <- readPAR(treesfile,"xmax","plot")
		Ymax <- readPAR(treesfile,"ymax","plot")
		addarrow(x0=X0 + 0.1*Xmax,y0=Y0 + 0.1*Ymax,
			 len=0.1*(Ymax-X0),bearing=Bearing)
  }
	if(xyaxes){
	  rgl::par3d(cex=axiscex)
	  rgl::axes3d(c('x-','y-'))
	  rgl::par3d(cex=labcex)
	  rgl::title3d(xlab="X",ylab="Y")
	}

}