YplantQMC : Plant modelling in R with Yplant - QuasiMC Remko Duursma

Mik Cieslak

December 12, 2014

Contents 1

Introduction 1.1 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Installing YplantQMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Installing the QuasiMC executable . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2 Quick start guide 2.1 Step 1 : Reading input files . . . . . . . . . . . . . . . . . 2.2 Step 2 : Summarizing and plotting a plant . . . . . . . . . 2.3 Step 3 : Setting the geographical location . . . . . . . . . 2.4 Step 4 : Generating weather data . . . . . . . . . . . . . 2.5 Step 5 : Setting up the leaf photosynthesis model . . . . 2.6 Step 6 (optional) : Reading a hemispherical photograph 2.7 Step 7 : Running the daily model . . . . . . . . . . . . . . 2.8 Step 8 : Generating a report . . . . . . . . . . . . . . . .

2 2 2 3

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

3 3 4 4 5 5 5 6 6

3 Technical manual 3.1 About this section . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Overview of important functions . . . . . . . . . . . . . . . . . . . . 3.3 Input file types . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 constructplant : Reading input . . . . . . . . . . . . . . . . . . . . 3.5 randomplant : Generating plants with randomly distributed leaves 3.6 summary : Summarize the plant . . . . . . . . . . . . . . . . . . . . . 3.7 Plot the plant in 3D . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7.1 Side and top views . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Modifying virtual plants . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 STARbar : Calculating the projected and displayed leaf area . . . . 3.9.1 Projecting a plant . . . . . . . . . . . . . . . . . . . . . . . . 3.9.2 Displayed and projected leaf area . . . . . . . . . . . . . . . 3.10 setMet : Generating weather diurnal . . . . . . . . . . . . . . . . . 3.11 setPhy : Setting a leaf photosynthesis model . . . . . . . . . . . . . 3.11.1 Writing your own leaf gas exchange model . . . . . . . . . . 3.12 setHemi : Reading a hemispherical photograph . . . . . . . . . . . 3.13 YplantDay : Running a daily simulation . . . . . . . . . . . . . . . . 3.14 runYplant : A single YplantQMC run for custom simulations . . . . 3.15 QuasiMC options and settings . . . . . . . . . . . . . . . . . . . . . . 3.16 Batch processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.16.1 Reading and storing multiple plants . . . . . . . . . . . . . . 3.16.2 Plotting multiple plants . . . . . . . . . . . . . . . . . . . . . 3.16.3 Summarizing multiple plants . . . . . . . . . . . . . . . . . . 3.16.4 Calculating displayed area for multiple plants . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . . . . . . . . . .

6 6 6 7 8 8 9 11 12 12 13 13 14 14 14 15 16 17 20 21 21 21 22 22 22

1

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

. . . . . . . .

3.16.5 Running YplantDay for multiple plants . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

4 Unimplemented features of YplantQMC

23

5 Frequently Asked Questions (FAQ)

23

6 Acknowledgments

24

7 References

24

1

Introduction

1.1 Overview This document describes the YplantQMC package for R, a toolkit to analyse and visualize plant architecture. The YplantQMC model is a new implementation of the Yplant model developed by Pearcy and Yang (1996), in combination with the ray-tracer QuasiMC. A number of functions from the YpTools package have also been included (Duursma et al. 2012).

