10. Distance between plants, objects length, and removing objects (plot, cloud, weed, etc.)

The function fieldObject can be used to take relative measurement from objects (e.g., area, "x" distance, "y" distance, number, extent, etc.) in the entire mosaic or per plot using the fieldShape file. Download example here: EX_Obj.jpg.

# Uploading file (EX_Obj.tif)
EX.Obj <- rast("EX_Obj.jpg")
plotRGB(EX.Obj)
EX.Obj <- aggregate(EX.Obj,4)

# Removing the background
EX.Obj<-imgLAB(EX.Obj)
EX.Obj.M<- fieldMask(mosaic = EX.Obj, index = "BGI",cropValue = 0.7,cropAbove = T)

# Taking measurements:
EX.Obj.D<-fieldCount(mosaic = EX.Obj.M$mask,
                     plot = T)

fieldView(EX.Obj,EX.Obj.D)
plotRGB(EX.Obj)
plot(EX.Obj.D$geometry, add=T, border="red")
plot(EX.Obj.D[1,], add=T, col="yellow")
EX.Obj.I<- fieldIndex(mosaic = EX.Obj,index = c("SI","BGI","BI"))
EX.Obj.Data<-fieldInfo_extra(mosaic = EX.Obj.I[[c("SI","BGI","BI")]], 
                             fieldShape = EX.Obj.D)
plot(EX.Obj.Data)

# Data visualization: 
library(reshape2)
Data.Obj1<-melt(EX.Obj.Data,measure.vars = c("area","perimeter","width","SI","BGI","BI"))

library(ggplot2)
ggplot(Data.Obj1, aes(x=value, fill=variable)) +
  geom_histogram(aes(y=..density..), colour="black")+
  geom_density(alpha=.2)+
  facet_wrap(~variable, scales = "free")

The function fieldCrop can be used to remove objects from the field image. For instance, the parameter remove=TRUE and nPoint should be used to select the object boundaries to be removed. Download EX_RemObj.tif).

# Uploading file (EX_RemObj.tif)
EX.RemObj <- rast("EX_RemObj.tif")

# Selecting the object boundaries to be removed (nPoint = 10)
EX.RemObj.Crop <- fieldCrop(mosaic = EX.RemObj, remove = T, nPoint = 10) # Selecting the plant in plot 13

# Removing the soil
EX.RemObj.RemSoil<- fieldMask(mosaic = EX.RemObj.Crop,index = "HUE")

# Building the plot shapefile (ncols = 8 and nrows = 4)
EX.RemObj.Shape<-fieldShape_render(mosaic = EX.RemObj.RemSoil,ncols = 8, nrows = 4)

# Building indice (NGRDI)
EX.RemObj.Indices<- fieldIndex(mosaic = EX.RemObj.RemSoil$newMosaic,index = c("NGRDI"))

# Extracting data (NGRDI)
EX.RemObj.Info<- fieldInfo_extra(mosaic = EX.RemObj.Indices$NGRDI,
                      fieldShape = EX.RemObj.Shape)

# Comparing plots values (the plant in plot 13 was removed and its value must be lower than plot 12 and 14)                      
EX.RemObj.Info[c(12,13,14),]

The function fieldDraw can be used to draw lines or polygons in the field image. This function allows to extract information of specific positions in the field (x, y, value). Also, this function can be used to evaluate distances between objects (for example: distance between plants in a line) or either objects length (for example: seed length, threes diameter, etc.). Let’s use the image above to evaluate the distance between potato plots and later to extract NGRDI values from a line with soil and vegetation to observe the profile in a distance plot.

# Uploading file (EX_RemObj.tif)
EX.Dist <- rast("EX_RemObj.tif")

# vegetation indices
EX.Dist.Ind <- fieldIndex(mosaic = EX.Dist,index = c("NGRDI"))

# Removing the soil
EX.Dist.RemSoil <- fieldMask(mosaic = EX.Dist.Ind)

# Evaluating distance between plants. Remember to press ESC when finished to draw the line.
EX.Dist.Draw <- fieldDraw(mosaic = EX.Dist.RemSoil$mask, 
                          dist = T,
                          value = 1, # Use the number to identify the object to be measured (1 for space and 0 for plants)
                          distSel = 0.4) 
EX.Dist.Draw$drawDist # Distance between plants

