Image phenotyping platforms

Image phenotyping platforms are a fast and non-invasive tool for phenotyping plant under lab and field conditions. The most common are sensors for imaging. The reflected light captured in these images can be used to draw inference about many traits, including:

• Geometric traits (i.e. plant height, leaf area index, lodging, crop canopy cover)
• Canopy spectral texture (spectral features)
• Physiological traits (i.e., chlorophyll, biomass, pigment content, photosynthesis)
• Abiotic/biotic stress indicators (i.e., stomatal conductance, canopy temperature difference, leaf water potential, senescence index)
• Nutrients (nitrogen concentration, protein content)
• Yield

Image segmentation

Image segmentation is the process of partitioning a digital image into multiple segments (sets of pixels or image objects).

Thresholding

Thresholding is the simplest method of image segmentation. This method is based on a clip-level (or a threshold value) to turn a gray-scale image into a binary image (Mask).

FIELDimageR

FIELDimageR is a R package to analyze images from research fields, and allows to:

• Crop the image
• Remove the background
• Build vegetation indices
• Rotate the image
• Build the plot shapefile
• Extract information for each plot

Required packages

• FIELDimageR

# Installing 
install.packages("sp")
install.packages("raster")
install.packages("rgdal")
install.packages("ggplot2")
install.packages("agricolae")
install.packages("reshape2")
install.packages("devtools")
install.packages("parallel")
install.packages("foreach")
install.packages("doParallel")
install.packages("lme4")
devtools::install_github("filipematias23/FIELDimageR")

# Necessary packages 
library(FIELDimageR)
library(raster)
library(ggplot2)
library(agricolae)
library(reshape2)
library(lme4)

A) Counting green seeds in Soybean:

• 4 pictures: soybean samples with different green seeds percentage

• Number of green seeds (0, 1, 6, and 11)

• Donwload: soybean.zip

 # Uploading one image as example and decreasing the resolution
EX1<-stack("soybean/11.jpg")
EX1<-aggregate(EX1, fact= 4)
plotRGB(EX1)

# extent=T (The whole image area will be the shapefile)
EX.shapeFile<-fieldPolygon(EX1,extent = T)

# Select one index to identify leaves and remove the background
EX1.I1<- fieldIndex(mosaic = EX1,index = c("SI","BGI","BI"))

# Thresholding
dev.off()
par(mfrow=c(1,2))
hist(EX1.I1$BGI)
plot(EX1.I1$BGI)

# Removing the background
EX1.R1<- fieldMask(mosaic = EX1, index = "BGI",
                   cropValue = 0.7,
                   cropAbove = T)

# Counting the total number of seeds
EX.P.Total<-fieldCount(mosaic = EX1.R1$mask, 
                       fieldShape = EX.shapeFile$fieldShape, 
                       minSize = 0.1,
                       cex = 1.5,
                       na.rm = T)

# Select one index to identify green seeds
EX1.I2<- fieldIndex(mosaic = EX1.R1$newMosaic,index = c("SI","BGI","BI"))

#BI index
plot(EX1.I2$BI)

# Selecting green seeds
EX1.R2<- fieldMask(mosaic = EX1.R1$newMosaic, 
                   index = "BI",
                   #myIndex = "Blue",
                   cropValue = 130,
                   cropAbove = T)

# Counting the number of green seeds
EX.Green<-fieldCount(mosaic = EX1.R2$mask, 
                     fieldShape = EX.shapeFile$fieldShape, 
                     minSize = 0.07,
                     cex = 1.5,
                     na.rm = T)

# Joying information
data.frame(Total=EX.P.Total$fieldCount,
           Green=EX.Green$fieldCount,
           Percentage=round(EX.Green$fieldCount/EX.P.Total$fieldCount,2))

################
### Parallel ###
################

# Required packages
library(parallel)
library(foreach)
library(doParallel)

# Images names (folder directory: "./soybean/")
pics<-list.files("./soybean/")

# Number of cores
n.core<-2