QuasiMC developed by Cieslak et al. (2008), which is part of the Virtual Laboratory plant modelling software created at the University of Calgary (http://www.algorithmicbotany.org). The YplantQMC model has the following capabilities:

• Reading Yplant input files to construct virtual plants in 3D. • Generating single plant canopies with random distribution of leaves. • Visualizing virtual plants in 3D. • Detailed plant architecture summary of virtual plants. • Ray-tracing : light absorption by individual leaves, based on transmission, reflection and scattering. • An extendable leaf physiology model, to estimate photosynthesis, transpiration and stomatal conductance for each leaf. • Diurnal simulations of radiation absorption and photosynthesis, based on standard simulated weather data, or input weather data. • A flexible toolkit for batch analyses, sensitivity analyses, and modifying virtual plants.

1.2

Installing YplantQMC

There are two parts to installing YplantQMC . First, a number of R packages need to be installed (including YplantQMC ). Then, the QuasiMC raytracer needs to be installed. The YplantQMC package is hosted on CRAN, so just do :

install.packages("YplantQMC") This should also install the required packages rgl, geometry, LeafAngle and devtools. The package is loaded the usual way.

library(YplantQMC)

2

## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##

Loading Loading Loading Loading Loading Loading

required required required required required required

package: package: package: package: package: package:

Attaching package:

rgl geometry magic abind LeafAngle gplots

’gplots’

The following object is masked from ’package:stats’: lowess

-- Plant modelling with Yplant - QuasiMC for R (version 0.6-3) Refer to the manual for detailed instructions, or see the example in ?YplantDay.

If you see error messages, this is because some of the above packages were not installed successfully.

1.2.1

Installing the QuasiMC executable

Simply execute the following command to download the QuasiMC files needed for YplantQMC . This creates a new directory on your C: drive (C:/QuasiMC). Don’t worry, QuasiMC is only 7MB.

installQuasiMC() If successful, you will see the message:

'QuasiMC.exe' and 'enviro.e' have been placed in the folder 'c:/QuasiMC'. !- Please do not rename or move! If you are connecting to the internet through a proxy (as at some universities), use the following command;

installQuasiMC(proxy=TRUE) Note that the QuasiMC executable is currently available for Windows and MacOS only. If you are on a different system, everything in YplantQMC works except runYplant and YplantDay (but those are the most important functions).

2 Quick start guide In this section, the functionality of YplantQMC is shown, in the most basic way. We will read one plant, and show basic analyses: visualizing, summarizing, running YplantQMC and plotting the simulation results.

2.1 Step 1 : Reading input files To read in a plant (.P) and leaf (.L or .LF) file, make sure you have set the working directory in R to the directory that contains your files, and use the command below. If you like to skip this step, no problem : sugarmaple is a built-in plant in YplantQMC .

3

sugarmaple <- constructplant("sugarmaple.p","sugarmaple.l") See section 3.4 on page 8 for more options, and section 3.3 on page 7 for more information on the format of the input files. The following three plants are shipped with YplantQMC (and auto-loaded) :

sugarmaple A shaded sugar maple with near-horizontal leaves (Acer saccharum) pilularis A blackbutt sapling, glasshouse-grown (Eucalyptus pilularis) toona An Australian red cedar seedling (Toona australis)

2.2 Step 2 : Summarizing and plotting a plant The next step in any analysis is to make sure the 3D plant looks like it should, has the correct dimensions, leaf area, and so on. To make a 3D plot,

plot(sugarmaple)

To summarize a plant,

summary(sugarmaple) Try also the pilularis and toona plants. See section 3.6 on page 9 for example output and further details. To calculate displayed and projected areas, and STAR, see section 3.9.

2.3 Step 3 : Setting the geographical location A geographical location is simply a latitude/longitude, used in calculating the solar path. Set a location like this,

4

sydney <- setLocation(lat=-33.8, long=151.2) See the help page (?setLocation) for more details. Note that the longitude is not strictly necessary. To double check the location, make a simple plot of the world (output not shown, try this yourself):

plot(sydney)

2.4 Step 4 : Generating weather data A simple built-in weather generator simulates the diurnal course of temperature, solar position (altitude and azimuth), photosynthetically active radiation (PAR), vapor pressure deficit (VPD), relative humidity (RH) and the fraction direct radiation (fbeam). An example with fairly minimum settings:

sunnyday <- setMet(sydney, month=12, day=12, nsteps=12, Tmin=12, Tmax=28, PARday=24) A simple plot can be generated like this (not shown).

plot(sunnyday) See section 3.10 on page 14 for more details and settings.

2.5 Step 5 : Setting up the leaf photosynthesis model The following example sets up a (non-rectangular) lightresponse curve to be used in estimating leaf photosynthesis from light absorption:

clidlrc <- setPhy("lightresponse", leafpars=list(Amax=8.5, Rd=0.7, phi=0.045, theta=0.95, reflec=0.1, transmit=0.05)) where:

Amax Light-saturated rate of photosynthesis (µmol m-2 s-1 ). Rd Dark respiration (µmol m-2 s-1 ) phi Quantum yield (mol mol-1 ) theta Shape of the light response curve (0-1). reflec Reflectance (default = 0.1) transmit Transmittance (default = 0.1) See section 3.11 on page 14 for more details and settings. Note that unlike the original Yplant model, YplantQMC does not read leaf physiology parameters from the leaf file. They are stored, but not used at the moment in any leaf gas exchange calculations.

2.6 Step 6 (optional) : Reading a hemispherical photograph YplantQMC does not actually (yet) read in the actual hemispherical photograph, but it does read in a table with gap fractions, averaged by azimuth and altitude bins. An example file is included with the package. To read a hemiphoto in that format, use this command: 5

smallgap <- setHemi("canview1.csv") If you don’t read a hemiphoto, YplantQMC will assume an open sky. Two built-in hemiphotos are available (as yphemi objects) in YplantQMC : smallgap and largegap. These are used in the examples below. See section 3.12 on page 16 for more details and explanations.

2.7 Step 7 : Running the daily model A diurnal simulation is performed with the YplantDay function, like so:

maplerun <- YplantDay(sugarmaple, met=sunnyday, phy=clidlrc, hemi=largegap) All information needed to run the simulation is contained in the four objects: plant, met, phy and hemi.

2.8 Step 8 : Generating a report Finally, a standard PDF report is generated, using any combination of YplantQMC objects:

ypreport(plant=sugarmaple, phy=clidlrc, met=sunnyday, hemi=largegap, ypsim=maplerun) See ?ypreport for more details.

3 Technical manual 3.1 About this section The built-in help pages in R provide much more detail than what is given in the following technical manual. The purpose here is to give an overview of the available tools, give a few examples and common plots of plants and simulation results.

3.2 Overview of important functions The table below summarizes the most important functions in YplantQMC . For each function, the most crucial arguments are given (for example, the constructplant needs a pfile and an lfile). However, each function has many more options, some of which are described in later sections in the manual, but the reader is encouraged to read the individual help files for all settings (in R, type ?constructplant to read its help page, for example).

6

Function

constructplant(pfile, lfile, ...)

randomplant(...)

setLocation(lat, long, ...)

setMet(location, ...)

setPhy(leafmodel, ...)

setHemi(canfile, ...)

YplantDay(plant,met,phy,hemi, ...)

ypreport(...)

3.3

Description Reads input files, constructs a 3D plant in the YplantQMC format. Generates a plant (without stems) with leaves randomly distributed within a geometric shape, including cones, cylinders, and a few others. Set the location for the simulation, used to calculate the solar path. A daily weather generator, with a few settings. Can also read user-input weather variables. Construct a leaf physiology object, that summarizes the dependency of CO2 assimilation rate on light and other environmental drivers. Read a hemispherical photograph, to take into account canopy shading on light availability Run the YplantQMC model for a whole day, using the four objects: plant, weather, physiology and hemi-photo. Generate a PDF report with various standard graphs and tables, using any combination of the YplantQMC objects.

Input file types

I assume here that you are familiar with the two Yplant input files: the ’leaf file’ (.L or .LF) which describes the leaf outline, and the ’plant file’ (.P) which describes the plant branching structure, leaf lengths, and angles. See Pearcy et al. (2011) (Prometheus Wiki : http://tinyurl.com/6r2wn2q) for an in-depth description of the plant and leaf files, and how to make your own. Warning: the current version does not read leaf physiology parameters from the leaf file! A number of example files *will be* available on www.remkoduursma.com/yplantqmc (or on request). As an alternative to the standard Yplant input files, YplantQMC can read Q files, which contain no information on the branches or stems. As such, it is a simple alternative to P files, and particularly useful for simulations. This is an example Q file, with only three leaves:

X Y Z ang az or L -340.9 103.3 156.0 18.6 -73.8 121.6 26.9 -342.7 105.6 156.7 56.5 -90.2 357.8 25.9 -335.6 81.1 167.9 66.0 63.2 33.8 41.5 where:

7

X,Y,Z The coordinates of the leaf base ang The leaf angle (0 = flat, 90 = vertical). az Azimuth angle (of the normal to the leaf surface) or Orientation : azimuth of the midrib L Leaf length (mm) All angles are in degrees.

3.4 constructplant : Reading input The first step in any analysis is to read in the leaf and plant files that together describe the plant architecture, and store the result in a special object. Let’s read in two files, constructing a Eucalyptus pilularis.

pilularis <- constructplant("epilularis.P","epilularis.L") The resulting object is of the class plant3d. The print method gives some basic information:

pilularis ## ## ## ## ## ## ## ## ##

Yplant - plant object (class 'plant3d') -----------------------------Constructed from P file format Input files used : test/epilularis.p test/epilularis.l

Plant has 284 leaves on 378 nodes. To generate a detailed summary of your plant, use the summary() function.

There are a few options to customize the plant construction. First, the centre of the plant can be moved to a new x,y,z coordinate (in case the digitized plant was not properly centred). For example,

pilularis2 <- constructplant("epilularis.P","epilularis.L",X0=100,Y0=200,Z0=0) Another option handles the case where the plant file is not in units of mm (the default), but in cm (as can be the case when the plant was manually digitized). In that case, we need to mulitply all length units by 10 :

pilularis2 <- constructplant("epilularis.P","epilularis.L",multiplier=10) Alternatively, to select the input files with a menu, use this command:

myplant <- constructplant() Finally, we can read in a .Q file; see section 3.3 on page 7. (Note: write your own .Q file first)

myplant <- constructplant(qfile="someplant.Q", lfile="leaf.l")

3.5 randomplant : Generating plants with randomly distributed leaves This function generates a plant with leaves randomly distributed within a crown, which has one of six shapes (see ?randomplant for details). The output of the function is an object of class plant3d, just like that from constructplant.

8

There are two basic usages : either specify the number of leaves and leaf length, or the total leaf area and the leaf length. In both cases, the final result will not have exactly the number of leaves or total leaf area as specified, but it will be very close. The reason for this discrepancy is that it is difficult to place leaves inside a convex hull with exactly the dimensions as specified. The below example makes a cone-shaped plant. Note that if the leaf file (here, ’clidemia.l’ is omitted, constructplant will use a built-in triangle-shaped leaf). Also note: randomplant needs to know which leaf file you are using when you specify the total leaf area, because it determines the leaf shape.

coneplantQ <- randomplant(radius=250, height=500, crownshape="CONE", leaflen=65, LA=0.5, crownbase=100, lfile="Clidemia.l") ## Aiming for leaf area = 0.5 m2 ## Actual leaf area = 0.3047 m2 ## Number of leaves = 193 coneplant <- constructplant(coneplantQ, lfile="Clidemia.l") plot(coneplant)

3.6 summary : Summarize the plant More detailed information is given in the output of the summary.plant3d function. These include a number of advanced metrics that are described in later sections (see Details section). Here, I will give a brief outline of each of the variables calculated. First, let’s print the summary:

summary(pilularis) ## ## ## ## ## ## ## ## ## ##

Total leaf area (m2) Mean leaf size (cm2) Number of leaves Mean leaf length (cm) Mean leaf angle (deg) Mean leaf angle weighted by leaf area (deg) Ellipsoidal leaf angle dist. par.

= = = = = = =

0.537 18.9 284 9.36 55.9 57.3 1.01

Crown volume (convex hull) (dm3) = 392 Crown surface area (convex hull) (m2) = 2.83

9

## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##

Crown density (AL/AC) (m2 m-2) Crown width (m) Crown length (m) Crown projected area (2D convex hull) (m2) Total height(m) Crown shape index (-)

= = = = = =

0.19 0.944 1.14 0.56 1.06 2.64

Stem + branch surface area (cm2) Stem + branch volume (cm3) Stem base diameter (mm) Mean pipe length (mm) Standard deviation of pipe length (mm) Total woody segment length (mm)

= = = = = =

2210 379 12.7 730 288 13100

= = = = =

0.0991 0.0994 0.0843 0.851 0.848

Expected distance to 5 nearest leaves (no edge Expected distance to 5 nearest leaves (with edge Observed distance to 5 nearest Dispersion parameter (no edge Dispersion parameter (with edge

corr.) corr.) leaves corr.) corr.)

The summary displays the following variables:

• Total leaf area Total plant leaf area (m2 ). • Mean leaf size Average area of individual leaves (cm2 ). • Number of leaves Total number of leaves. • Mean leaf length Average length of leaves (cm). • Mean leaf angle Average leaf inclination angle (0 = horizontal, 90 = vertical). • Mean leaf angle, weighted by leaf area Same, but weighted by leaf size (more appropriate in some, but not all cases!). • Ellipsoidal distribution parameter The parameter of the ellipsoidal leaf angle distribution fitted to this plant. • Crown volume The volume of the crown (dm3 ), as determined by the convex hull around the leaves. • Crown surface area The surface area (m2 ) of the same convex hull. • Stem and branch volume and surface area Total volume (m3 ) and surface area (m2 ) of woody segments in the crown. • Stem base diameter Diameter of the main stem (diameter of the first node) (mm). • Mean pipe length Mean length of the path from base of petiole to the base of the first node (mm). • Standard deviation of pipe length SD across leaves of the pipe length (mm). • Total woody segment length Total length of branch and stem sections (mm). • Expected distance to 5 nearest leaves Distance to 5 nearest leaves for randomly distributed leaves, for a square box with the same leaf number density as the current plant (with or without edge correction, see Details below). • Observed distance to 5 nearest leaves Actual mean distance between leaves and their 5 nearest neighbors.

10

• Dispersion parameter (with or without edge correction A parameter describing how random the leaf distribution within the crown hull is: less than 1 indicates more clumped, larger than 1 more uniform. • Crown width Twice the mean crown radius, where crown radius is calculated for four quadrants as the distance from the middle of the plant to the furthest leaf in that quadrant. • Crown length Vertical distance between lowest and highest leaf on the plant. • Crown shape index A parameter describing the shape (this may change in future versions!).

3.7 Plot the plant in 3D To plot a plant in 3D is easy:

plot(pilularis)

Shown here is a rather low resolution output of the plot. In an R session, This opens up a window with the plant plotted in 3D. The plant can be rotated with the mouse (left-click the plant and drag), zoomed in and out (with the scroll wheel on the mouse), and the window can be resized. The plot can be saved using the snapshot3d function in the rgl package. You can also rotate the plant around its z-axis (note that this continues until Escape is pressed):

play3d(spin3d()) See spin3d in the rgl package for more information, as well as movie3d to convert the spinning scene to a movie (an animated GIF, for example). The following code plots the maple with the convex hull included.

plot(pilularis, addcrownhull=TRUE)

11

See the help page plot.plant3d for more options. 3.7.1 Side and top views Simple 2D plots of the plant are produced with viewplot, like so:

viewplot(pilularis) View from above

400 200

800

0

Y

600

Z 0

0

−400

200

200

−200

400

400

Z

600

800

1000

View from South

1000

View from East

−600

−400

−200

0

200

400

600

Y

−600

−400

−200

0

X

200

400

−400

−200

0

200

400

X

This is a rather lame way to plot the leaves and branches together in 2D, because the branches are always plotted on top (whether or not they are visible). To plot leaves only, do:

viewplot(pilularis, stems=FALSE) See also section ( 3.9.1) on a more general function for projecting a plant from any direction.

3.8 Modifying virtual plants There are a number of tools to modify existing plant files. First, we can replace the leaf angles in a plant file using the replaceangles function. This function generates a new plant file. There are multiple ways to provide the leaf angles, for example replace all angles with the same value, say 45 degrees:

replaceangles(whichangle="An.3", pfile="SugarMaple.p", outputfile="Maple45.p", newangles=45) The whichangle argument specifies which of the angles to replace, ’An.3’ corresponds the leaf angle. The other angles are coded as: 12

• Az, An Azimuth and angle of stem sections • Az.1, An.1 Branch sections • Az.2, An.2 Petioles • Az.3, An.3 Leaves Another option is to draw a sample from a given leaf angle distribution (LAD). For this, the LeafAngle package is used. First, set up a LAD. A common distribution if the ellipsoidal LAD, which uses one parameter.

ellipsLAD <- angledist("ellipsoidal", distpars=2.2) This sets up a leaf angle distribution from which we can now draw samples (see ?drawsample). The replaceangles function can take the argument distobj, a leaf angle distribution object, like so:

replaceangles(whichangle="An.3", pfile="SugarMaple.p", outputfile="Maple45.p", distobj=ellipsLAD) Note that we can only replace leaf and petiole angles this way, not azimuth angles (because leaf angles are between 0 and 90, azimuth angles between 0 and 360). For more modifications of plant files, see these functions:

modifypfile changeinternodes

3.9 STARbar : Calculating the projected and displayed leaf area 3.9.1 Projecting a plant We can generate two-dimensional projections of the three-dimensional virtual plant canopy by using the projectplant function. Note: toona is a built-in plant, you don’t have to construct it.

toona_fromabove <- projectplant(toona, azimuth=0, altitude=90) This generates a projected plant object, in this case the Toona australis plant viewed from above (an altitude of 90 degrees). Typing the name of the object reminds you the plant files it was based on, and the viewing angles:

toona_fromabove ## ## ## ## ##

Yplant - object of class 'projectedplant3d'. -----------------------------Plant has 49 leaves. Projection from azimuth = 0 and altitude = 90

And, we can also make a simple plot of the projected plant. This plotting method has no special options at this time.

plot(toona_fromabove)

13

300 250 200 Y

150 100 50 0 −200

−150

−100

−50

0

50

100

X

3.9.2

Displayed and projected leaf area

The displayed and projected areas for a plant are calculated with the STARbar function. By default, this function calculates the hemi-spherically averages ratio of displayed area to total leaf area (known as STAR, the ’silhouette to total area ratio’). This average value is known as STAR (see Duursma et al. 2012). It is also possible to calculate STAR for one or many specified viewing angles (see ?STARbar for details).

STARbar(sugarmaple) ## Calculating STARbar with method:

gridtracer, integration scheme:

## ## Average displayed and projected leaf area. ## -----------------------------## Total leaf area (m2) ## Spherically averaged displayed leaf area (m2) ## Spherically averaged STAR (-) ## ## Spherically averaged projection (-)

Turtlesky.

= = =

0.1400225 0.05725343 0.4088874

=

0.4985828 (Should be 0.5).

The STARbar function also takes a list of plants as input. See section 3.16.4 for batch analysis.

3.10 setMet : Generating weather diurnal Details and examples will be added at a later date. See ?setMet for full details.

3.11 setPhy : Setting a leaf photosynthesis model The YplantQMC package has two built-in leaf gas exchange models, the coupled Farquhar-Ball-Berry model (with various stomatal conductance models), and a simple non-rectangular hyperbola. We already saw an example of how to use the latter model in section 2.5 on page 5.

14

The default stomatal conductance model is that of Medlyn et al. (2011) in its simplified form (their Figure 3). To use the Farquhar model, there are four required parameters, and more than two dozen optional parameters. See ?setPhy for the full list of parameters. More background can also be found on the MAESTRA website (http://bio.mq.edu.au/maestra/manual.htm) ; the Farquhar model as implemented in YplantQMC was taken directly from MAESTRA, and all optional parameters have the same name in both models. These are the required parameters: where: Vcmax Maximum carboxylation rate (mu mol m-2 s-1) Jmax Maximum electron transport rate (mu mol m-2 s-1) Rd0 Dark respiration rate at 25 deg C (or other reference temperature set by RTEMP) G1 Slope of the stomatal conductance model, see Medlyn et al. (2011) To use the Farquhar model, use this example:

eucphy <- setPhy("Farquhar", leafpars=list(Vcmax=50, Jmax=100, G1=8, G0=0.01, Rd0=1)) Note that we also set one of the optional parameters (G0, the stomatal conductance when PAR = 0). 3.11.1 Writing your own leaf gas exchange model Although YplantQMC has two built-in leaf gas exchange models, it is straightforward to implement your own. All you have to do is write an R function that follows a simple (flexible) template. The template is shown below:

leafmodel <- function(PAR, Tair, parameter1, parameter2, ...){ # Calculations here. return(data.frame(A=A)) } Where PAR is a required input, Tair is air temperature, one of the optional environmental inputs (see below), and parameter1 and parameter2 are parameters for this leaf gas exchange models (there is no limit to the number of parameters, and names can be changed to anything). Note the ... - it is crucial that this is included in the list of arguments. The parameters are set through setPhy, if the parameters do not have default values (as in the example above), each parameter has to be set. To use the above template, for example:

myphy <- setPhy("leafmodel", leafpars=list(parameter1=100, parameter2=0.5)) The following environmental variables are all passed to the leaf gas exchange model, it is up to you to decide if you use them (but they have exactly these names!): Tair Air temperature (deg C) Ca Atmospheric CO2 concentration (ppm) VPD Vapor pressure deficit (kPa) RH Relative humidity ( Patm Atmospheric pressure (kPa)

15

SWP Soil water potential (MPa) The output of the leaf gas exchange model has to be a data.frame, with at most (as of now), the following three variables: A Net CO2 assimilation rate (mu mol m-2 s-1) E Transpiration rate (mmol m-2 s-1) gs Stomatal conductance (mol m-2 s-1) Finally, here is the non-rectangular hyperbola as implemented in YplantQMC :

# Non-rectangular hyperbola lightresponse <- function(PAR, Amax, phi, theta, Rd, ...){ A <- (phi*PAR+Amax-((phi*PAR+Amax)^2-4*theta*phi*PAR*Amax)^0.5)/(2*theta) - Rd return(data.frame(A=A)) }

3.12

setHemi : Reading a hemispherical photograph

The YplantQMC model can read a processed hemispherical photograph (’hemiphoto’), where the canopy gap fraction is averaged by viewing altitude (or zenith) angle, and azimuth (compass) angle. The following box shows an example hemiphoto (this can be a space-delimited, or comma delimited file).

Azimuth Zenith 0 45 90 135 180 225 270 2.25 0.375 0.701 0.408 0.338 0.355 0.041 0.003 0.205 6.75 0.011 0.641 1.000 0.473 0.023 0.000 0.000 0.061 11.25 0.010 0.460 0.936 0.400 0.008 0.005 0.000 0.077 15.75 0.038 0.282 0.335 0.371 0.048 0.023 0.040 0.058 20.25 0.090 0.160 0.045 0.064 0.019 0.000 0.005 0.086 24.75 0.025 0.073 0.014 0.108 0.000 0.000 0.000 0.124 29.25 0.000 0.045 0.005 0.048 0.007 0.000 0.002 0.099 33.75 0.000 0.009 0.004 0.020 0.058 0.000 0.004 0.018 38.25 0.001 0.001 0.019 0.044 0.063 0.000 0.005 0.079 42.75 0.000 0.038 0.071 0.025 0.027 0.000 0.000 0.060 47.25 0.000 0.061 0.046 0.001 0.000 0.000 0.001 0.021 51.75 0.000 0.014 0.048 0.060 0.001 0.000 0.000 0.000 56.25 0.000 0.000 0.022 0.039 0.000 0.000 0.000 0.000 60.75 0.013 0.002 0.000 0.000 0.000 0.000 0.000 0.000 65.25 0.000 0.004 0.000 0.000 0.000 0.000 0.000 0.000 69.75 0.000 0.005 0.001 0.000 0.000 0.000 0.000 0.000 74.25 0.000 0.001 0.012 0.000 0.000 0.000 0.000 0.000 78.75 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 83.25 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 87.75 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000

315

The first row of the file will not be read (so you can add a label ’Azimuth’ as shown above, although it is not

16

necessary). The second row lists the bins of azimuth angles for which gap fractions are measured. These bin values are not read from the file - is is simply assumed that the bins are evenly spaced, and start at 0 (North), and proceed clockwise (N, E, S, W). Recall that the azimuth angle is defined as North=0, East=90, South=180. In the first column, the first cell is empty, followed by either ’Zenith’ or ’Altitude’. This label will determine whether the bins are in order of zenith angle, or altitude angle. As with the azimuth angle, the bin ranges are not read from the file (the values in the first column), it is instead assumed that they are evenly spaced (starting at 0, ending at 90 degrees). The following code reads an example file, and makes a standard plot. In the second plot, the solar path is added from a yphemi object.

# A built-in hemiphoto (you don't need to read this in). # smallgap <- setHemi("canopygap.csv") plot(smallgap) # A larger gap, with solar path: southernFrance <- setLocation(44) sunnyday <- setMet(southernFrance, nsteps=12, Tmin=9, Tmax=25, month=6, day=21) plot(largegap, sunnyday)

N

W

N

E

W

Gap fraction 0 0.25 0.5 0.75 1

3.13

E

Gap fraction 0 0.25 0.5 0.75 1

S

S

YplantDay : Running a daily simulation

This is the main function of the YplantQMC package, and also the most user-friendly. For more customizable simulations, it is recommended to use the runYplant function. The YplantDay function runs the ray-tracer and the leaf gas exchange model for every timestep as defined in the daily weather object (the met object). For example,

# Set location, southernfrance <- setLocation(lat=44) # A daily weather object, use a constant beam fraction of 0.4. sunnyday <- setMet(southernfrance, month=6, day=21, nsteps=12, Tmin=9, Tmax=29, PARday=22, fbeamday=0.4, fbeammethod="constant") # Light response curve:

17

toonalrc <- setPhy("lightresponse", leafpars=list(Amax=14.5, Rd=1.4, phi=0.05, theta=0.5, reflec=0.1, transmit=0.05)) # Run YplantQMC for a day. Use the built-in 'largegap' hemiphoto. toonarun <- YplantDay(toona, sunnyday, toonalrc, largegap) The simulation will take less than a minute to complete, and a fraction of that time when only diffuse radiation is calculated for every timestep (because the relative diffuse radiation absorption does not change with time of day, so QuasiMC is only run once). To generate a standard PDF report, which summarizes all YplantQMC objects, use this command:

ypreport(toona, phy=toonalrc, met=sunnyday, hemi=largegap, ypsim=toonarun) ## Warning in ypreport(toona, phy = toonalrc, met = sunnyday, hemi = largegap, : object and plant used for simulation are not the same!

**-Plant

## Report generation probably successful. ## See file ’YplantQMC report 2014-12-11.pdf’ in directory : ## G:/Work/Projects/YPLANT/Yplant in R/YplantQMC Vignette Note that this might fail if you have a PDF open with the same name (as you may have generated this report once before on the same day). Make sure to close the PDF before running ypreport. If a PDF is generated that can’t be opened, use this command:

dev.off() And try again. By default, a .PSR file is written to disk, which contains totals and averages by timestep for the simulation. The PSR data can also be viewed with the command:

psrdata(toonarun) ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ## ##

timeofday PAR0 PARinc PARleaf PARdir 4.934376 109.3673 25.93051 43.11843 34.18778 6.224785 538.8108 127.74976 104.35947 60.36155 7.515195 1034.9948 456.26583 306.99106 149.84989 8.805604 1483.5091 781.64084 471.47934 202.27673 10.096013 1818.4194 1082.41389 613.58905 240.79809 11.386422 1996.8433 1223.96506 681.80524 260.26305 12.676831 1996.8433 857.17521 467.20977 171.99259 13.967241 1818.4194 719.92781 397.54553 149.59732 15.257650 1483.5091 581.23879 328.55901 128.37627 16.548059 1034.9948 285.81357 169.68808 71.25187 17.838468 538.8108 127.74976 93.28365 49.28573 19.128877 109.3673 25.93051 43.03357 34.10292 timestep LAplant LAproj LAsunlit 1 4645.473 0.0266781 0.01202210 0.009003645 2 4645.473 0.0266781 0.01329398 0.011204353 3 4645.473 0.0266781 0.01447754 0.013253675 4 4645.473 0.0266781 0.01559476 0.014189962 5 4645.473 0.0266781 0.01655462 0.014580186 6 4645.473 0.0266781 0.01668066 0.015121874 7 4645.473 0.0266781 0.01575962 0.014272768 8 4645.473 0.0266781 0.01443521 0.013578982 9 4645.473 0.0266781 0.01278929 0.012085962 10 4645.473 0.0266781 0.01132061 0.010038018 1 2 3 4 5 6 7 8 9 10 11 12

18

PARdiff 8.930651 43.997917 157.141167 269.202610 372.790968 421.542187 295.217181 247.948207 200.182732 98.436205 43.997917 8.930651

A 0.5430344 2.8313835 7.0626325 8.7254833 9.6046084 9.9018687 8.6846681 8.1035212 7.3886929 4.6402706 2.4613369 0.5379535

A0 -0.1613238 3.6429304 8.8818974 10.4968407 11.1912418 11.4056766 10.7136132 10.2893018 9.6835416 7.0320564 3.6429304 -0.1613238

## 11 4645.473 0.0266781 0.01073435 0.009136927 ## 12 4645.473 0.0266781 0.01040365 0.008989465 This returns a dataframe, with the following variables:

timeofday Time of day (hours) of timestep. PAR0 Incident PAR above canopy (mu mol m-2 s-1) PARinc Incident PAR below canopy (mu mol m-2 s-1) PARleaf Absorbed total PAR (mu mol m-2 s-1) PARdir Absorbed direct PAR (mu mol m-2 s-1) PARdiff Absorbed diffuse PAR (mu mol m-2 s-1) A CO2 assimilation rate (mu mol m-2 s-1) A0 A by horizontal leaf above canopy (mu mol m-2 s-1) E Transpiration rate (mmol m-2 s-1) timestep Length of timestep (seconds) LAplant Total plant leaf area (m2) LAproj Total projected leaf area (m2) (’NA’ if diffuse PAR only). LAsunlit Total displayed (sunlit) leaf area (m2) (’NA’ if diffuse PAR only). Note that if E is missing, the leaf gas exchange model used in the simulation did not calculate transpiration rate, and when LAproj and LAsunlit are missing, there was zero direct radiation for every timestep (because these variables are meaningful for direct radiation only). We can make a standard plot of the simulation. Note that in this example, photosynthesis is not plotted because we did not calculate it (because in the above example code, we did not include a phy object).

4

6

8

10

Plant Hor. leaf

2

A plant (µmol m−2(leaf) s−1)

12

plot(toonarun)

0

5

10

15

20

Time of day (hours)

19

400 300 200 0

100

APAR (µmol m−2 s−1)

500

1000

1500

2000

Direct Diffuse

0

APAR (µmol m−2 s−1)

Plant Hor. leaf Above canopy

0

5

10

15

20

0

Time of day (hours)

5

10

15

20

Time of day (hours)

Finally, a dataframe with all leaf-by-leaf simulation results is stored in the simulation object as well. This dataframe is stored as toonarun$outdata, and can be saved to disk either at simulation time,

toonarun <- YplantDay(toona, sunnyday, toonalrc, largegap, writeOUT=TRUE) which produces a TAB-delimited text file (with the name, in this case, ”Toona-YplantQMC.OUT”). Or using the standard R function write.table, for example:

write.csv(toonarun$outdata, "Toona.OUT", sep=",") which produces a comma-delimited output file.

3.14

runYplant : A single YplantQMC run for custom simulations

The runYplant function runs one simulation of the YplantQMC model. There are no summary, plot or print methods for the result, the function simply returns a dataframe with results for each leaf. Useful function, because the solar position (altitude, azimuth), beam fraction, and individual weather variables can be set directly as input arguments, rather than through defining objects as in YplantDay. This makes runYplant ideal for any custom simulation, sensitivity analyses, response curves, and so on. The only exception is the leaf physiology object, which still has to be passed as an object defined by setPhy. To run multiple simulations, and storing the results in a list, consider the following general setup:

# Define an empty list ypresults <- list() # Define a dataframe where each row is a separate simulation. # In this case, a number fo combinations of fbeam and solar altitude values. rundfr <- expand.grid(altitude=seq(30,90,by=15), fbeam=seq(0,1,by=0.2)) # Run the model for each row of the rundfr defined above: for(i in 1:nrow(rundfr)){ ypresults[[i]] <- runYplant(myplant, myphy, altitude=rundfr$altitude[i], fbeam=rundfr$altitude[i]) }

20

The resulting list ypresults is quite long, and best further analyzed with the lapply and/or sapply functions. See ?runYplant for full details on the arguments and options.

3.15

QuasiMC options and settings

Most of the settings of the ray-tracer QuasiMC can be set through YplantQMC as well. See ?runYplant for a few options; the full list is set in the invisible function writecfg, which can be displayed with page(YplantQMC:::writecfg). A description of the various options will appear in this section.

3.16

Batch processing

For the keen R programmer, it should be obvious already that YplantQMC can be used quite easily to perform all kinds of batch analyses and custom simulations, using for loops, lists, and other standard tools. For convenience, a few built-in functions make batch analyses very straightforward.

3.16.1 Reading and storing multiple plants The readplantlist function is used to construct a list of plant3d objects, which can then be used in several functions for convenient analysis of many plants in one go.

myplants <- readplantlist(pfiles=c("Epilularis.p","Sugarmaple.p","Toona.p"), lfiles=c("Epilularis.l","Sugarmaple.l","Toona.l")) ## Plant 1 of 3 ( Epilularis.p ) constructed in 0.46 sec. ## Plant 2 of 3 ( Sugarmaple.p ) constructed in 0.05 sec. ## Plant 3 of 3 ( Toona.p ) constructed in 0.08 sec. This object myplants is a list of plants, so that you can access the first plant by myplants[[1]], and so on. For most analyses, I recommend the following setup. Create a directory with all your plant and leaf files, and a file called leafplantkey.txt, which links the plant to the leaf files. An example file looks like this:

acaflo1.p,acaflo.l acaflo2.p,acaflo.l acaflo3.p,acaflo.l acamyr1.p,acamyr.l acamyr2.p,acamyr.l acamyr3.p,acamyr.l Then, read in all your plants with this command:

myplants <- readplantlist(lpk="leafplantkey.txt") This setup guarantees you have the leaf files correctly linked to the corresponding plant files.

21

3.16.2 Plotting multiple plants For example, plot the first of your plants:

plot(myplants[[1]]) To plot each of the plants in the myplants object, and save a PNG snapshot, use this command:

plot(myplants, png=TRUE) This produces a series of .PNG files in the current working directory (see getwd() if you can’t find your files).

3.16.3 Summarizing multiple plants To summarize all the plants, and write a text file with the output,

summary(myplants, writefile=TRUE) See section 3.6 for details on the output, or the help page ?summary.plant3d.

3.16.4 Calculating displayed area for multiple plants Using the myplants object constructed above, we first estimate STAR for each plant, and then convert the results to a dataframe.

starmyplants <- STARbar(myplants) ## Plant 1 of 3 done in 2.06 sec. ## Plant 2 of 3 done in 0.3 sec. ## Plant 3 of 3 done in 0.36 sec. as.data.frame(starmyplants) ## STARbar DAbar PAbar LA PALAbar pfile lfile ## 1 0.3470306 186496.00 269177.70 537405.0 0.5008842 Epilularis.p Epilularis.l ## 2 0.3493017 48910.08 69812.79 140022.5 0.4985828 Sugarmaple.p Sugarmaple.l ## 3 0.3834896 10230.77 13314.23 26678.1 0.4990696 Toona.p Toona.l Where STARbar is STAR, DAbar is the spherically averaged displayed area, PAbar the projected area (, LA the plant leaf area, PALAbar the average projected area fraction (which should be 0.5, +/- some numerical error), and pfile and lfile the files used to construct the plants. Note that if you are only interested in STAR, it can also be calculated with,

summary(myplants, calcSTARbar=TRUE, writefile=TRUE)

3.16.5 Running YplantDay for multiple plants The YplantDay function has methods for both single plants (objects of class plant3d) and lists of plants objects of class plant3dlist). To run one single plant of your plant list object, do this:

run1 <- YplantDay(myplants[[2]], mymet) This will run a simulation for the second plant in your list. To simulate all your plants, simply run the following command:

22

allruns <- YplantDay(myplants, sunnyday) The result of this function (which we stored in the object allruns) is a list of YplantDay result objects. Each of the functions described in section 3.13 can thus be used in the same way, for example:

psrdata(allruns[[1]]) will display the results for the first plant. There is also a convenient function that produces daily totals of APAR and photosynthesis, transpiration (if they were calculated), like this:

allsumm <- summary(allruns) This function returns a dataframe, use the following command to save it to also write a text (TAB-delimited) file to the current working directory, including a header that explains the variables.

summary(allruns, writefile=TRUE) To combine the results of the summary of your plants, and the YplantQMC simulation results, simply use cbind :

combresults <- cbind(summary(myplants), summary(allruns)) and use write.csv to produce text output, for example.

4 Unimplemented features of YplantQMC This is a list of the features that are not yet implemented, but are likely to be implemented at some point, provided there is sufficient need:

• The QuasiMC ray-tracer separates incident PAR by lower and upper surfaces of the leaf, but this is not yet available through YplantQMC . (High priority). • Leaf temperature is not yet calculated. (High priority). • Penumbra effects : direct solar radiation from a solar disk, rather than a point source (High priority). • Plasticity in leaf physiology : ability to specify leaf physiology parameters by leaf, by top/bottom of leaf, or as a function of some other variable (height in the plant, relative shading, etc.) (High priority). • Hydraulics (i.e., calculation of stem and leaf water potential at various points). (Medium/high priority) • Sunfleck information read from hemiphoto analyses (as in Yplant). (Medium priority). • Currently, stem and branch sections do not cast shadows. (Low priority). (Very low priority).

5 Frequently Asked Questions (FAQ) Why is YplantQMC slower than the original Yplant? The ray-tracer that is used, QuasiMC, is much more advanced than the ray-tracer in the original Yplant, and also more accurate. The QuasiMC ray-tracer calculates multiple scattering and re-absorption of reflected and transmitted solar rays, which is computationally intensive (Cieslak et al. 2008). Also, for very large plants, YplantQMC is actually faster than the original Yplant.

How do I write any output to a text file?

23

Some functions have a built-in option to write (commented) text output files, but perhaps you want to write any text output to a file. Use the sink function, for example:

sink("Summary output.txt") summary(toona) sink() Will write the output of summary(toona) to a text file, rather than to the console.

Why is there no graphical user interface? Because I have not implemented one, nor do I plan on doing so. It is in principle possible to write a GUI with the Tcl/Tk toolkit, or more simply with the winMenuAdd utility (Windows only).

6 Acknowledgments Thanks to Bob Pearcy for sharing the original code for Y-plant, and encouraging a new implementation of the model. Also thanks to Daniel Falster for many useful discussions, and to Andrew Haigh for developing the ypConvert python toolkit.

7 References Cieslak, M., C. Lemieux, J. Hanan and P. Prusinkiewicz. 2008. Quasi-Monte Carlo simulation of the light environment of plants. Functional Plant Biology. 35:837-849. Duursma, R.A., D.S. Falster, F. Valladares, F.J. Sterck, R.W. Pearcy, C.H. Lusk, K.M. Sendall, M. Nordenstahl, N.C. Houter, B.J. Atwell, N. Kelly, J.W.G. Kelly, M. Liberloo, D.T. Tissue, B.E. Medlyn and D.S. Ellsworth. 2012. Light interception efficiency explained by two simple variables: a test using a diversity of small- to medium-sized woody plants. New Phytologist. 193:397-408. Medlyn, B.E., R.A. Duursma, D. Eamus, D.S. Ellsworth, I.C. Prentice, C.V.M. Barton, K.Y. Crous, P. De Angelis, M. Freeman and L. Wingate. 2011. Reconciling the optimal and empirical approaches to modelling stomatal conductance. Global Change Biology. 17:2134-2144. Pearcy, R.W. and W. Yang. 1996. A three-dimensional crown architecture model for assessment of light capture and carbon gain by understory plants. Oecologia. 108:1-12. Pearcy, R.W., R.A. Duursma, D.S. Falster. 2011. Studying plant architecture with Y-plant and 3D digitising. Prometheus Wiki. http://tinyurl.com/6r2wn2q

24