In chgigot/canis: Canopy File Investigator and Shreder

Import a `.can` file

library(knitr)
library(rgl)
knitr::opts_chunk$set(echo = TRUE, cache = TRUE)
knit_hooks$set(rgl = hook_rgl)

library(canis)

We need to use the function readCan to import a .can file in the proper way.

library(canis)
field1 <- readCan(system.file("extdata", "field1.can", package = "canis"))

field1 contains r nrow(field1) rows. Each row corresponds to a polygon (i.e. a triangle here) of the virtual canopy. In the rest of this document, we will use the plants data set instead of field1. plants contains a bit less triangles (r nrow(plants)) and represents a sparse wheat canopy with random plant locations.

head(plants, 3)
str(plants)

We see that the first r ncol(plants) - 1 columns correspond to different IDs or labels. Contrary to the other ones, the last column (vertices) is actually a so-called list-column. A key advantage of list-columns is a structural independance from the number of vertices each polygon is made of. The package canis could therefore be extented to deal with many sorts of polygons (and not only triangles) without the need to totally redesign the interface.

If we glance at an element of the column vertices, we can see that it is a matrix. Each row of this matrix corresponds to the coordinates of a vertex of the polygon (in a Cartesian coordinate system). The column names (x, y and z) refer to the three spatial dimensions.

is(plants$vertices[[36]])
plants$vertices[[36]]

Display in 3D

The plants data set can be displayed in 3D using display3d. The color parameter may refer to a color (a color name, an output of the rgb function,...) or the name of a column of plants (except vertices). Type ?canis::display3d for more information.

view3d(theta = 0, phi = -75, zoom = 0.65)

display3d(plants, col = "green4", alpha = 0.5)

view3d( theta = 0, phi = -75, zoom = 0.65)

display3d(plants, col = "plantID")

view3d(theta = 0, phi = -75, zoom = 0.65)

As the 3D plotting process is especially time-consuming, the use of the parameter fraction may be relevant when working with a large number of polygons. A random fraction (corresponding to the value of fraction, between 0 and 1) of all the polygons is selected to be displayed.

display3d(plants, col = "green4", fraction = 0.2)

Split the canopy in accordance with a set of volumes

Design the set of volumes to play with

We need to create a set of volumes whose all sides will define the boundaries to chop the canopy accordingly. To do so, we use the dedicated function (i.e. makeVolumeSet) which allows to build a set of volumes based on a regular grid.

rangeDim <- function(polygon, dim) {
    range(sapply(polygon$vertices, function(i) i[, dim]))
}
x <- rangeDim(plants, "x")
y <- rangeDim(plants, "y")
z <- rangeDim(plants, "z")
volumes <- makeVolumeSet(x, y, z, intervals = c(4, 1, 4), label = "auto")
head(volumes)
str(volumes)

Let's display this set of volumes above the canopy:

view3d(theta = 0, phi = -75, zoom = 0.65)

display3d(plants, col = "green4")
display3d(volumes, col = "volumeID", alpha = 0.3, label = "levels")

Allocate a volume to each polygon (and shred polygons if needed)

Once in possession of this set of volumes, we can allocate each polygon of the virtual canopy to a volume using the function splitCanopy. To shred each polygon crossing a side of a volume, the parameter shred must be set to TRUE (default).

Note: In its current implementation, canis does not take into account overlapping of volumes. This means that if, for a given polygon, the first volume handled by the algorithm contains the polygon, it's done: canis will not check if this polygon is also included in another volume. So be careful, and remember: Visualize, visualize and visualize again your scenes before going further ahead!

newPlants <- splitCanopy(plants, volumes, shred = TRUE,
                         parallel = TRUE, nCores = 8)

That's it. Let's check that we have a new column named volumeID.

head(newPlants, 3)
str(newPlants)

view3d( theta = 0, phi = -75, zoom = 0.65)

