DOI: 10.7910/dvn/cw20dz
title: Linear workflow for relating multispectral remote sensing imagery with on-farm plant measurements

description: &lt;b&gt;Linear workflow for relating multispectral remote sensing imagery with on-farm plant measurements&lt;/b&gt;
&lt;br&gt;&lt;br&gt;
1. &lt;i&gt;Automated spectral index computation in ArcPy&lt;/i&gt;
&lt;br&gt;
Spectral_Index_List.xlsx
&lt;br&gt;
Spectral_Index_Computation_ArcPy.txt
&lt;br&gt;
sUAS_Spectral_Index.jpg
&lt;br&gt;&lt;br&gt;

2. &lt;i&gt;Pearson Correlation Coefficient (PCC) plotting in R&lt;/i&gt;
&lt;br&gt;
PCC_Plotting_R.txt
&lt;br&gt;
PCC_ZStats_Median.jpg
&lt;br&gt;&lt;br&gt;

3. &lt;i&gt;Multiple linear regression (MLR) permutations in R & variance inflation factor (VIF) calculation&lt;/i&gt;
&lt;br&gt;
MLR_Permutations_R.txt
&lt;br&gt;
SPAD_Spectral_Index_Relation.jpg
&lt;br&gt;&lt;br&gt;

4. &lt;i&gt;LASSO (least absolute shrinkage and selection operator) regression in R&lt;/i&gt;
&lt;br&gt;
LASSO_Regression_R.txt
&lt;br&gt;&lt;br&gt;
---PART I-----------------------------------------------------------------------------------------------------------
&lt;br&gt;&lt;br&gt;
This dataverse entry was originally created to supplement a publication by Peter et al. (2020) that evaluated effective spatial resolutions for monitoring crop productivity from sUAS, commercial, and national government imaging platforms. The original entry included a workflow for automating spectral index computation in ArcPy and conducting multiple linear regressions (R) in R. 
&lt;br&gt;&lt;br&gt;
&lt;b&gt;Associated manuscript&lt;/b&gt;:
&lt;br&gt;
Peter, B. G., Messina, J. P., Carroll, J. W., Zhi, J., Chimonyo, V., Lin, S., and Snapp, S. S. 2020. Multi-spatial resolution satellite and sUAS imagery for precision agriculture on smallholder farms in Malawi. &lt;i&gt;Photogrammetric Engineering & Remote Sensing (PE&RS)&lt;/i&gt; 86(2):107–119. doi: &lt;a href="https://doi.org/10.14358/PERS.86.2.107"&gt;10.14358/PERS.86.2.107&lt;/a&gt;
&lt;br&gt;&lt;br&gt;
---PART II----------------------------------------------------------------------------------------------------------
&lt;br&gt;&lt;br&gt;
In follow-up research, the focus shifted toward evaluating a comprehensive collection of multispectral indices for monitoring agricultural productivity (Peter et al. &lt;i&gt;in review&lt;/i&gt;). The automated spectral index computation in ArcPy and the MLR regression code were updated enhance their utility. A variance inflation factor (VIF) calculation was added to the MLR code to further evaluate collinearity among predictor variables. Pearson Correlation Coefficient (PCC) plotting was also improved and code for LASSO (least absolute shrinkage and selection operator) regression in R was added. 
&lt;br&gt;&lt;br&gt; Three figures and a table were also added to this entry to supplement the Peter et al. (&lt;i&gt;in review&lt;/i&gt;) manuscript: 
&lt;li&gt;PCC_ZStats_Median.jpg, 
&lt;li&gt;SPAD_Spectral_Index_Relation.jpg
&lt;li&gt;sUAS_Spectral_Index.jpg, 
&lt;li&gt;and Spectral_Index_List.tab. 
&lt;br&gt;These are expanded versions of the figures available in the manuscript, which had to be shortened to meet length requirements. The table lists equations, abbreviations, and references of the 20 spectral indices used in the Peter et al. (&lt;i&gt;in review&lt;/i&gt;) manuscript.
&lt;br&gt;&lt;br&gt;
&lt;b&gt;Associated manuscript&lt;/b&gt;:
&lt;br&gt;
Peter, B. G., Messina, J. P., Carroll, J. W., and Chikowo, R. A case for green-based spectral indices: Plot-scale sUAS imagery relation to crop chlorophyll content at two maize farms in Malawi. &lt;b&gt;[Unpublished, submitted to &lt;i&gt;Remote Sensing Letters&lt;/i&gt; January 2021]&lt;/b&gt;.
&lt;br&gt;&lt;br&gt;
---CODE------------------------------------------------------------------------------------------------------------
&lt;br&gt;&lt;br&gt;
1. &lt;i&gt;Automated spectral index computation in ArcPy&lt;/i&gt;
&lt;br&gt;
&lt;pre&gt;
import arcpy, os
from arcpy.sa import *
arcpy.env.overwriteOutput = True