# Making plots                         
plot(EX.Dist.Ind$NGRDI)
points(EX.Dist.Draw$drawData$x,EX.Dist.Draw$drawData$y, col="red",pch=16,cex=0.7)
points(EX.Dist.Draw$drawSegments$x,EX.Dist.Draw$drawSegments$y, col="blue",pch=16,cex=0.7)
lines(EX.Dist.Draw$drawDist[1,c("x1","x2")],EX.Dist.Draw$drawDist[1,c("y1","y2")], col="green",lwd=5)

# Evaluating an specific layer profile (e.g. NGRDI)
EX.Dist.Draw.2 <- fieldDraw(mosaic = EX.Dist.Ind,
                            ndraw = 4, # Making 4 lines (press ESC to conclude each line)
                            lwd = 5)
EX.Data<-EX.Dist.Draw.2$Draw1$drawData
dev.off()
plot(x = EX.Data$x, y = EX.Data$NGRDI, type="l", col="red",lwd=2,xlab="Distance (m)", ylab="NGRDI")
abline(h=0.0,col="blue", lty=2, lwd=3)

# Making polygons and extracting data per cell
EX.Dist.Draw.3 <- fieldDraw(mosaic = EX.Dist.Ind,
                            line = F, # Making 2 polygons (press ESC to conclude each polygon)
                            ndraw = 2,
                            lwd = 5)
plotRGB(EX.Dist.RemSoil$newMosaic)
plot(EX.Dist.Draw.3$Draw1$drawObject, col="red",add=T)

Menu

### Image and resolution decrease 

RES_1<-rast("EX1_RGB_HighResolution.tif")
RES_2<-aggregate(RES_1, fact=2)
RES_4<-aggregate(RES_1, fact=4)

res(RES_1)
res(RES_2)
res(RES_4)

par(mfrow=c(1,3))
plotRGB(RES_1)
plotRGB(RES_2)
plotRGB(RES_4)

### Removing Soil 

system.time({RES_1_S <- fieldMask(RES_1,index="HUE")})
system.time({RES_2_S <- fieldMask(RES_2,index="HUE")})
system.time({RES_4_S <- fieldMask(RES_4,index="HUE")})

### Indices

system.time({RES_1_I <- fieldIndex(RES_1_S$newMosaic,index=c("BGI"))})
system.time({RES_2_I <- fieldIndex(RES_2_S$newMosaic,index=c("BGI"))})
system.time({RES_4_I <- fieldIndex(RES_4_S$newMosaic,index=c("BGI"))})

### Get Information (1 Band)

system.time({RES_1_Info <- fieldInfo_extra(RES_1_I$BGI,fieldShape = EX1.Shape)})
system.time({RES_2_Info <- fieldInfo_extra(RES_2_I$BGI,fieldShape = EX1.Shape)})
system.time({RES_4_Info <- fieldInfo_extra(RES_4_I$BGI,fieldShape = EX1.Shape)})

### Get Information (3 Bands)

system.time({RES_1_Info2 <- fieldInfo_extra(RES_1_I[[c(1,2,3)]],fieldShape = EX1.Shape)})
system.time({RES_2_Info2 <- fieldInfo_extra(RES_2_I[[c(1,2,3)]],fieldShape = EX1.Shape)})
system.time({RES_4_Info2 <- fieldInfo_extra(RES_4_I[[c(1,2,3)]],fieldShape = EX1.Shape)})

### Correlation

DataBGI <- data.frame(R1=RES_1_Info$BGI,
                    R2=RES_2_Info$BGI,
                    R4=RES_4_Info$BGI)
DataBlue <- data.frame(R1=RES_1_Info2$Blue,
                       R2=RES_2_Info2$Blue,
                       R4=RES_4_Info2$Blue)
DataGreen <- data.frame(R1=RES_1_Info2$Green,
                       R2=RES_2_Info2$Green,
                       R4=RES_4_Info2$Green)
DataRed <- data.frame(R1=RES_1_Info2$Red,
                       R2=RES_2_Info2$Red,
                       R4=RES_4_Info2$Red)
cor(DataBGI)
cor(DataBlue)
cor(DataGreen)
cor(DataRed)

# Uploading Flowering (EX2_RGB.tif) and Senescence (EX3_RGB.tif) files:
EX2 <- rast("EX2_RGB.tif")
EX3 <- rast("EX3_RGB.tif")

# Cropping the image using the previous shape from step 2:

EX2.Crop <- fieldCrop(mosaic = EX2,fieldShape = EX1.Crop, plot = T)
EX3.Crop <- fieldCrop(mosaic = EX3,fieldShape = EX1.Crop, plot = T)

