Several major modifications. Everything seems to be working fine now.

NAMESPACE

@@ -2,7 +2,7 @@
 
 export(dgvm3d.options)
 export(drawEllipsoid)
-export(establishPatch)
+export(establishVegetation)
 export(gapless.rank)
 export(getCone)
 export(getHexagon)

R/classes.R

@@ -115,25 +115,19 @@ dgvm3d.options <- function(x=NULL,
   dgvm3d.options("default")
 }
 
-#' Year 2000 snapshot of a LPJ-GUESS simulation
-#'
-#'  Simulation with 12 patches at 13 locations (default test gridlist in LPJ-GUESS)
-#'
-#'  @name dgvm3d.snapshots
-#'  @docType data
-#'  @author Joerg Steinkamp \email{steinkamp.joerg@@gmail.com}
-#'  @keywords data
-"dgvm3d.snapshots"
-
 #' LPJ-GUESS gridlist of 13 locations
 #'
+#' This is the default test gridlist shipped with the LPJ-GUESS source code.
+#'
 #'  @name dgvm3d.locations
 #'  @docType data
 #'  @author Joerg Steinkamp \email{steinkamp.joerg@@gmail.com}
 #'  @keywords data
 "dgvm3d.locations"
 
-#' timeseries data from 1865-2005 is 5 year steps
+#' timeseries data from 1865-2005 in 5 year steps
+#'
+#' A list of 13 data.frames with simulation results of a model run without random patch disturbance at the 13 locations defined in \code{\link{dgvm3d.locations}}.
 #'
 #'  @name dgvm3d.succession
 #'  @docType data

R/geometry.R