folder = "C:\\imagery\\folder\\"
name = "naming_prefix_"
region = "name"
R = Float(folder + "\\" + name + "R.tif")
NIR = Float(folder + "\\" + name + "NIR.tif")
G = Float(folder + "\\" + name + "G.tif")
shape = "C:\\zone\\shapes\\zone_shapes.shp"

output = arcpy.CreateFolder_management("C:\\imagery\\folder\\", "indices_"+region)
intFolder = arcpy.CreateFolder_management("C:\\imagery\\folder\\", "integer_"+region)
stats = arcpy.CreateFolder_management("C:\\imagery\\folder\\", "zstats_"+region)

outFile = G/(NIR+R+G)
outFile.save(str(output) + "\\" + name + "ngi_float.tif")

outFile = R/(NIR+R+G)
outFile.save(str(output) + "\\" + name + "nri_float.tif")

outFile = NIR/R
outFile.save(str(output) + "\\" + name + "rvi_float.tif")

outFile = NIR/G
outFile.save(str(output) + "\\" + name + "grvi_float.tif")

sq = SquareRoot((NIR/R)+1)
outFile = ((NIR/R)-1)/sq
outFile.save(str(output) + "\\" + name + "msr_float.tif")

outFile = NIR-R
outFile.save(str(output) + "\\" + name + "dvi_float.tif")

outFile = NIR-G
outFile.save(str(output) + "\\" + name + "gdvi_float.tif")

sq = SquareRoot(NIR+R)
outFile = (NIR-R)/sq
outFile.save(str(output) + "\\" + name + "rdvi_float.tif")

outFile = (NIR-R)/(NIR+R)
outFile.save(str(output) + "\\" + name + "ndvi_float.tif")

outFile = (NIR-G)/(NIR+G)
outFile.save(str(output) + "\\" + name + "gndvi_float.tif")

outFile = 1.5*((NIR-R)/(NIR+R+0.5))
outFile.save(str(output) + "\\" + name + "savi_float.tif")

outFile = 1.5*((NIR-G)/(NIR+G+0.5))
outFile.save(str(output) + "\\" + name + "gsavi_float.tif")

outFile = (NIR-R)/(NIR+R+0.16)
outFile.save(str(output) + "\\" + name + "osavi_float.tif")

outFile = (NIR-G)/(NIR+G+0.16)
outFile.save(str(output) + "\\" + name + "gosavi_float.tif")

sq = SquareRoot(((2*NIR)+1)*((2*NIR)+1)-(8*(NIR-R)))
outFile = (((2*NIR)+1)-sq)/2
outFile.save(str(output) + "\\" + name + "msavi2_float.tif")

outFile = ((0.2*NIR)-R)/((0.2*NIR)+R)
outFile.save(str(output) + "\\" + name + "wdrvi_float.tif")

outFile = 0.5*((120*(NIR-G))-(200*(R-G)))
outFile.save(str(output) + "\\" + name + "tvi_float.tif")

outFile = 1.2*((2.5*(NIR-R))-(1.3*(NIR-G)))
outFile.save(str(output) + "\\" + name + "mcari1_float.tif")

outFile = (NIR/G)*(R/G)
outFile.save(str(output) + "\\" + name + "cvi_float.tif")

outFile = ((NIR*NIR)-R)/((NIR*NIR)+R)
outFile.save(str(output) + "\\" + name + "nli_float.tif")

arcpy.env.workspace = str(output)

for raster in arcpy.ListRasters():
 name = os.path.basename(raster).rstrip(os.path.splitext(raster)[1])
 temp = arcpy.Times_3d(raster, 100000, str(intFolder) + "\\" + "temp.tif")
 outInt = Int(temp)
 outInt.save(str(intFolder) + "\\" + name + ".tif")

arcpy.env.workspace = str(intFolder)

for raster in arcpy.ListRasters():
 name = os.path.basename(raster).rstrip(os.path.splitext(raster)[1])
 table = str(stats) + "\\" + name + ".dbf"
 arcpy.gp.ZonalStatisticsAsTable(shape, "id_rt", raster, table, "DATA", "ALL")
 arcpy.TableToTable_conversion(table, str(stats), name + ".csv")