# Rotating the image using the same theta from step 3:

EX2.Rotated<-fieldRotate(EX2.Crop,theta = 2.3, plot = T)
EX3.Rotated<-fieldRotate(EX3.Crop,theta = 2.3, plot = T)

# Removing the soil using index and mask from step 4:

EX2.RemSoil<-fieldMask(EX2.Rotated,index="HUE",cropValue=0,cropAbove=T,plot=T)
EX3.RS<-fieldMask(EX3.Rotated,index="HUE",cropValue=0,cropAbove=T,plot=T) # Removing soil at senescence stage
EX3.RemSoil<-fieldMask(EX3.RS$newMosaic,mask = EX2.RemSoil$mask ,cropValue=0,cropAbove=T,plot=T) # Removing weeds from senescence stage with flowering mask 

# Building indices

EX2.Indices <- fieldIndex(EX2.RemSoil$newMosaic,Red=1,Green=2,Blue=3,
                 index = c("NGRDI","BGI"), myIndex = c("(Red-Blue)/Green"))
EX3.Indices <- fieldIndex(EX3.RemSoil$newMosaic,Red=1,Green=2,Blue=3,
                 index = c("NGRDI","BGI"), myIndex = c("(Red-Blue)/Green"))

# Extracting data using the same fieldShape file from step 5:

EX2.Info<- fieldInfo_extra(mosaic = EX2.Indices$myIndex,fieldShape = EX1.Shape)
EX3.Info<- fieldInfo_extra(mosaic = EX3.Indices$myIndex,fieldShape = EX1.Shape)

Data.Cycle<-data.frame(EX1=EX1.Info$myIndex,
      EX2=EX2.Info$myIndex,
      EX3=EX3.Info$myIndex)

Data.Cycle

#####################
### Multispectral ###
#####################

# Uploading multispectral mosaic:
EX1.5b <- rast("EX1_5Band.tif")

# Cropping the image using the previous shape from step 2:
EX1.5b.Crop <- fieldCrop(mosaic = EX1.5b,fieldShape = EX1.Crop, plot = T)

# Rotating the image using the same theta from step 3:
EX1.5b.Rotated<-fieldRotate(EX1.5b.Crop,theta = 2.3, plot = T)

# Removing the soil using index and mask from step 4:
EX1.5b.RemSoil<-fieldMask(EX1.5b.Rotated,Red=1,Green=2,Blue=3,index="HUE",cropValue=0,cropAbove=T,plot=T)

# Building indices (NDVI and NDRE)
EX1.5b.Indices <- fieldIndex(EX1.5b.RemSoil$newMosaic,Red=1,Green=2,Blue=3,RedEdge=4,NIR=5,
                 index = c("NDVI","NDRE"))

# Extracting data using the same fieldShape file from step 5:
EX1.5b.Info<- fieldInfo_extra(mosaic = EX1.5b.Indices$NDVI,fieldShape = EX1.Shape)

#####################
### Hyperspectral ###
#####################

# Uploading hyperspectral file with 474 bands (EX_HYP.tif)
EX.HYP<-rast("EX_HYP.tif")

# Wavelengths (namesHYP.csv)
NamesHYP<-as.character(read.csv("namesHYP.csv")$NameHYP)

# Building RGB image 
R<-EX.HYP[[78]] # 651nm (Red)
G<-EX.HYP[[46]] # 549nm (Green)
B<-EX.HYP[[15]] # 450nm (Blue)
RGB<-rast(c(R,G,B))
plotRGB(RGB, stretch="lin")

# Removing soil using RGB (index NGRDI)
RGB.S<-fieldMask(RGB,index="NGRDI",cropValue = 0.0, cropAbove = F)

# Data frame with field information to make the Map
Data<-read.csv("DataHYP.csv")
Map<-fieldMap(fieldPlot = as.character(Data$Plot),fieldRow = as.character(Data$Range),fieldColumn = as.character(Data$Row),decreasing = T)

# Building plot shapefile using RGB as base
plotFile<-fieldShape_render(RGB.S,ncols = 14, nrows = 14, fieldMap = Map,fieldData = Data, PlotID = "Plot")

# Removing soil using the RGB mask
EX.HYP.S<-fieldMask(EX.HYP,mask = RGB.S$mask, plot = F)

# Extracting data (474 bands)
EX.HYP.I<-fieldInfo_extra(EX.HYP.S$newMosaic,fieldShape = plotFile)