display3d(newPlants, col = "volumeID")
display3d(volumes, col = "grey", alpha = 0.1)

Note that depending on your computer and the scene complexity (i.e. the number of polygons and volumes), 3D display may take a few tens of seconds.

And now?

Save the output of `splitCanopy` in a file

You may want to save the data frame returned by splitCanopy in a text file (e.g. a .csv file) for backup reasons or to use it in another software. To do so, you first need to flatten the data frame in order to get rid of the list-column structure. The function flatten does the job and return a regular data frame that you can seamlessly save in a text file.

flatNewPlants <- flatten(newPlants)
filename <- tempfile(fileext = ".csv")
write.csv(flatNewPlants, filename, row.names = FALSE)

file.remove(filename)

Perform subsequent analyses in R

As we are already in a nice and rich working environment, let's try to plot a couple of graphes of our data set. We first need to import a few helpful R packages and write some functions.

library(dplyr)
library(tidyr)
library(ggplot2)

## A set of useful functions
# Area of a triangle computed using Heron's formula
heron <- function(triangle) {
    distances <- as.vector(dist(triangle))
    s <- sum(distances) / 2
    m <- s * (s - distances[1]) * (s - distances[2]) * (s - distances[3])
    # Even if it's not supposed to occur in theory, we need to check that
    # m > 0. Indeed, due to machine's numerical precision, we can get negative
    # numbers sometimes.
    if (m < 0) return(0)
    else       return(sqrt(m))
}
# Sum of all the areas of a list of triangles
area <- function(triangles) sapply(triangles, function(x) heron(x))

Now, let's go:

myData <- newPlants %>%
    filter(leafID != 0) %>% # Get rid of the internodes
    mutate(area = area(vertices))
mySubData1 <- myData %>%
    group_by(plantID, leafID, volumeID) %>%
    summarise(leafAreaInVol = sum(area))
mySubData2 <- myData %>%
    group_by(plantID, leafID) %>%
    summarise(leafArea = sum(area))
myData <- left_join(mySubData1, mySubData2, by = c("plantID", "leafID")) %>%
    mutate(pcentLeafAreaInVol = leafAreaInVol / leafArea) %>%
    select(-leafArea) %>%
    separate(volumeID, into = c("x", "y", "z"), sep = "-") %>%
    as.data.frame()
rm(mySubData1, mySubData2)
head(myData)

# Leaf area density in a 2D x-z space
myDataA <- myData %>%
    group_by(x, z) %>%
    summarise(leafAreaInVol = sum(leafAreaInVol))

myDataA$x <- as.factor(myDataA$x)
levels(myDataA$x) <- c("left", "middle-left", "middle-right", "right")
myDataA$z <- as.factor(myDataA$z)
levels(myDataA$z) <- c("low", "middle-low", "middle-high", "high")

g <- ggplot(myDataA, aes(x = x, y = z)) +
    geom_raster(aes(fill = leafAreaInVol)) +
    scale_fill_gradient(low = "white", high = "green4",
                        name = expression("Leaf area (" * cm^2 * ")")) +
    theme_classic() +
    labs(x = "Canopy side", y = "Canopy height",
         title = "Leaf area density in a 2D x-z space")
print(g)

# Leaf area density in a 2D "leaf level-z space"
test6 <- myData %>%
    group_by(leafID, z) %>%
    summarise(leafAreaInVol = sum(leafAreaInVol))

test6$z <- as.factor(test6$z);
levels(test6$z) <- c("low", "middle-low", "middle-high", "high")

g <- ggplot(test6, aes(x = leafID, y = z)) +
    geom_raster(aes(fill = leafAreaInVol)) +
    scale_fill_gradient(low = "white", high = "green4",
                        name = expression("Leaf area (" * cm^2 * ")")) +
    theme_classic() +
    labs(x = "Leaf level", y = "Canopy height",
         title = "Leaf area density in a 2D \"leaf level-z space\"")
print(g)