# Starting parallel
cl <- makeCluster(n.core, output = "")
registerDoParallel(cl)
EX.Table.Parallel <- foreach(i = 1:length(pics), .packages = c("raster","FIELDimageR"), 
                               .combine = rbind) %dopar% {
EX1<-stack(paste("./soybean/",pics[i],sep = ""))
EX1<-aggregate(EX1, fact= 4)
EX.shapeFile<-fieldPolygon(EX1,extent = T,
                           plot = F)
EX1.R1<- fieldMask(mosaic = EX1, index = "BGI",
                   cropValue = 0.7,
                   cropAbove = T,
                   plot = F)
EX.P.Total<-fieldCount(mosaic = EX1.R1$mask, 
                       fieldShape = EX.shapeFile$fieldShape, 
                       minSize = 0.1,
                       cex = 1.5,
                       na.rm = T)
EX1.R2<- fieldMask(mosaic = EX1.R1$newMosaic, 
                   index = "BI",
                   cropValue = 130,
                   cropAbove = T,
                   plot = F)
EX.Green<-fieldCount(mosaic = EX1.R2$mask, 
                     fieldShape = EX.shapeFile$fieldShape, 
                     minSize = 0.07,
                     cex = 1.5,
                     na.rm = T)
data.frame(Total=EX.P.Total$fieldCount,
           Green=EX.Green$fieldCount,
           Percentage=round(EX.Green$fieldCount/EX.P.Total$fieldCount,2))
}
rownames(EX.Table.Parallel)<-pics
EX.Table.Parallel 

B) Vegetation indices and plant height in sorghum:

• 256 sorghum genotypes: granifero, biomassa, forrageiro
• 288 total plots (16 columns and 18 rows)
• Trait: plant height (PH)

• Donwload: orthophoto.tif

# Uploading file (orthophoto.tif):
EX1<-stack("orthophoto.tif")
plotRGB(EX1)

 # Rotating the image with theta -20:
EX1.Rotated<-fieldRotate(mosaic = EX1,theta = -20) # theta = -20

# Removing the NAs:
EX1.Crop <- fieldCrop(mosaic = EX1.Rotated) 

# Choosing an index to remove soil:
EX1.Indices<- fieldIndex(mosaic = EX1.Crop, index = c("BGI","NGRDI","SCI"))
dev.off()
par(mfrow=c(1,2))
hist(EX1.Indices$NGRDI)
plot(EX1.Indices$NGRDI)

# Removing soil (example using NGRDI):
EX1.RemSoil<- fieldMask(mosaic = EX1.Crop, 
                        index = "NGRDI",
                        cropValue = -0.1,
                        cropAbove = F)

Field data table with traits:

• Donwload: Data.txt

# Dataset - Field notes - phenotype (.exel or .txt):
Data<-read.csv("Data.txt",header = T, sep = "\t")

# Field map ID identification (ID="Plot"):
Map<-as.matrix(fieldMap(fieldPlot = Data$Plot,fieldColumn = Data$Column, fieldRow = Data$Row,decreasing = T))

# Building the plot shapefile (ncols = 45 and nrows = 5)
x11()
EX1.Shape<-fieldShape(mosaic = EX1.RemSoil$newMosaic,
                      ncols = 16, 
                      nrows = 18,
                      fieldData = Data,
                      ID = "Plot",
                      fieldMap = Map)

plotRGB(RGB.rescale(EX1.RemSoil$newMosaic,3))
plot(EX1.Shape$fieldShape,add=T)

# Building indices ("NGRDI","BGI","GLI","SCI")
EX1.Indices<- fieldIndex(mosaic = EX1.Shape$cropField, 
                         Blue = 3, Green = 2, Red = 1, 
                         index = c("NGRDI","BGI"), 
                         myIndex = c("(Red-Blue)/Green"))

EX1.Info<- fieldInfo(mosaic = EX1.Indices[[c("NGRDI","BGI","myIndex")]],
                     fieldShape = EX1.Shape$fieldShape,
                     buffer = -0.05,
                     n.core = 3)
