FIELDimageR pipeline:
Image analysis applied to plant breeding

Introduction
- Plant Breeding
- FIELDimageR
FIELDimageR pipeline
- 1) Example 01
- 2) Example 02
Reference

Introduction

Plant Breeding

Images can be used in plant breeding to draw inference about many traits:

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

FIELDimageR

FIELDimageR is a R package to analyze images from plant breeding programs, and allows to:

Crop the image
Remove soil effect
Build vegetation indices
Rotate the image
Build the plot shapefile
Extract information for each plot
Evaluate stand count, canopy percentage, and plant height

FIELDimageR pipeline

1) Example 01

Remote sensing (Potato Breeding)

Jeffrey Endelman and Filipe Matias (UW-Madison)

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("lme4")
# install.packages("plyr")
# install.packages("DescTools")
# devtools::install_github("filipematias23/FIELDimageR")

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

Data (R/FIELDimageR)

Potato breeding: 140 plots (3x15 feet)
Population: 97 genotypes
Six UAV flights: 30 (soil base), 40, 55, 75, 90, and 100 days after planting (DAP)

Download EX1: RGB orthomosaics.tif MOSAIC and dsm.tif DSM


# EX1 data files: 
MOSAIC<-list.files("./MOSAIC/")
DSM<-list.files("./DSM/")

Basic steps

# Uploading an example mosaic
Test <- stack(paste("./MOSAIC/",MOSAIC[4],sep = ""))
plotRGB(Test)

# Cropping and reducing size 
Test.Crop <- fieldCrop(mosaic = Test)

# Rotating image with theta 2.3:
Test.Rotate<-fieldRotate(Test.Crop, theta = 2.3)
#Test.Rotate<-fieldRotate(Test.Crop, clockwise = F)

# Removing soil (example using HUE):
Test.RemSoil<-fieldMask(Test.Rotate)

Building the plot shapefile

Dataset - Field notes - phenotype (Download: ‘EX1_Data.csv’)

# Reading FieldData.csv
Data<-read.csv(file = "EX1_Data.csv",header = T)
Data

# Making the field Map
Map<-fieldMap(fieldPlot = Data$Plot,
              fieldColumn = Data$Column,
              fieldRow = Data$Row,
              decreasing = T)
Map

EX1 field design (MAP) should start from bottom to top and left to the right. In other words, plot 1 should be where the blue arrow is showing in the image below.

# Rotating the MAP to correct plots position:
rotate <- function(x) t(apply(x, 2, rev))
Map<-rotate(rotate(Map))
Map

# Building the plot shapefile (ncols = 14 and nrows = 10)
# x11()
plotShape<-fieldShape(mosaic = Test.RemSoil$newMosaic, 
                      ncols = 14, 
                      nrows = 10, 
                      fieldData = Data, 
                      ID = "Plot", 
                      fieldMap = Map)

# Building indices ("NGRDI","GLI")
Test.Indices<- fieldIndex(mosaic = Test.RemSoil$newMosaic, 
                         Red = 1, Green = 2, Blue = 3, 
                         index = c("NGRDI","GLI"), 
                         myIndex = c("(Red-Blue)/Green"))

plot(Test.Indices$myIndex)
plot(plotShape$fieldShape,add=T)

Extracting Data

Test.Info<- fieldInfo(mosaic = Test.Indices[[c("NGRDI","GLI","myIndex")]],
                      fieldShape = plotShape$fieldShape,
                      buffer = -0.10,
                      n.core = 3)
Test.Info$fieldShape@data

Estimation of plant height using digital surface model (DSM).

# Uploading files from soil base and vegetative growth:
DSM0 <- stack(paste("./DSM/",DSM[1],sep = ""))
DSM1 <- stack(paste("./DSM/",DSM[4],sep = ""))
plot(DSM1)

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

# Canopy Height Model (CHM):
DSM0.R <- resample(DSM0.Crop, DSM1.Crop)
CHM <- DSM1.Crop-DSM0.R
plot(CHM)

# Rotating the image using the same theta from RGB mosaic:
CHM.Rotate<-fieldRotate(CHM, theta = 2.3)

# Removing the soil using the mask from RGB mosaic:
CHM.RemSoil <- fieldMask(CHM.Rotate, mask = Test.RemSoil$mask)

# Extracting the estimate plant height average (EPH):
EPH <- fieldInfo(CHM.RemSoil$newMosaic, 
                 fieldShape = Test.Info$fieldShape, 
                 fun = "mean") #fun="quantile"
EPH$plotValue

Evaluating all mosaics in a loop