# Saving the new csv with hyperspectral information per plot
DataHYP<-EX.HYP.I[,-dim(EX.HYP.I)[2]]
colnames(DataHYP)<-c(colnames(DataHYP)[1:9],NamesHYP)
write.csv(DataHYP,"DataHypNew.csv",col.names = T,row.names = F)

###############
### Graphic ###
###############

dev.off()
DataHYP1<-EX.HYP.I$plotValue[,-1]

plot(x=as.numeric(NamesHYP),y=as.numeric(DataHYP1[1,]),type = "l",xlab = "Wavelength (nm)",ylab = "Reflectance", col="black",lwd=2,cex.lab=1.2)
for(i in 2:dim(DataHYP1)[2]){
  lines(x=as.numeric(NamesHYP),y=as.numeric(DataHYP1[i,]),type = "l",col=i,lwd=2)
}
abline(v=445,col="blue",lwd=2,lty=2)
abline(v=545,col="green",lwd=2,lty=2)
abline(v=650,col="red",lwd=2,lty=2)
abline(v=720,col="red",lwd=2,lty=3)
abline(v=840,col="red",lwd=2,lty=4)
legend(list(x = 2000,y = 0.5),c("Blue (445nm)","Green (545nm)","Red (650nm)","RedEdge (720nm)","NIR (840nm)"),
       col =c("blue","green","red","red","red"),lty=c(2,2,2,3,4),box.lty=0)

# Uploading file (EX_polygonShape.tif)
EX.polygon<-rast("EX_polygonShape.tif")
plotRGB(EX.polygon, r = 1, g = 2, b = 3)

# Removing soil
EX.polygon.RemSoil<- fieldMask(mosaic = EX.polygon)

# Data frame with polygons information
polygonData<-data.frame(ID=c("Polygon1","Polygon2","Polygon3"),
                        FlowerColor=c("white","white","white"),
                        FlowerPercent=c(20,40,50),
                        LeafColor=c("dark","light","dark"))
polygonData

# Building plot shapefile with 3 polygons (select 4 points around the polygon area)
EX.polygon.Shape<-fieldPolygon(mosaic = EX.polygon.RemSoil,
                               nPolygon = 3,nPoint = 4,ID = "ID",
                               polygonData = polygonData,cropPolygon = T,
                               polygonID = c("Polygon1","Polygon2","Polygon3"))
plotRGB(EX.polygon.Shape$cropField)

# Building indice (NGRDI and BGI)
EX.polygon.Indices<- fieldIndex(mosaic = EX.polygon.RemSoil$newMosaic, Red = 1, Green = 2, Blue = 3, 
                             index = c("NGRDI","BGI"))

# Extracting data (NGRDI and BGI)
EX.polygon.Info<- fieldInfo_extra(mosaic = EX.polygon.Indices[[c("NGRDI","BGI")]],
                   fieldShape = EX.polygon.Shape)
EX.polygon.Info

# Making graphics (BGI)
fieldPlot(fieldShape=EX.polygon.Info,
          fieldAttribute="BGI",
          mosaic=EX.polygon, color=c("red","blue"), alpha = 0.5)

### Interpolating colors: c("white","black")
fieldPlot(fieldShape=EX1.Info$fieldShape,fieldAttribute="Yield", mosaic=EX1.Indices, color=c("white","black"), alpha = 0.5)

### Interpolating colors: c("red","blue")
fieldPlot(fieldShape=EX1.Info$fieldShape,fieldAttribute="myIndex", mosaic=EX1.Indices, color=c("red","blue"), alpha = 0.5)

### Images (single and multi layers)
writeRaster(EX1.Indices, filename="EX1.Indices.tif", overwrite=TRUE)
# EX1.Indices.2 <- rast("EX1.Indices.tif") # Reading the saved image.

### FieldShape file
st_write(EX1.Info$fieldShape, ".", "EX1.fieldShape", driver="ESRI Shapefile")
# EX1.fieldShape.2 <- st_read("EX1.fieldShape.shp") # Reading the saved shapefile option 01.

### CSV file (table)
write.csv(EX1.Info,file = "EX1.Info.csv",col.names = T,row.names = F)
# Data.EX1.Info<-read.csv("EX1.Info.csv",header = T,check.names = F) # Reading the saved data table.

# Uploading file (odm_orthophoto.tif):
EX.ODM<-rast("odm_orthophoto.tif")
plotRGB(EX.ODM, r = 1, g = 2, b = 3)