NewData<-EX1.Info$fieldShape@data
NewData

# Making map plots
fieldPlot(fieldShape=EX1.Info$fieldShape,
          fieldAttribute="NGRDI", 
          mosaic=EX1.Indices, 
          color=c("white","black"), 
          alpha = 0.6,
          round = 2)

fieldPlot(fieldShape=EX1.Info$fieldShape,
          fieldAttribute="PH", 
          mosaic=EX1.Indices, 
          color=c("red","blue"), 
          alpha = 0.4,
          round = 2)

Estimated plant height using digital surface model (DSM)

• Donwload: dsm_0.tif
• Donwload: dsm_1.tif

# Uploading files from soil base (dsm_0.tif) and vegetative growth (dsm_1.tif):
DSM0 <- stack("dsm_0.tif")
DSM1 <- stack("dsm_1.tif")

dev.off()
par(mfrow=c(1,2))
plot(DSM0)
plot(DSM1)

# Rotating the image using the same theta (-20):
DSM0.R<-fieldRotate(DSM0, theta = -20)
DSM1.R<-fieldRotate(DSM1, theta = -20)

# Cropping the image using the previous shape:
DSM0.C <- fieldCrop(mosaic = DSM0.R,fieldShape = EX1.Crop)
DSM1.C <- fieldCrop(mosaic = DSM1.R,fieldShape = EX1.Crop)

# Canopy Height Model (CHM):
DSM0.New <- resample(DSM0.C, DSM1.C)
CHM <- DSM1.C-DSM0.New
dev.off()
plot(CHM)

# Removing the soil using mask from step 4:
CHM.S <- fieldMask(CHM, mask = EX1.RemSoil$mask)

# Extracting the estimate plant height average (EPH):
EPH <- fieldInfo(CHM.S$newMosaic, fieldShape = EX1.Info$fieldShape, fun = "mean", n.core=3)
colnames(EPH$fieldShape@data)[14]<-"EPH" # Changing name "layer" to "EPH"
EPH$plotValue

# Extracting the estimate plant height quantile (EPH.Q=0%,25%,50%,75%,100%):
EPH.Q <- fieldInfo(CHM.S$newMosaic, fieldShape = EPH$fieldShape, fun = "quantile", n.core=3)
colnames(EPH.Q$fieldShape@data)[c(16:20)]<-c("0%","25%","50%","75%","100%") # Changing names
EPH.Q$fieldShape@data

# Extracting the estimate plant height without soil base correction:
DSM1.S <- fieldMask(DSM1.C, mask = EX1.RemSoil$mask)
EPH.T <- fieldInfo(DSM1.S$newMosaic, fieldShape = EPH.Q$fieldShape, fun = "mean", n.core=3)
colnames(EPH.T$fieldShape@data)[c(22)]<-c("EPH.DSM1") # Changing names
EPH.T$fieldShape@data

NewData2<-EPH.T$fieldShape@data[,c("PH","EPH","0%","25%","50%","75%","100%","EPH.DSM1","NGRDI","BGI","myIndex")]
cor1<-correlation(NewData2)
cor1$correlation
cor1$pvalue

# Regression
NewData3<-melt(EPH.T$fieldShape@data[,c("Plot","Genotype","PH","EPH","75%","EPH.DSM1","NGRDI")],
               value.name = "Index",
               measure.vars = c("EPH","NGRDI","EPH.DSM1","75%"))
NewData3$variable<-factor(NewData3$variable,level=c("EPH","75%","EPH.DSM1","NGRDI"))
NewData3$PH<-as.numeric(NewData3$PH)
NewData3$Index<-as.numeric(NewData3$Index)

ggplot(NewData3,aes(x=PH,y=Index,col=variable))+
  facet_wrap(~variable,scales = "free")+
  geom_point() +
  geom_smooth(method=lm)+
  theme_bw()