DataTotal<-NULL
for(i in 2:length(MOSAIC)){
  EX1 <- stack(paste("./MOSAIC/",MOSAIC[i],sep = ""))
  EX1.Crop <- fieldCrop(mosaic = EX1,
                    fieldShape = Test.Crop,
                    plot = F)
  EX1.Rotate<-fieldRotate(EX1.Crop,
                    theta = 2.3,
                    plot = F)
  EX1.RemSoil<-fieldMask(EX1.Rotate, plot = F)
  EX1.Indices<- fieldIndex(mosaic = EX1.RemSoil$newMosaic, 
                            Red = 1, Green = 2, Blue = 3, 
                            index = c("NGRDI","BGI", "GLI","VARI"), 
                            myIndex = c("(Red-Blue)/Green"),
                            plot = F)
  EX1.Info<- fieldInfo(mosaic = EX1.Indices[[c("NGRDI","BGI", "GLI","VARI","myIndex")]],
                        fieldShape = plotShape$fieldShape, 
                        n.core = 3)
  DSM0 <- stack(paste("./DSM/",DSM[1],sep = ""))
  DSM1 <- stack(paste("./DSM/",DSM[i],sep = ""))
  DSM0.Crop <- fieldCrop(mosaic = DSM0,fieldShape = Test.Crop,plot = F)
  DSM1.Crop <- fieldCrop(mosaic = DSM1,fieldShape = Test.Crop,plot = F)
  DSM0.R <- resample(DSM0.Crop, DSM1.Crop)
  CHM <- DSM1.Crop-DSM0.R
  CHM.Rotate<-fieldRotate(CHM, theta = 2.3,plot = F)
  CHM.RemSoil <- fieldMask(CHM.Rotate, mask = EX1.RemSoil$mask,plot = F)
  EPH <- fieldInfo(CHM.RemSoil$newMosaic, 
                   fieldShape = EX1.Info$fieldShape, 
                   fun = "mean")
  DataTotal<-rbind(DataTotal,
                   data.frame(DAP=as.character(do.call(c,strsplit(MOSAIC[i],split = "_"))[2]),
                              EPH$fieldShape@data))
  # Making map plots
  fieldPlot(fieldShape=EPH$fieldShape,
            fieldAttribute="NGRDI", 
            mosaic=EX1.Rotate, 
            color=c("red","green"), 
            alpha = 0.5,
            round = 2)
  print(paste("### Completed: ", "Mosaic_",i," ###",sep=""))
  }


# Organizing the data table:
colnames(DataTotal)<-c(colnames(DataTotal)[-c(dim(DataTotal)[2])],"EPH") # layer=EPH
DataTotal<-DataTotal[,!colnames(DataTotal)%in%c("ID","ID.1","PlotName")] # Removing column 12 ("ID.1")
DataTotal

# Saving extracted data "DataTotal.csv":
write.csv(DataTotal,"DataTotal.csv",row.names = F,col.names = T)

Graphics

DataTotal$Name<-as.factor(as.character(DataTotal$Name))
DataTotal$Row<-as.factor(as.character(DataTotal$Row))
DataTotal$Column<-as.factor(as.character(DataTotal$Column))
DataTotal$DAP<-as.numeric(as.character(DataTotal$DAP))
DataTotal$NGRDI<-as.numeric(as.character(DataTotal$NGRDI))
DataTotal$BGI<-as.numeric(as.character(DataTotal$BGI))
DataTotal$EPH<-as.numeric(as.character(DataTotal$EPH))

ggplot(DataTotal, aes(x = NGRDI,fill=as.factor(DAP))) +
  geom_density(alpha=.5,position = 'identity') +
  facet_wrap(~DAP,ncol = 1)+
  scale_fill_grey(start=1, end=0)+
  labs(y="#genotypes",x="NGRDI", fill="DAP") +
  theme_bw()

Heritability

# NGRDI  (mixed model):

DAP<-unique(DataTotal$DAP)

H2<-NULL
for(h in 1:length(DAP)){
  mod<-lmer(NGRDI~Row+Column+(1|Name),DataTotal[as.character(DataTotal$DAP)==DAP[h],])
  H2<-c(H2,as.data.frame(VarCorr(mod))$vcov[1]/sum(as.data.frame(VarCorr(mod))$vcov))
}

NewData<-data.frame(DAP,H2)

ggplot(NewData,aes(x=as.factor(DAP),y=H2,fill=as.factor(DAP)))+
  geom_bar(stat="identity")+
  scale_fill_grey(start=0.8, end=0.2)+
  labs(title="NGRDI",x="Days After Planting (DAP)", fill="")+
  geom_text(aes(label=round(H2,2)), vjust=1.6, color="white", size=6)+
  theme_bw()+
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