# Cropping the image to select only one trial (Selecting the same trial as EX.2 from step.13):
EX.ODM.Crop <- fieldCrop(mosaic = EX.ODM)

# Rotating the image using the same theta from step 3:
EX.ODM.Rotated<-fieldRotate(EX.ODM.Crop,theta = 2.3)

# Removing soil
EX.ODM.RemSoil<- fieldMask(mosaic = EX.ODM.Rotated)

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

# Extracting data using the same fieldShape file from step 5:
EX.ODM.Info<- fieldInfo_extra(mosaic = EX.ODM.Indices$myIndex,fieldShape = EX1.Shape)

EX.ODM.Info$myIndex

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

# Vegetation indices
index<- c("BGI","VARI","SCI")

############
### Loop ###
############

system.time({ # system.time: used to compare the processing time using loop and parallel
EX.Table.Loop<-NULL
for(i in 1:length(pics)){
  EX.L1<-rast(paste("./images/",pics[i],sep = ""))
  plotRGB(EX.L1)
  EX.L.Shape<-fieldPolygon(mosaic=EX.L1, extent=T, plot=F) # extent=T (The whole image area will be the shapefile)
  EX.L2<-fieldMask(mosaic=EX.L1, index="BGI", cropValue=0.8, cropAbove=T, plot=F) # Select one index to identify leaves and remove the background
  EX.L3<-fieldMask(mosaic=EX.L2$newMosaic, index="VARI", cropValue=0.1, cropAbove=T, plot=F) # Select one index to identify demaged area in the leaves  
  EX.L4<-fieldIndex(mosaic=EX.L2$newMosaic, index=index, plot=F) # Indices
  EX.L5<-rast(EX.L3$mask, EX.L4[[index]]) # Making a new rast raster with new layers (demage area and indices)
  EX.L.Info<- fieldInfo_extra(mosaic=EX.L5, fieldShape=EX.L.Shape, projection=F) # projection=F (Ignore projection. Normally used only with remote sensing images)
  plot(EX.L5,col = grey(1:100/100))
  EX.Table.Loop<-rbind(EX.Table.Loop, EX.L.Info$plotValue) # Combine information from all images in one table
}})
rownames(EX.Table.Loop)<-pics
EX.Table.Loop

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

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

# Number of cores
n.core<-detectCores()-1

# Starting parallel
cl <- makeCluster(n.core, output = "")
registerDoParallel(cl)
system.time({
EX.Table.Parallel <- foreach(i = 1:length(pics), .packages = c("raster","FIELDimageR"), 
                     .combine = rbind) %dopar% {
                       EX.L1<-rast(paste("./images/",pics[i],sep = ""))
                       EX.L.Shape<-fieldPolygon(mosaic=EX.L1, extent=T, plot=F) # extent=T (The whole image area will be the shapefile)
                       EX.L2<-fieldMask(mosaic=EX.L1, index="BGI", cropValue=0.8, cropAbove=T, plot=F) # Select one index to identify leaves and remove the background
                       EX.L3<-fieldMask(mosaic=EX.L2$newMosaic, index="VARI", cropValue=0.1, cropAbove=T, plot=F) # Select one index to identify demaged area in the leaves  
                       EX.L4<-fieldIndex(mosaic=EX.L2$newMosaic, index=index, plot=F) # Indices
                       EX.L5<-rast(EX.L3$mask, EX.L4[[index]]) # Making a new rast raster with new layers (demage area and indices)
                       EX.L.Info<- fieldInfo_extra(mosaic=EX.L5, fieldShape=EX.L.Shape, projection=F) # projection=F (Ignore projection. Normally used only with remote sensing images)
                       EX.L.Info$plotValue # Combine information from all images in one table
                     }})
stopCluster(cl)
rownames(EX.Table.Parallel)<-pics
EX.Table.Parallel 

# The rasterOptions() allows you to customize your R session (raster package): 
rasterOptions()
rasterOptions(chunksize = 1e+09)
rasterOptions(maxmemory = 1e+09)

new_mosaic<-aggregate(previous_mosaic, fact=4)

rm(previous_mosaic)

ShapeFile <- rgdal::readOGR("Other_Software_ShapeFile.shp") # Option 01
ShapeFile <- terra::vect("Other_Software_ShapeFile.shp") # Option 02

ShapeFile <- rbind(ShapeFile1, ShapeFile2, ShapeFile3, ...)

# Example in macOS (type the following code before running FIELDimageR functions):
options(device = "quartz")