&lt;/pre&gt;
&lt;br&gt;&lt;br&gt;
2. &lt;i&gt;Pearson Correlation Coefficient (PCC) plotting in R&lt;/i&gt;
&lt;br&gt;
&lt;pre&gt;
library(GGally)

name &lt;- "naming_prefix"

cor &lt;- ggcorr(dataset, method = c("pairwise", "pearson"), digits = 2, name = "PCC", midpoint = 0, 
geom = "tile", label = TRUE, label_size = 5, label_round = 3, hjust = 0.9, layout.exp = 3, legend.size = 12)
png(filename=paste(name, "_PRS.jpg", sep=""), width=20, height=20, units="in", res=300)
pcor &lt;- cor(dataset, method = c("pearson"))
write.csv(pcor, file=paste(name, "_PRS.csv", sep=""))
cor
dev.off()
&lt;/pre&gt;
&lt;br&gt;&lt;br&gt;
3. &lt;i&gt;Multiple linear regression (MLR) permutations in R & variance inflation factor (VIF) calculation&lt;/i&gt;
&lt;br&gt;
&lt;pre&gt;
library(MuMIn)
library(car)

model &lt;- lm(SPAD ~ NGI + NRI + RVI + GRVI + MSR + DVI + GDVI + RDVI + NDVI + GNDVI + SAVI + GSAVI + 
OSAVI + GOSAVI + MSAVI2 + WDRVI + TVI + MCARI1 + CVI + NLI, dataset, na.action="na.fail")
var_count &lt;- 2^(length(coef(model))-1)

name &lt;- "naming_prefix"

combinations &lt;- pdredge(model, extra = list(
 "R^2", "AIC", "BIC", "*" = function(x) {
 p &lt;- 999
 vif &lt;- 999
 s &lt;- summary(x)
 f &lt;- s$fstatistic
 if( s$r.squared != 0 ) {
 p &lt;- pf(f[1],f[2],f[3],lower.tail=F)[[1]]
 }
 if( length(coef(x))-1 &gt; 1 ) {
 vif &lt;- max(vif(x))
 }
 c("F-stat" = f[[1]], 
 "p-value" = p, 
 "p-value_BC" = p*var_count, 
 "RSE" = s$sigma, 
 "adj-R^2" = s$adj.r.squared, 
 "maxVIF" = vif)
 })
)

write.csv(combinations, file=paste(name, "_MR_SPD.csv", sep=""))
&lt;/pre&gt;
&lt;br&gt;&lt;br&gt;
4. &lt;i&gt;LASSO (least absolute shrinkage and selection operator) regression in R&lt;/i&gt;
&lt;br&gt;
&lt;pre&gt;
library(readr)
library(rsample)
library(glmnet)

name &lt;- "naming_prefix"

datalist = list()

for (i in 1:10000) {
 d_split &lt;- initial_split(dataset, prop = 0.7)
 d_train &lt;- training(d_split)
 d_test &lt;- testing(d_split)

 train_x &lt;- model.matrix(SPAD ~ ., d_train)[, -1]
 train_y &lt;- d_train$SPAD

 test_x &lt;- model.matrix(SPAD ~ ., d_test)[, -1]
 test_y &lt;- d_test$SPAD

 cv_lasso &lt;- cv.glmnet(x = train_x, y = train_y, alpha = 1)

 mse_min &lt;- min(cv_lasso$cvm) 
 lambda_min &lt;- cv_lasso$lambda.min

 coefs &lt;- data.frame(as.matrix(coef(cv_lasso)))

 pred &lt;- predict(cv_lasso, s = lambda_min, test_x)
 mse &lt;- mean((test_y - pred)^2)

 merge &lt;- rbind(coefs, mse, sqrt(mse), mse_min, sqrt(mse_min), lambda_min)
 row.names(merge)[22] &lt;- "MSE_pred"
 row.names(merge)[23] &lt;- "RMSE_pred"
 row.names(merge)[24] &lt;- "MSE_min"
 row.names(merge)[25] &lt;- "RMSE_min"
 row.names(merge)[26] &lt;- "lambda_min"

 datalist[[i]] &lt;- merge
}

append = do.call(cbind, datalist)
write.csv(append, file=paste(name, "_LASSO.csv", sep=""))
&lt;/pre&gt;
&lt;br&gt;
&lt;i&gt;This content is made possible by the support of the American People provided to the Feed the Future Innovation Lab for Sustainable Intensification through the United States Agency for International Development (USAID). The contents are the sole responsibility of the authors and do not necessarily reflect the views of USAID or the United States Government. Program activities are funded by USAID under Cooperative Agreement No. AID-OAA-L-14-00006.&lt;/i&gt;