# Estimated plant height (mixed model):

H2<-NULL
for(h in 1:length(DAP)){
  mod<-lmer(EPH~Row+Column+(1|Name),DataTotal[as.character(DataTotal$DAP)==DAP[h],])
  H2<-c(H2,as.data.frame(VarCorr(mod))$vcov[1]/sum(as.data.frame(VarCorr(mod))$vcov))
}

NewData<-data.frame(DAP,H2)

ggplot(NewData,aes(x=as.factor(DAP),y=H2,fill=as.factor(DAP)))+
  geom_bar(stat="identity")+
  scale_fill_grey(start=0.8, end=0.2)+
  labs(title="EPH",x="Days After Planting (DAP)", fill="")+
  geom_text(aes(label=round(H2,2)), vjust=1.6, color="white", size=6)+
  theme_bw()+
  theme(legend.position = "none",
        plot.title = element_text(hjust = 0.5))

Area under the curve (AUC)

DataTotal1<-DataTotal[as.character(DataTotal$Name)%in%c("G43","G44","G45"),]

ggplot(data=DataTotal1, aes(x=as.numeric(DAP), y= NGRDI, col= Name, group=Name)) +
  geom_point(size=6)+
  geom_line(size=1.2) +
  scale_color_grey(start=0.8, end=0.2)+
  labs(x="Days After Planting (DAP)", fill="", col="")+
  theme_linedraw()

Trait<-c("NGRDI") # c("GLI","EPH")
Plot<-as.character(unique(DataTotal$Plot))

DataAUC<-NULL
for(a1 in 1:length(Plot)){
  D1<-DataTotal[as.character(DataTotal$Plot)==Plot[a1],]
  x1<-c(0,as.numeric(D1$DAP))
  y1<-c(0,as.numeric(D1[,Trait]))
  DataAUC <- rbind(DataAUC,
                   c(Name=unique(as.character(D1$Name)),
                     Row=unique(D1$Row),
                     Column=unique(D1$Column),
                     NGRDI_AUC=AUC(x = x1[!is.na(y1)], y = y1[!is.na(y1)])))}
                     
DataAUC<-as.data.frame(DataAUC)
DataAUC$NGRDI_AUC<-as.numeric(as.character(DataAUC$NGRDI_AUC))
DataAUC$Name<-as.factor(DataAUC$Name)
DataAUC$Row<-as.factor(DataAUC$Row)
DataAUC$Column<-as.factor(DataAUC$Column)
DataAUC

AUC Heritability

# Basic mixed model:
mod<-lmer(NGRDI_AUC~Row+Column+(1|Name),DataAUC)
H2<-as.data.frame(VarCorr(mod))$vcov[1]/sum(as.data.frame(VarCorr(mod))$vcov)

ggplot(DataAUC, aes(x = NGRDI_AUC)) +
  geom_histogram(aes(y=..density..), colour="black", fill="white")+
  geom_density(alpha=.5,position = 'identity', fill="#FF6666") +
  labs(title=paste("NGRDI_AUC (H2=",round(H2*100,2),"%)",sep=""), y="#genotypes",x="Area under the curve (AUC)") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5))

2) Example 02

Evaluating disease damage (Tomato Breeding)

Fernando Piotto and Jéssica Nogueira (ESALQ/USP)

Donwload: EX2.zip
Pictures of 10 genotypes with Xanthomonas perforans

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

# Vegetation indices
index<- c("BI","GLI","NGRDI")

#Choose one image to prepare the pipeline
EX.L1<-stack(paste("./EX2/",pics[1],sep = ""))
plotRGB(EX.L1)

# Shapefile with extent=T (The whole image area will be the shapefile)
EX.L.Shape<-fieldPolygon(mosaic=EX.L1, extent=T)

# Select one index to identify leaves and remove the background
EX.L2<-fieldMask(mosaic=EX.L1, index="BGI", cropValue=0.8, cropAbove=T)

plotRGB(EX.L2$newMosaic)

# Select one index to identify damaged area in the leaves  
EX.L3<-fieldMask(mosaic=EX.L2$newMosaic, index="VARI", cropValue=0.1, cropAbove=T)

plotRGB(EX.L3$newMosaic)

# Indices
EX.L4<-fieldIndex(mosaic=EX.L2$newMosaic, index=index)

# Making a new stack raster with new layers (damage area and indices)
EX.L5<-stack(EX.L4[[index]],(1-EX.L2$mask),(1-EX.L3$mask)) 
names(EX.L5)<-c(index,"Area","Damage")
plot(EX.L5)

