####Script for calculating the extincton coefficient (k) for different wheat genotypes
#####Joost van Heerwaarden, Wageningen University, 20/10/2015

##read data ##do complicated things because the input data has a complicated format

#####Part1#########

# read crop data
crop.data<-read.csv("DATA/RAW/Crop data.csv", na.strings=c("-",""),as.is=T,skip=3,header=F)

#read information rows of data seperately,  extract name row and set names for data columns
names<-read.csv("DATA/RAW/Crop data.csv", na.strings=c("-",""),as.is=T, nrows=3)
new.names<-as.character(names[1,])
colnames(crop.data)<-new.names

##read PAR data
par.data<-read.csv("DATA/RAW/AccuPAR data.csv", na.strings=c("-",""),as.is=T,skip=3,header=F)

#read information rows of data seperately,  extract name row and set names for data columns
names<-read.csv("DATA/RAW/AccuPAR data.csv", na.strings=c("-",""),as.is=T, nrows=3)
new.names<-as.character(names[1,])
colnames(par.data)<-new.names


###subset data to correspond to dates in report?
dates<-c("31/03/14","22/04/14","07/05/14","26/05/14","12/06/14")
crop.data<-crop.data[crop.data$Date%in%dates,]
par.data<-par.data[par.data$'Harvest Date'%in%dates,]


#subset correct interception data
par.data<-par.data[par.data$Measurement=="Light interception",]

#####Part1##########end


################################################

#####Part2##########calculate required variables from RAW data


###caclulate dry weight (g) of sample (by summing individual dry weights)
crop.data$DWS1<-rowSums(crop.data[,c("DWDL","DWGL","DWST","DWEA")],na.rm=T) #

#caclulate dry weight per 0.5m2 ground (by summing total dry weight of sample and of 0.5m2 plot)
crop.data$DWT05m<-rowSums(crop.data[,c("DWS1","DWOD")],na.rm=T) ##DWOD = dry weight per plot - sample

##in some cases a different calculation was made because of use of subsample for the calculation of dry weights...!!!!!!
#identify cases where different calculation is made
sel<-which(crop.data$DWOD!=crop.data$DWMP) # cases where dry weight is not the same as dry weight major part (i.e. where a sample was taken)

#do alternative calculation of dry weight per 0.5m2
crop.data$DWT05m.alt<-crop.data$FWT*(crop.data$DWOD/crop.data$FWOD) ##dry weight of total 0.5m2 calculated based on relation sample fresh/dry weight

##set appropriate values to alternative calculation
#crop.data$DWT05m[sel]<-crop.data$DWT05m.alt[sel]


#calculate LAI based on leaf area per sample (LA, cm-2 per sample),  DWT (g) per sample en  en DWT (g) per 0.5 m-2

##calculate ground area per sample
ground.area.per.sample<-(crop.data$DWS1/crop.data$DWT05m)*(0.5*10000)##calculate 1/ground area per sample (taken from) in cm-2

##calculate LAI based on leaf area and ground area
crop.data$LAIobs<-(crop.data$LA/ground.area.per.sample)

##set 0 LAI values to missing data 
crop.data$LAIobs[crop.data$LAI==0]<-NA #set 0 to NA

##make temporary dataframe with accuPAR readings (µmol m-2 m-1) per segment 
temp<-par.data[,grep("Segment",colnames(par.data))]

#calculate means over segments for each row in the dataframe
par.data$I.av<-apply(temp,1,mean,na.rm=T)

#remove temporary object
rm(temp)

##calculate averages per date per variety

##make vector with variety/date combinations
grouping.vector<-paste(crop.data$Variety,crop.data$'Date',sep=" ")

##make vector of means per variety/date combination for LAI
LAI.mean<-tapply(crop.data$LAIobs,grouping.vector,mean,na.rm=T)

##make vector with variety/date combinations
grouping.vector<-paste(par.data$Variety,par.data$'Harvest Date',sep=" ")

##make vector of means per variety/date combination for, I0, I (µmol m-2 m-1) and LAI
I0.mean<-tapply(par.data$'External Sensor',grouping.vector,mean,na.rm=T)
I.mean<-tapply(par.data$I.av,grouping.vector,mean,na.rm=T)

##match order of elements in vector of LAI with other vectors
LAI.mean<-LAI.mean[match(names(I0.mean),names(LAI.mean))]

##make new data frame with variety, date, I0,I,LAI en k

###make vector of varieties and dates (by splitting  names)
temp<-strsplit(names(I0.mean)," ")
temp<-unlist(temp)

###convert into a matrix containing varieties and dates
temp<-matrix(temp,,2,byrow=T)
colnames(temp)<-c("variety","date")


###make new data frame that holds mean values of I0, I and LAI for each variety/date
new.dat<-data.frame(temp,I0.mean,I.mean,LAI.mean)

#write data to file
write.csv(new.dat,"DATA/PROCESSED/new.dat.csv")

#####Part2##########end



#####Part3##########do calculations and generate results


#read the data back again

new.dat<-read.csv("DATA/PROCESSED/new.dat.csv",row.names=1)

##########Do regression (with 0 intercept) between -ln(I/I0) en LAI to calculate k for each variety

##make seperate file for each variety
temp.Julius<-new.dat[new.dat$variety=="Julius",]
temp.Ritmo<-new.dat[new.dat$variety=="Ritmo",]
temp.Tabasco<-new.dat[new.dat$variety=="Tabasco",]

#perform regression analysis for each variety seperately 
lm1<-lm(-log(I.mean/I0.mean)~0+LAI.mean,data=temp.Julius)
lm2<-lm(-log(I.mean/I0.mean)~0+LAI.mean,data=temp.Ritmo)
lm3<-lm(-log(I.mean/I0.mean)~0+LAI.mean,data=temp.Tabasco)

##extract regression coefficients and round to two decimals
coef1<-round(coef(lm1),2)
coef2<-round(coef(lm2),2)
coef3<-round(coef(lm3),2)

####write regression results to file for use in report

##open a text file to write to
sink(file="RESULTS/k.regression.txt",append=FALSE)

###print regression results
print(summary(lm1))
print(summary(lm2))
print(summary(lm3))

#close the file
sink() 


##now do pretty plot for k 
##make vectors of (varieties and corresponding colors 
vts<-new.dat$variety
cols<-c("black","green","red")[as.numeric(vts)]
pchs<-as.numeric(vts)

###open a pdf file to write the plot to
out<-paste("RESULTS/k_regression.pdf",sep="")
pdf(file=out, onefile=TRUE)

###plot -ln(I/I0) as function of LAI (with variety-specific colors)
plot(-log(I.mean/I0.mean)~LAI.mean,data=new.dat,ylim=c(0,5),pch=pchs,col=cols,xlim=c(0,6))

###add regression lines (with variety-specific colors) 
abline(c(0,coef(lm1)),col=unique(cols)[1])
abline(c(0,coef(lm2)),col=unique(cols)[2])
abline(c(0,coef(lm3)),col=unique(cols)[3])

#add value of k for each variety (with variety-specific colors)
text(5,1.5,label=paste(unique(vts)[1]," = ", coef1),col=unique(cols)[1])
text(5,1,label=paste(unique(vts)[2]," = ", coef2),col=unique(cols)[2])
text(5,0.5,label=paste(unique(vts)[3]," = ", coef3),col=unique(cols)[3])

####close pdf file
dev.off()

########
#####Part3##########ens


####THE END!
#(now run it all in one go)


