## Date: 9-25-2019 ## Personal: Jagdeep Singh Sidhu (jagdeeproots@gmail.com) ## Purpose: This code: ## 1. Calculates voronoi and Rustic values based on 6 coring locations using the SoilCoreTool developed in "An analysis of soil coring strategies to estimate root depth in maize (Zea mays) and common bean (Phaseolus vulgaris)" ## which is available on Gitlab (https://gitlab.com/ericnord/soilCoreTools). ## ## Data needed: The data used here is comes for field season 2015 and 2016. File structure: Make two main folders named "2015" and "2016" and in each folder add corresponding ## field data files namely - "2015_core_length_JB_format11june18_jgdp_extraction_1172020.csv" and "2016__core_length_JB_format11june18_jgdp_extraction_1172020.csv". ####################################################### ## Needed packages can be loaded using load functions## load_funcs <- function (){ ## Source files below can be found on Gitlab (https://gitlab.com/ericnord/soilCoreTools) source("D:/DES/Penn State/roots lab/Coring_study_Jimmy/gitlab/soilCoreTools-master/R/CheckCoreData.R") source("D:/DES/Penn State/roots lab/Coring_study_Jimmy/gitlab/soilCoreTools-master/R/CoreLayout.R") source("D:/DES/Penn State/roots lab/Coring_study_Jimmy/gitlab/soilCoreTools-master/R/PlotCores.R") source("D:/DES/Penn State/roots lab/Coring_study_Jimmy/gitlab/soilCoreTools-master/R/CoreVoronoiWeights.R") source("D:/DES/Penn State/roots lab/Coring_study_Jimmy/gitlab/soilCoreTools-master/R/VoronoiWeightedReferences.R") source("D:/DES/Penn State/roots lab/Coring_study_Jimmy/gitlab/soilCoreTools-master/R/CompareCoreCombinations.R") source("D:/DES/Penn State/roots lab/Coring_study_Jimmy/gitlab/soilCoreTools-master/R/sad.R") source("D:/DES/Penn State/roots lab/Coring_study_Jimmy/gitlab/soilCoreTools-master/R/DX.R") source("D:/DES/Penn State/roots lab/Coring_study_Jimmy/gitlab/soilCoreTools-master/R/CoreVisualizationTools.R") #install.packages("tripack") #install.packages("tidyverse") #install.packages("dplyr") #install.packages("plotrix") #install.packages("Rmisc") #install.packages("equivalence) library(dplyr) library(plotrix) library("tripack") library("tidyverse") library(ggplot2) library(reshape2) library(Rmisc) ## for summarySE library(ggrepel) ## for avoiding overlapping labels library(equivalence) ## for TOST } ############################################################################# #################Editable variables - list ################################## seasons <- c("2015","2016") crops <- c("Bean","Maize") ###############one time inititation of the combined results dataframe TOST_Combined_meth = data.frame("Depth", "Location", "comparedto", "P-value", "decision","Type", stringsAsFactors=FALSE) ### dataframe for combined results RVW_combined <- data.frame("Depth", "Length.cm", "Rep", "method","type", "crop", "Lengthcmpercm3", stringsAsFactors=FALSE) names(RVW_combined) <- c("Depth", "Length.cm", "Rep", "method","type", "crop","Lengthcmpercm3" ) RVW_combined <- RVW_combined[-1,] voronoi_field_combined <- data.frame() ############################################################################################### ############################# 1. SET directory function ####################################### setdir <- function(){ directory <- paste("D:/DES/Penn State/roots lab/Coring_study_Jimmy/gitlab/",year, sep ="") #Set Working Directory setwd(paste(directory)) getwd() message(paste("/working directory - ",directory)) ### read file for voronoi and rustic averages file_name = read.csv(paste(year,"_core_length_JB_format11june18_jgdp_extraction_1172020.csv",sep="")) %>% mutate(Depth = (Depth) * (-10)) %>% dplyr::select(Location, Rep, Depth, crop, core.volume,Length.cm.) file_name$RLD = file_name$Length.cm./file_name$core.volume file_name <<- file_name } ############################################################################################################### ##################### 2. Core layout around the focal plant ################################################### layout <- function(){ if (Crop == "Bean" & year == "2016"){layout.simroot <- coreLayout(px=76.2,py=23,cores=c(1,2,3,5,6),x.locs="c(0.25*py,0,0.5*py,0.5*px,0)", y.locs="c(0,0.5*py,0,0.5*py,0.25*py)",ids=c("1","2","3","5","6"), type=c("e","c","e","c","e"))}else{ layout.simroot <- coreLayout(px=76.2,py=23,cores=1:6,x.locs="c(0.25*py,0,0.5*py,0.5*px,0.5*px,0)", y.locs="c(0,0.5*py,0,0,0.5*py,0.25*py)",ids=c("1","2","3","4","5","6"), type=c("e","c","e","c","c","e"))} layout.maize = layout.simroot plotCores(layout.simroot,size="whole") layout.maize <<- layout.maize Simroot.weights<-coreVoronoiWeights(layout.simroot,plot.cores="whole") Simroot.weights <<- Simroot.weights } # only some core locations #coreVoronoiWeights(layout.beanbean,plot.cores="whole",which=c(1,3,4,6)) # only some ################################################################################################################## ##### 4. Calculate Voronoi values ############################################################################# voronoi_by_rep = data.frame("Depth", "RLD", "Rep", "year", "crop", stringsAsFactors=FALSE) ### dataframe for saving results vor_by_rep <- function(){ colnames(voronoi_by_rep) <- c("Depth", "RLD", "Rep", "year", "crop") voronoi_by_rep <- voronoi_by_rep[-1,] i = 1 ## variable for increasing the row number in results dataframe r =1 if (year == "2015"){reps = 1:8} else {reps = 1:5} for (r in reps) { voronoi_temp <- voronoiWeightedReferenceValues(core.weights=Simroot.weights, core.data=subset(file_name,Rep == r & crop == Crop), data.cols=7) rep_year_crop <- data.frame("Rep" = r, "year" = year, "crop" = Crop) # creating a data frame for rep and type rep_year_crop1 <- rep_year_crop %>% slice(rep(1:n(), each = 6)) voronoi_temp <- cbind(voronoi_temp, rep_year_crop1) names(voronoi_temp)[2] <- "RLD" voronoi_by_rep <- rbind(voronoi_temp, voronoi_by_rep) i = i+6 r = r+1 voronoi_by_rep <<- voronoi_by_rep } voronoi_by_rep$Location <- "voronoi" #write.csv(voronoi_by_rep, paste(year, Crop,"voronoibyrep.csv", sep ="")) rep_year_crop1 <<- rep_year_crop1 voronoi_field_combined <<- voronoi_field_combined } ################################################################################################################# ####### 5. calculate rustic averages ############################################################################### #RUSTIC <- function(){ #rustic1 <- as.data.frame (aggregate(RLD ~ Depth + Rep, data = file_name, FUN= "mean" )) #rustic <- cbind(rustic1, rep_year_crop1$year, rep_year) #names(rustic)[4]<- "type" #rustic$method <- "rustic" #write.csv(rustic, paste(type, "rusticbyrep.csv", sep ="")) #} for (year in seasons){ for (Crop in crops){ load_funcs() setdir() layout() vor_by_rep() # RUSTIC() } } ############ Combining location data with the voronoi data ############## combined_2015_2016_raw <- data.frame() for (year in seasons){ setdir() file_name$year <- year combined_2015_2016_raw <- rbind(file_name, combined_2015_2016_raw) combined_2015_2016_raw1 <- combined_2015_2016_raw[,-c(5,6)] } voronoi_by_rep$Location <- "voronoi" All_data_field <- rbind(combined_2015_2016_raw1, voronoi_by_rep) All_data_field write.csv(All_data_field, paste("D:/DES/Penn State/roots lab/Coring_study_Jimmy/gitlab/TOST_172020/All_data_field_and_voronoibyrep.csv", sep ="")) summary(All_data_field) ############################################################################################# ## After running this code you should have: ## 1. A combined file (All_data_field_and_voronoibyrep.csv) which includes all of the field data and voronoi estimated averages. ##