# projection=F (Ignore projection. Normally used only with remote sensing images)
EX.L.Info<- fieldInfo(mosaic=EX.L5, 
                      fieldShape=EX.L.Shape$fieldShape, 
                      projection=F)

# Combine information from all images in one table
EX.L.Info$plotValue

* Parallel


# Installing 

# install.packages("foreach")
# install.packages("parallel")
# install.packages("doParallel")

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

# Number of cores
n.core<-3

# Starting parallel
cl <- makeCluster(n.core, output = "")
registerDoParallel(cl)
system.time({
  EX.Table <- foreach(i = 1:length(pics), .packages = c("raster","FIELDimageR"), 
                               .combine = rbind) %dopar% {
                                 EX.L1<-stack(paste("./EX2/",pics[i],sep = ""))
                                 EX.L.Shape<-fieldPolygon(mosaic=EX.L1, extent=T, plot=F)
                                 EX.L2<-fieldMask(mosaic=EX.L1, index="BGI", cropValue=0.8, cropAbove=T, plot=F) 
                                 EX.L3<-fieldMask(mosaic=EX.L2$newMosaic, index="VARI", cropValue=0.1, cropAbove=T, plot=F)   
                                 EX.L4<-fieldIndex(mosaic=EX.L2$newMosaic, index=index, plot=F)
                                 EX.L5<-stack(EX.L4[[index]],(1-EX.L2$mask),(1-EX.L3$mask)) 
                                 names(EX.L5)<-c(index,"Area","Damage")
                                 EX.L.Info<- fieldInfo(mosaic=EX.L5, fieldShape=EX.L.Shape$fieldShape, projection=F) 
                                 EX.L.Info$plotValue 
                               }})

EX.Table.2<-data.frame(Genotype=do.call(rbind,strsplit(pics,split = ".jpeg")),EX.Table[,-1])
EX.Table.2

Dataset: Visual Score
Horsfall-Barratt scale
4 Evaluators
(Download: ‘EX2_Data.txt’)

# Field data
DataEX2<-read.table("EX2_Data.txt",header = T)
DataEX2$Score_HB<-as.numeric(as.character(DataEX2$Score_HB))
DataEX2$Person<-as.factor(DataEX2$Person)
DataEX2$Genotype<-as.factor(DataEX2$Genotype)
DataEX2

ggplot(DataEX2,aes(x=Genotype,y=Score_HB,fill=Person))+
  facet_wrap(~Genotype,scales = "free_x",nrow = 2)+
  geom_bar(stat="identity",position=position_dodge())+
  scale_fill_grey()+
  theme_bw()+
  theme(axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

# Regression
DataEX2.1<-ddply(DataEX2,~Genotype,summarise,Score_HB=mean(Score_HB))
DataEX2.2<-merge(DataEX2.1,EX.Table.2,by="Genotype")
DataEX2.2

DataEX2.3<-melt(DataEX2.2,
               value.name = "Index",
               measure.vars = c("Damage","BI","GLI","NGRDI"))
DataEX2.3$Score_HB<-as.numeric(DataEX2.3$Score_HB)
DataEX2.3$Index<-as.numeric(DataEX2.3$Index)
DataEX2.3$variable<-as.factor(DataEX2.3$variable)
DataEX2.3

ggplot(DataEX2.3,aes(x=Index,y=Score_HB,col=variable))+
  facet_wrap(~variable,scales = "free_x")+
  geom_point() +
  geom_smooth(method=lm)+
  theme_bw()

Reference

Matias, F.I. (2019). FIELDimageR Pipeline (tutorial). https://github.com/filipematias23/FIELDimageR

Matias, F.I.; Caraza-Harter, M.; Endelman, J.B. (2020). FIELDimageR: A R Package to Analyze Orthomosaic Images from Agricultural Field Trials. The Plant Phenome Journal. DOI:10.1002/ppj2.20005

FIELDimageR pipeline: Image analysis applied to plant breeding

FIELDimageR pipeline: Image analysis applied to plant breeding

Introduction

Plant Breeding

FIELDimageR

FIELDimageR pipeline

1) Example 01

Remote sensing (Potato Breeding)

Required packages

Data (R/FIELDimageR)

Basic steps

Building the plot shapefile

Extracting Data

Estimation of plant height using digital surface model (DSM).

Evaluating all mosaics in a loop

Area under the curve (AUC)

2) Example 02

Evaluating disease damage (Tomato Breeding)

* Parallel

Reference

FIELDimageR pipeline:
Image analysis applied to plant breeding

FIELDimageR pipeline:
Image analysis applied to plant breeding