@@ -170,7 +170,7 @@ getCone <- function(radius=0.5, height=1, faces=72, close=FALSE) {
 
 #' draw a quadrilateral ellipsoid
 #'
-#' I copied and modified it from 'Code demonstations' rgl::shapes3d (search with "??ellipsoid" on commandline)
+#' I copied and modified it from 'Code demonstations' rgl::shapes3d (searched with "??ellipsoid" on commandline)
 #'
 #' @param rx radius in x direction
 #' @param ry radius in y direction

R/input.R

@@ -13,6 +13,15 @@
 #' @importFrom utils read.table
 #' @export
 #' @author Joerg Steinkamp \email{steinkamp.joerg@@gmail.com}
+#' @examples
+#' \dontrun{
+#' data(dgvm3d.locations)
+#' dgvm3d.succession = list()
+#' for (i in 1:nrow(dgvm3d.locations)) {
+#'   dgvm3d.succession[[i]] = read.LPJ("~/WIP/Establishment/output/disturb/vegstruct.out", lon=dgvm3d.locations$Lon[i], lat=dgvm3d.locations$Lat[i])
+#'   dgvm3d.succession[[i]] = dgvm3d.succession[[i]][!(dgvm3d.succession[[i]]$Year %% 5) & dgvm3d.succession[[i]]$Year > 1859, ]
+#' }
+#' }
 read.LPJ <- function(file="vegstruct.out", stand.id=1, patch.id=NULL, year=NULL, lon=NULL, lat=NULL, grass=FALSE) {
   ## location.names <- read.table("/Users/jsteinkamp/WIP/Establishment/output/gridlist.txt", sep="\t", col.names=c("Lon", "Lat", "Name"))
   ## file = "/Users/jsteinkamp/WIP/Establishment/output/vegstruct.out"

R/stand.R

@@ -4,6 +4,7 @@
 #' If soil is a matrix, the number of columns must be equal to npatch. In that way each patch can have its own soil depth.
 #' The patches represented as hexagons can either be arranged in a square or in a line. The later one for example to represent a time series (succession).
 #' @param npatch number of patches
+#' @param year the initialization year
 #' @param soil a vector or matrix of soil depths.
 #' @param z the height of each patch.
 #' @param layout patch layout ('square' or 'linear'), a two element vector with number of rows/colums. A matrix for layout (not yet ready).
@@ -36,10 +37,12 @@ initStand <- function(npatch=1, year=2000, soil=c(0, -0.5, -1.5), z=0, layout="s
   }
 
   hexagon <- getHexagon(area=dgvm3d.options("patch.area"), z=c(0, -1))
+
   if (typeof(layout) == "character") {
     if (layout=="square") {
       nxy.max = ceiling(sqrt(npatch))
       nxy.min = floor(sqrt(npatch))
+      ## with 3 patches the above got the values 2 and 1, so 3 was never reached in the loop below
       if (npatch==3)
         nxy.min=2
       layout = c(nxy.min, nxy.max)
@@ -89,6 +92,7 @@ initStand <- function(npatch=1, year=2000, soil=c(0, -0.5, -1.5), z=0, layout="s
       }
     }
   }
+
   return(new("Stand", area=dgvm3d.options("patch.area"), year=year, hexagon=hexagon, layout=layout, composition=composition, patch.pos=t(patch.pos), patches=patches))
 }
 
@@ -126,6 +130,7 @@ updateStand <- function(stand, vegetation, year=NULL) {
       new.patch.veg$dnn = NA
       old.vid = unique(stand@patches[[i]]@vegetation$VID)
       if (length(old.vid) > 0) {
+        ## removal of killed individuals
         for (j in 1:length(old.vid)) {
           remain = sum(new.patch.veg$VID==old.vid[j])
           old.trees = stand@patches[[i]]@vegetation[stand@patches[[i]]@vegetation$VID==old.vid[j], ]
@@ -134,7 +139,7 @@ updateStand <- function(stand, vegetation, year=NULL) {
           new.patch.veg[new.patch.veg$VID==old.vid[j], "y"] = old.trees$y[1:remain]
         }
       }
-      stand@patches[[i]]@vegetation = establishPatch(new.patch.veg, stand@hexagon@supp$inner.radius)
+      stand@patches[[i]]@vegetation = establishVegetation(new.patch.veg, stand@hexagon@supp$inner.radius)
     } else {
       stand@patches[[i]]@vegetation = data.frame()
     }

R/vegetation.R

@@ -17,8 +17,8 @@
 #' veg$Crownarea = veg$DBH * 10
 #' veg$LeafType  = sample(0:1, nrow(veg), replace=TRUE)
 #' veg$ShadeType = sample(0:1, nrow(veg), replace=TRUE)
-#' stand@patches[[1]]@vegetation = establishPatch(veg, stand@hexagon@supp[['inner.radius']])
-establishPatch <- function(vegetation=NULL, radius=1) {
+#' stand@patches[[1]]@vegetation = establishVegetation(veg, stand@hexagon@supp[['inner.radius']])
+establishVegetation <- function(vegetation=NULL, radius=1) {
   if (is.null(vegetation))
     stop("'vegetation' data.frame is missing!")
 
@@ -29,7 +29,7 @@ establishPatch <- function(vegetation=NULL, radius=1) {
   est.beta.param   <- dgvm3d.options("establish.beta.parameters")
 
   if (dgvm3d.options("verbose")) {
-    message("### establishPatch ###")
+    message("### establishVegetation ###")
     message(sprintf("Using %i samples in max. %i repetitions (max. crown radius overlap: %0.3f).", samples[1], samples[2], overlap))
     message(paste0("Sorting by '", sort.column[1], "' in '", sort.column[2], "' order."))
   }
@@ -40,11 +40,16 @@ establishPatch <- function(vegetation=NULL, radius=1) {
     } else {
       vegetation = eval(parse(text=paste0("vegetation[with(vegetation, order(-", sort.column[1],")), ]")))
     }
+  } else {
+    warning(paste0("Column '", sort.column[1], "' does not exist. No sorting perfomed."))
+    message(paste0("Column '", sort.column[1], "' does not exist. No sorting perfomed."))
   }
 
+  ## if vegetation table is empty, e.g. after disturbance
   if (nrow(vegetation)==0)
     return(vegetation)
 
+  ## check for present position columns
   if (!all(c("x","y") %in% colnames(vegetation))) {
     if (dgvm3d.options("verbose"))
       message("New establishment.")
@@ -54,12 +59,15 @@ establishPatch <- function(vegetation=NULL, radius=1) {
     message("Establishing new trees.")
   }
 
+  ## first tree
   if (all(is.na(vegetation$x)) || all(is.na(vegetation$x))) {
     phi <- runif(1) * 2 * pi
     r   <- runif(1) * radius
     vegetation$x[1] = sin(phi) * r
     vegetation$y[1] = cos(phi) * r
   }
+
+  ## any other tree
   for (i in which(is.na(vegetation$x) | is.na(vegetation$y))) {
     trees.with.xy = which(!is.na(vegetation$x) & !is.na(vegetation$y))
     nwhile = 0
@@ -74,6 +82,7 @@ establishPatch <- function(vegetation=NULL, radius=1) {
       new.x <- sin(phi) * r
       new.y <- cos(phi) * r
 
+      ## distance to all other trees
       dist <- matrix(NA, length(trees.with.xy), samples[1])
       for (j in 1:samples[1])
         dist[,j] <- sqrt((new.x[j] - vegetation$x[trees.with.xy])^2 + (new.y[j] -vegetation$y[trees.with.xy])^2)
@@ -84,7 +93,8 @@ establishPatch <- function(vegetation=NULL, radius=1) {
 
       dist[dist < min.dist] = NA
 
-      ##dist <- colSums(dist)
+      ## choose the nearest tree for each location and apply the desired 'method'
+      ## on those
       dist = apply(dist, 2, min)
       if (!all(is.na(dist))) {
         if (establish.method=="max") {
@@ -127,7 +137,7 @@ establishPatch <- function(vegetation=NULL, radius=1) {
 #' veg$Crownarea = veg$DBH * 5
 #' veg$LeafType  = sample(1:2, nrow(veg), replace=TRUE)
 #' veg$ShadeType = sample(1:2, nrow(veg), replace=TRUE)
-#' stand@patches[[1]]@vegetation = establishPatch(veg, stand@hexagon@supp[['inner.radius']])
+#' stand@patches[[1]]@vegetation = establishVegetation(veg, stand@hexagon@supp[['inner.radius']])
 #' dummy = plant3D(stand, 1)
 #'
 #' stand3D(stand, 2)
@@ -136,7 +146,7 @@ establishPatch <- function(vegetation=NULL, radius=1) {
 #' veg$Crownarea = veg$DBH * 5 * rnorm(nrow(veg), 1, 0.1)
 #' veg$LeafType  = sample(1:2, nrow(veg), replace=TRUE)
 #' veg$ShadeType = sample(1:2, nrow(veg), replace=TRUE)
-#' stand@patches[[2]]@vegetation = establishPatch(veg, stand@hexagon@supp[['inner.radius']])
+#' stand@patches[[2]]@vegetation = establishVegetation(veg, stand@hexagon@supp[['inner.radius']])
 #' dummy = plant3D(stand, 2)
 plant3D <- function(stand=NULL, patch.id=NULL, crown.opacity=1) {
   if (is.null(patch.id))
@@ -152,8 +162,10 @@ plant3D <- function(stand=NULL, patch.id=NULL, crown.opacity=1) {
   color.column = dgvm3d.options("color.column")
 
   for (i in patch.id) {
-    if (nrow(stand@patches[[i]]@vegetation)>0) {
+    ## only in those patches with vegetation
+    if (nrow(stand@patches[[i]]@vegetation) > 0) {
       offset = stand@patch.pos[i, ]
+      ## chose the canopy color
       if (!any(names(stand@patches[[i]]@color.table) == "crown")) {
         n = eval(parse(text=paste0("length(unique(stand@patches[[i]]@vegetation$", color.column, "))")))
         stand@patches[[i]]@color.table[['crown']] = crown.colors(n)
@@ -163,10 +175,8 @@ plant3D <- function(stand=NULL, patch.id=NULL, crown.opacity=1) {
         col = eval(parse(text=paste0("crown.colors(length(unique(stand@patches[[i]]@vegetation$", color.column,")))")))
         stand@patches[[i]]@color.table[['crown']] = col
       }
-      if (nrow(stand@patches[[i]]@vegetation) > 0) {
-        for (j in 1:nrow(stand@patches[[i]]@vegetation)) {
-          tree3D(stand@patches[[i]]@vegetation[j, ], offset, col, opacity=crown.opacity)
-        }
+      for (j in 1:nrow(stand@patches[[i]]@vegetation)) {
+        tree3D(stand@patches[[i]]@vegetation[j, ], offset, col, opacity=crown.opacity)
       }
     }
   }

R/visualize.R

@@ -8,7 +8,7 @@
 #' @return a stand object
 #' @export
 #' @examples
-#' stand=succession(dgvm3d.succession, init.year=1865, years=c(1865,seq(1875, 2000, 25)), patch.id=c(2,4,11))
+#' stand=succession(dgvm3d.succession[[3]], init.year=1865, years=c(1865, seq(1875, 2000, 25)), patch.id=sample(1:12, 3))
 #' stand3D(stand)
 #' stand=plant3D(stand)
 succession <- function(vegetation, stand.id=1, patch.id=NULL, init.year=1901, years=seq(1950, 2000, 10)) {
@@ -20,7 +20,7 @@ succession <- function(vegetation, stand.id=1, patch.id=NULL, init.year=1901, ye
   stand <- initStand(npatch=npatch, year=init.year)
 
   for ( i in 1:npatch) {
-    stand@patches[[i]]@vegetation = establishPatch(subset(init.vegetation, PID==patch.id[i]), stand@hexagon@supp[['inner.radius']])
+    stand@patches[[i]]@vegetation = establishVegetation(subset(init.vegetation, PID==patch.id[i]), stand@hexagon@supp[['inner.radius']])
     stand@patches[[i]]@pid = patch.id[i]
   }
 
@@ -60,7 +60,7 @@ succession <- function(vegetation, stand.id=1, patch.id=NULL, init.year=1901, ye
 #' @export
 #' @author Joerg Steinkamp \email{steinkamp.joerg@@gmail.com}
 #' @examples
-#' stand=snapshot(dgvm3d.succession)
+#' stand=snapshot(dgvm3d.succession[[1]])
 snapshot <- function(vegetation, stand.id=1, patch.id=NULL, year=2000) {
   SID=PID=Year=NULL
   patch.ids = unique(vegetation$PID)
@@ -77,13 +77,13 @@ snapshot <- function(vegetation, stand.id=1, patch.id=NULL, year=2000) {
     if (!all( c("x", "y") %in% colnames(vegetation))) {
       if (dgvm3d.options("verbose"))
         message("Randomly distributing trees in patch")
-      stand@patches[[i]]@vegetation = establishPatch(subset(vegetation, PID==patch.ids[i]), stand@hexagon@supp[['inner.radius']])
+      stand@patches[[i]]@vegetation = establishVegetation(subset(vegetation, PID==patch.ids[i]), stand@hexagon@supp[['inner.radius']])
       stand@patches[[i]]@pid = patch.ids[i]
     } else {
       if (any(is.na(vegetation$x)) || any(is.na(vegetation$y))) {
         warning("NAs in tree position. Random redistribution of missing positions.")
         message("NAs in tree position. Random redistribution of missing positions.")
-        stand@patches[[i]]@vegetation = establishPatch(subset(vegetation, PID==patch.ids[i]), stand@hexagon@supp[['inner.radius']])
+        stand@patches[[i]]@vegetation = establishVegetation(subset(vegetation, PID==patch.ids[i]), stand@hexagon@supp[['inner.radius']])
       } else {
         if (dgvm3d.options("verbose"))
           message("Using valid x/y values present in 'vegetation data.frame.")