###this script does the statistical analyses using the files created with the first script: processing NDVI ad thermal imagery

####this script is meant to analyse the pilot field data collected in south africa in 2018

#final files for use: dat_all (contains the final data for use in analyses)
#hemi: canopy lai at plot level: hemi_results_meanValues.csv

rm(list=ls())


# ####(1) Install packages ------------------------------------------------
necessary.packages<-c("doBy","reshape2","MASS","pixmap","ggplot2","MuMIn","rimage","biOps","tseries","raster","graphics","rtiff","jpeg","gamlss","Kendall","R.utils","plyr","reshape")
already.installed<-necessary.packages%in%installed.packages()[, 'Package']    #asks if the necessary packages are already installed in the library?
if (length(necessary.packages[!already.installed])>=1) {    #if any are NOT installed, download them now.
  install.packages(necessary.packages[!already.installed],dep=T)   #are the dependencies really necessary (there are lots!)?
}
sapply(necessary.packages,function(p) {require(p,quietly=T,character.only=T)})

length1 <- function (x) {
  length(!is.na(x)[!is.na(x)=="TRUE"])
}

stderr <- function(x, na.rm=FALSE) {
  if(is.matrix(x))
    apply(x,2,stderr,na.rm=rm)
  else if(is.vector(x))
    sqrt(var(x,na.rm=na.rm)/length1(x))
  else if(is.data.frame(x))
    sapply(x,stderr,na.rm=na.rm)
  else sqrt(var(as.vector(x),na.rm=na.rm)/length1(x))
}

remove_outliers <- function(x, na.rm = TRUE, ...) {
  qnt <- quantile(x, probs=c(.25, .75), na.rm = na.rm, ...)
  H <- 1.5 * IQR(x, na.rm = na.rm)
  y <- x
  y[x < (qnt[1] - H)] <- NA
  y[x > (qnt[2] + H)] <- NA
  y
}

summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  library(plyr)
  
  # New version of length which can handle NA's: if na.rm==T, don't count them
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  
  # This does the summary. For each group's data frame, return a vector with
  # N, mean, and sd
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm)
                   )
                 },
                 measurevar
  )
  
  # Rename the "mean" column    
  datac <- rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  # Confidence interval multiplier for standard error
  # Calculate t-statistic for confidence interval: 
  # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}

library(psych)
library(mgcv)

# ###(2) link to working directory and create output directories-----------------------------------------

wdname<-"D:/data_peerj"
setwd(wdname)

dirG<-paste(wdname,"R_Outputs",sep="/")
if (file.exists(dirG)){
  cat("\n check for directory Results")
} else {
  dir.create(file.path(dirG))
}

dirGT<-paste(wdname,"R_Outputs/Tables",sep="/")
if (file.exists(dirGT)){
  cat("\n check for directory Results")
} else {
  dir.create(file.path(dirGT))
}

dirGG<-paste(wdname,"R_Outputs/Graphs",sep="/")
if (file.exists(dirGG)){
  cat("\n check for directory Results")
} else {
  dir.create(file.path(dirGG))
}

# #####(3) Read data files ------------------------------------------------
dat_all<-read.csv(paste(wdname,"dat_all_PeerJ.csv",sep="/"))
nhem<-read.csv(paste(wdname,"hemi_results_meanValues.csv",sep="/"))
ov<-read.csv(paste(wdname,"SensorData_overview.csv",sep="/"))

# ####(4) analyse habitat quality patterns -------------------------------

##ndvi statistics, aggregate across habitats
#nav1: NDVI of the ground vegetation
#nav_s: NDVI of the canopy as seen from the ground
t1<-data.frame(nav1=dat_all$nav1, Habitat=dat_all$Habitat,plot=dat_all$plot)
t1<-na.omit(t1)

t2<-data.frame(nav_s=dat_all$nav_s, Habitat=dat_all$Habitat,plot=dat_all$plot)
t2<-na.omit(t2)

tb<-na.omit(data.frame(nav_s=dat_all$nav_s,nav1=dat_all$nav1,Habitat=dat_all$Habitat))
describe.by(tb$nav1, tb$Habitat)
length(rownames(na.omit(data.frame(dat_all$nav_s,dat_all$nav1)))) #108

tgc <- summarySE(t2, measurevar="nav_s", groupvars=c("Habitat","plot"))
tgc1<- summarySE(t1, measurevar="nav1", groupvars=c("Habitat","plot"))
tgc$ndvi_g_se<-tgc1$se[match(tgc$plot,tgc1$plot)]
tgc$ndvi_g_av<-tgc1$nav1[match(tgc$plot,tgc1$plot)]
names(tgc)[4]<-"ndvi_s_av"

#test models linking canopy NDVI to ground NDVI
describe.by(tgc$ndvi_g_av, tgc$Habitat)
summary(gam(ndvi_g_av~s(ndvi_s_av),data=tgc))

Colourfig2<- c("#E69F00", "#56B4E9","#009E73", "#CC79A7")
fv<-ggplot(tgc, aes(ndvi_s_av, ndvi_g_av,col=Habitat))+ geom_point(size=4)+
  geom_errorbarh(aes(xmin=ndvi_s_av-se, xmax=ndvi_s_av+se))+geom_errorbar(aes(ymin=ndvi_g_av-ndvi_g_se,ymax=ndvi_g_av+ndvi_g_se, width=.1)) 
fv<-fv+ ylim(-0.6, 1)
fv<-fv+geom_point(data=dat_all, aes(nav_s,nav1),col="grey")
fv<-fv+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fv<-fv+ guides(color=guide_legend(override.aes=list(fill=NA)))
fv<-fv+scale_color_manual(values=Colourfig2)
#fvc<-fv+stat_smooth(aes(size=1.5),method="gam",formula = y ~ s(x))
fvc<-fv+stat_smooth(data=tgc,fill="grey",col="darkgrey",aes(x=ndvi_s_av,y=ndvi_g_av),method="gam",formula = y ~ s(x))
fvc<-fvc+xlab("mean NDVI of the canopy (NDVIup)") +ylab("mean NDVI of the ground (NDVIdown)")
fvc<-fvc+theme(axis.title.x=element_text(vjust=-0.35,size=26),axis.title.y=element_text(vjust= +0.9,size=26),
               axis.text.x=element_text(size=26),axis.text.y=element_text(size=24),
               legend.title = element_text(size=26, face="bold"),
               legend.text = element_text(size = 26))
fvc<-fvc + theme(legend.position = c(0.8, 0.2))
fvc<-fvc+guides(size = FALSE)
nam.p<-paste("NDVI_up_down","Acrosshabitats1","tiff",sep=".")
ggsave(fvc,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)
summary(gam(ndvi_g_av~s(ndvi_s_av),data=tgc))
summary(lm(ndvi_g_av~ndvi_s_av,data=tgc))

#use ndvi ground and fractional vegetation cover (%) = fcover data
#tgc <- summarySE(hq1, measurevar="ndvi_s_av", groupvars=c("Habitat","plot"))
t3<-data.frame(FCover=dat_all$FCover, Habitat=dat_all$Habitat,plot=dat_all$plot)
t3<-na.omit(t3)
describe.by(t3$FCover, t3$Habitat)
tgc2<- summarySE(t3, measurevar="FCover", groupvars=c("Habitat","plot"))
names(tgc1)[4]<-"ndvi_g_av"
names(tgc1)[6]<-"ndvi_g_se"
names(tgc2)
tgc1$fco_se<-tgc2$se[match(tgc1$plot,tgc2$plot)]
tgc1$fco_av<-tgc2$FCover[match(tgc1$plot,tgc2$plot)]
#describe.by(dat_all$FCover,dat_all$Habitat)
tgc1<-droplevels(na.omit(tgc1))


Colourfig2<- c("#E69F00", "#56B4E9","#009E73", "#CC79A7")
fv<-ggplot(tgc1, aes(fco_av, ndvi_g_av,col=Habitat))+ geom_point(size=4)+
  geom_errorbar(aes(ymin=ndvi_g_av-ndvi_g_se, ymax=ndvi_g_av+ndvi_g_se,width=.1))+geom_errorbarh(aes(xmin=fco_av-fco_se,xmax=fco_av+fco_se, width=.1)) 
fv<-fv+geom_point(data=dat_all, aes(FCover, nav1),col="grey")
fv<-fv+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fv<-fv+ guides(color=guide_legend(override.aes=list(fill=NA)))
fv<-fv+scale_color_manual(values=Colourfig2)
#fvc<-fv+stat_smooth(aes(size=1.5),method="gam",formula = y ~ s(x,k=2))
fvc<-fv+stat_smooth(data=tgc1,fill="grey",col="darkgrey",aes(x=fco_av,y=ndvi_g_av),method="glm",formula = y ~ x)
fvc<-fvc+xlab("mean fractional vegetation cover (FCover)") +ylab("mean NDVI of the ground (NDVIdown)")
fvc<-fvc+theme(axis.title.x=element_text(vjust=-0.35,size=26),axis.title.y=element_text(vjust= +0.9,size=26),
               axis.text.x=element_text(size=26),axis.text.y=element_text(size=26),
               legend.title = element_text(size=26, face="bold"),
               legend.text = element_text(size = 26))
fvc<-fvc + theme(legend.position = c(0.2, 0.2))
fvc<-fvc+guides(size = FALSE)
nam.p<-paste("FCo_NDVI_Ground","Acrosshabitats1","tiff",sep=".")
ggsave(fvc,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)

write.csv(tgc1,paste(dirGT,"tgc1_fcover_ndvi.csv",sep="/"))
summary(loess(ndvi_g_av~fco_av,tgc1))
summary(lm(ndvi_g_av~fco_av,tgc1))
summary(gam(ndvi_g_av~s(fco_av),data=tgc1))

###use thermal data = ground surface temperature
describe.by(dat_all$tav, dat_all$Habitat)
t4<-data.frame(tav=dat_all$tav, Habitat=dat_all$Habitat,plot=dat_all$plot)
t4<-na.omit(t4)
describe.by(t4$tav, t4$Habitat)

#tgc <- summarySE(hq1, measurevar="ndvi_s_av", groupvars=c("Habitat","plot"))
tgc3<- summarySE(t4, measurevar="tav", groupvars=c("Habitat","plot"))
tgc2$tav_se<-tgc3$se[match(tgc2$plot,tgc3$plot)]
tgc2$tav_av<-tgc3$tav[match(tgc2$plot,tgc3$plot)]
tgc2<-na.omit(tgc2)
#remove the crop plot
tgc2<-droplevels(na.omit(tgc2))

Colourfig2<- c("#E69F00", "#56B4E9","#009E73", "#CC79A7")
fv<-ggplot(tgc2, aes(FCover, tav_av,col=Habitat))+ geom_point(size=4)+
  geom_errorbar(aes(ymin=tav_av-tav_se, ymax=tav_av+tav_se,width=.1))+geom_errorbarh(aes(xmin=FCover-se,xmax=FCover+se, width=.1)) 
fv<-fv+geom_point(data=dat_all, aes(FCover, tav),col="grey")
fv<-fv+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fv<-fv+ guides(color=guide_legend(override.aes=list(fill=NA)))
fv<-fv+scale_color_manual(values=Colourfig2)
#fvc<-fv+stat_smooth(aes(size=1.5),method="gam",formula = y ~ s(x,k=2))
fvc<-fv+stat_smooth(data=tgc2,fill="grey",col="darkgrey",aes(x=FCover,y=tav_av),method="glm",formula = y ~ x)
fvc<-fvc+xlab("mean fractional vegetation cover (FCover)") +ylab("Surface temperature (Thermalground)")
fvc<-fvc+theme(axis.title.x=element_text(vjust=-0.35,size=26),axis.title.y=element_text(vjust= +0.9,size=26),
               axis.text.x=element_text(size=26),axis.text.y=element_text(size=26),
               legend.title = element_text(size=26, face="bold"),
               legend.text = element_text(size = 26))
fvc<-fvc + theme(legend.position = c(0.2, 0.8))
fvc<-fvc+guides(size = FALSE)
nam.p<-paste("Thermal_Fcover","Acrosshabitats1","tiff",sep=".")
ggsave(fvc,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)

write.csv(tgc2,paste(dirGT,"tgc2_fcover_tav.csv",sep="/"))

summary(loess(tav_av~FCover,tgc2))
summary(lm(tav_av~FCover,tgc2))
summary(gam(tav_av~s(FCover),data=tgc2))


##attach thermal data to ndvi ground data

tgc4<- summarySE(t1, measurevar="nav1", groupvars=c("Habitat","plot"))
names(tgc4)[4]<-"ndvi_g_av"
names(tgc4)[6]<-"ndvi_g_se"
tgc4$tav_se<-tgc3$se[match(tgc4$plot,tgc3$plot)]
tgc4$tav_av<-tgc3$tav[match(tgc4$plot,tgc3$plot)]
tgc4<-na.omit(tgc4)

Colourfig2<- c("#E69F00", "#56B4E9","#009E73", "#CC79A7")
Colourfig2a<- c("red","#E69F00", "#56B4E9","#009E73", "#CC79A7")
tgc4$Habitat<-factor(tgc4$Habitat, levels = c("Grassland", "Bush", "Edge","Forest","Plantation"))

fv<-ggplot(tgc4, aes(ndvi_g_av, tav_av,col=Habitat))+ geom_point(size=4)+
  geom_errorbar(aes(ymin=tav_av-tav_se, ymax=tav_av+tav_se,width=.1))+geom_errorbarh(aes(xmin=ndvi_g_av-ndvi_g_se,xmax=ndvi_g_av+ndvi_g_se)) 
fv<-fv+geom_point(data=dat_all, aes(nav1, tav),col="grey")
fv<-fv+xlim(-1,1)
fv<-fv+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fv<-fv+ guides(color=guide_legend(override.aes=list(fill=NA)))
fv<-fv+scale_color_manual(values=Colourfig2a)
#fvc<-fv+stat_smooth(aes(size=1.5),method="gam",formula = y ~ s(x,k=2))
fvc<-fv+stat_smooth(data=tgc4,fill="grey",col="darkgrey",aes(x=ndvi_g_av,y=tav_av),method="glm",formula = y ~ x)
fvc<-fvc+xlab("mean NDVI of the ground (NDVIdown)") +ylab("Surface temperature (Thermalground)")
fvc<-fvc+theme(axis.title.x=element_text(vjust=-0.35,size=26),axis.title.y=element_text(vjust= +0.9,size=26),
               axis.text.x=element_text(size=26),axis.text.y=element_text(size=26),
               legend.title = element_text(size=26, face="bold"),
               legend.text = element_text(size = 26))
fvc<-fvc + theme(legend.position = c(0.2, 0.8))
fvc<-fvc+guides(size = FALSE)
nam.p<-paste("Thermal_NDVI_ground","Acrosshabitats1","tiff",sep=".")
ggsave(fvc,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)

summary(gam(tav_av~s(ndvi_g_av),data=tgc4))
summary(lm(tav_av~ndvi_g_av,data=tgc4))

tgc5 <- summarySE(t2, measurevar="nav_s", groupvars=c("Habitat","plot"))
tgc3<- summarySE(t4, measurevar="tav", groupvars=c("Habitat","plot"))
names(tgc5)[4]<-"ndvi_s_av"
names(tgc5)[6]<-"ndvi_s_se"
names(tgc3)[4]<-"tav_av"
names(tgc3)[6]<-"tav_se"
tgc5$tav_se<-tgc3$tav_se[match(tgc5$plot,tgc3$plot)]
tgc5$tav_av<-tgc3$tav_av[match(tgc5$plot,tgc3$plot)]
tgc5<-na.omit(tgc5)

Colourfig2<- c("#E69F00", "#56B4E9","#009E73", "#CC79A7")
fv<-ggplot(tgc5, aes(ndvi_s_av, tav_av,col=Habitat))+ geom_point(size=4)+
  geom_errorbar(aes(ymin=tav_av-tav_se, ymax=tav_av+tav_se,width=.1))+geom_errorbarh(aes(xmin=ndvi_s_av-ndvi_s_se,xmax=ndvi_s_av+ndvi_s_se)) 
fv<-fv+geom_point(data=dat_all, aes(nav_s, tav),col="grey")
#fv<-fv+xlim(-1,1)
fv<-fv+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fv<-fv+ guides(color=guide_legend(override.aes=list(fill=NA)))
fv<-fv+scale_color_manual(values=Colourfig2)
#fvc<-fv+stat_smooth(aes(size=1.5),method="gam",formula = y ~ s(x,k=2))
fvc<-fv+stat_smooth(data=tgc5,fill="grey",col="darkgrey",aes(x=ndvi_s_av,y=tav_av),method="gam",formula = y ~ x)
fvc<-fvc+xlab("mean NDVI of the canopy (NDVIup)") +ylab("Surface temperature (Thermalground)")
fvc<-fvc+theme(axis.title.x=element_text(vjust=-0.35,size=26),axis.title.y=element_text(vjust= +0.9,size=26),
               axis.text.x=element_text(size=26),axis.text.y=element_text(size=26),
               legend.title = element_text(size=26, face="bold"),
               legend.text = element_text(size = 26))
fvc<-fvc + theme(legend.position = c(0.2, 0.8))
fvc<-fvc+guides(size = FALSE)
nam.p<-paste("Thermal_NDVI_sky","Acrosshabitats1","tiff",sep=".")
ggsave(fvc,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)
summary(gam(tav_av~s(ndvi_s_av,k=2),data=tgc5))
summary(lm(tav_av~ndvi_s_av,data=tgc5))

#attach canopy lai data: these come at plot level only
tgc3$lai<-nhem$LAIt[match(tgc3$plot,nhem$folder)]
tgc6<-droplevels(na.omit(tgc3))
tgc7 <- summarySE(t2, measurevar="nav_s", groupvars=c("Habitat","plot"))
names(tgc7)[4]<-"ndvi_s_av"
names(tgc7)[6]<-"ndvi_s_se"
tgc7$lai<-nhem$LAIt[match(tgc7$plot,nhem$folder)]
tgc7<-droplevels(na.omit(tgc7))


#Colourfig2<- c("#E69F00", "#56B4E9","#009E73", "#CC79A7")
Colourfig2b<- c("#E69F00", "#56B4E9","#009E73", "#CC79A7","yellow")
tgc6$Habitat<-factor(tgc6$Habitat, levels = c("Bush", "Edge","Forest","Plantation","Crop"))
fv<-ggplot(tgc6, aes(lai, tav_av,col=Habitat))+ geom_point(size=4)+
  geom_errorbar(aes(ymin=tav_av-tav_se, ymax=tav_av+tav_se,width=.1))
#fv<-fv+geom_point(data=hq1, aes(ndvi_s_av, tav),col="grey")
fv<-fv+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fv<-fv+ guides(color=guide_legend(override.aes=list(fill=NA)))
fv<-fv+scale_color_manual(values=Colourfig2b)
#fvc<-fv+stat_smooth(aes(size=1.5),method="gam",formula = y ~ s(x,k=2))
fvc<-fv+stat_smooth(data=tgc6,fill="grey",col="darkgrey",aes(x=lai,y=tav_av),method="glm",formula = y ~ x)
fvc<-fvc+xlab("Leaf Area Index") +ylab("Surface temperature (Thermalground)")
fvc<-fvc+theme(axis.title.x=element_text(vjust=-0.35,size=26),axis.title.y=element_text(vjust= +0.9,size=26),
               axis.text.x=element_text(size=26),axis.text.y=element_text(size=26),
               legend.title = element_text(size=26, face="bold"),
               legend.text = element_text(size = 26))
fvc<-fvc + theme(legend.position = c(0.2, 0.2))
fvc<-fvc+guides(size = FALSE)
nam.p<-paste("Thermal_LAI","Acrosshabitats1","tiff",sep=".")
ggsave(fvc,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)


summary(gam(tav_av~s(lai),data=tgc6))
summary(lm(tav_av~lai,data=tgc6))
write.csv(tgc6,paste(dirGT,"tgc6_lai_thermal.csv",sep="/"))


plot(tgc6$lai,tgc6$tav_av)
#identify(tgc6$lai,tgc6$tav_av,cex=0.7,tgc6$plot)

Colourfig2<- c("#E69F00", "#56B4E9","#009E73", "#CC79A7")
fv<-ggplot(tgc7, aes(lai, ndvi_s_av,col=Habitat))+ geom_point(size=4)+
  geom_errorbar(aes(ymin=ndvi_s_av-ndvi_s_se, ymax=ndvi_s_av+ndvi_s_se,width=.1))
#fv<-fv+geom_point(data=hq1, aes(ndvi_s_av, tav),col="grey")
fv<-fv+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fv<-fv+ guides(color=guide_legend(override.aes=list(fill=NA)))
fv<-fv+scale_color_manual(values=Colourfig2)
#fvc<-fv+stat_smooth(aes(size=1.5),method="gam",formula = y ~ s(x,k=2))
fvc<-fv+stat_smooth(data=tgc7,fill="grey",col="darkgrey",aes(x=lai,y=ndvi_s_av),method="glm",formula = y ~ x)
fvc<-fvc+xlab("Leaf Area Index") +ylab("NDVIup")
fvc<-fvc+theme(axis.title.x=element_text(vjust=-0.35,size=20),axis.title.y=element_text(vjust= +0.9,size=26),
               axis.text.x=element_text(size=26),axis.text.y=element_text(size=26),
               legend.title = element_text(size=26, face="bold"),
               legend.text = element_text(size = 26))
fvc<-fvc + theme(legend.position = c(0.2, 0.2))
fvc<-fvc+guides(size = FALSE)
nam.p<-paste("NDVIup_LAI","Acrosshabitats1","tiff",sep=".")
ggsave(fvc,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)

summary(gam(ndvi_s_av~s(lai),data=tgc7))
cor.test(tgc7$lai,tgc7$ndvi_s_av,method="spearman")    
cor.test(tgc7$lai,tgc7$ndvi_s_av,method="spearman")    


tgc8<- summarySE(t1, measurevar="nav1", groupvars=c("Habitat","plot"))
names(tgc8)[4]<-"ndvi_g_av"
names(tgc8)[6]<-"ndvi_g_se"
tgc8$lai<-nhem$LAIt[match(tgc8$plot,nhem$folder)]
tgc8<-droplevels(na.omit(tgc8))

Colourfig2<- c("#E69F00", "#56B4E9","#009E73", "#CC79A7")
fv<-ggplot(tgc8, aes(lai, ndvi_g_av,col=Habitat))+ geom_point(size=4)+
  geom_errorbar(aes(ymin=ndvi_g_av-ndvi_g_se, ymax=ndvi_g_av+ndvi_g_se,width=.1))
#fv<-fv+geom_point(data=hq1, aes(ndvi_s_av, tav),col="grey")
fv<-fv+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fv<-fv+ guides(color=guide_legend(override.aes=list(fill=NA)))
fv<-fv+scale_color_manual(values=Colourfig2)
#fvc<-fv+stat_smooth(aes(size=1.5),method="gam",formula = y ~ s(x,k=2))
fvc<-fv+stat_smooth(data=tgc8,fill="grey",col="darkgrey",aes(x=lai,y=ndvi_g_av),method="glm",formula = y ~ x)
fvc<-fvc+xlab("Leaf Area Index") +ylab("NDVIground")
fvc<-fvc+theme(axis.title.x=element_text(vjust=-0.35,size=26),axis.title.y=element_text(vjust= +0.9,size=26),
               axis.text.x=element_text(size=26),axis.text.y=element_text(size=26),
               legend.title = element_text(size=26, face="bold"),
               legend.text = element_text(size = 26))
fvc<-fvc + theme(legend.position = c(0.2, 0.2))
fvc<-fvc+guides(size = FALSE)
nam.p<-paste("NDVIground_LAI","Acrosshabitats1","tiff",sep=".")
ggsave(fvc,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)

write.csv(tgc8,paste(dirGT,"tgc8_lai_ndvig.csv",sep="/"))

summary(gam(ndvi_g_av~s(lai),data=tgc8))
summary(lm(ndvi_g_av~lai,data=tgc8))
plot(tgc8$lai,tgc8$ndvi_g_av)
#identify(tgc8$lai,tgc8$ndvi_g_av,cex=0.7,tgc8$plot)

# ##(5) Test for correlation between the many variables -------------------
big.d<-dat_all[,c(7,8,10,11,14,15,17,18,21,22,24,25,34)]
#at level of sample points
resa<-cor(big.d,method="spearman",use="complete.obs")
resa<-round(resa,2)

# pairwise sample size
# Gabor G - 11/23/2004 R-help List

corr.test(big.d,method="spearman",use="complete.obs")

cor.test(tgc8$lai,tgc8$ndvi_g_av,method="spearman")    
cor.test(tgc7$lai,tgc7$ndvi_s_av,method="spearman")    
cor.test(tgc6$lai,tgc6$tav_av,method="spearman")    

# [1] "nav1"     "nmedian1" "nmax1"    "nmin1"    "nav_s"    "nmed_s"   "nmax_s"   "nmin_s"   "tav"      "tmedian"  "tmax"     "tmin"    
#[13] "FCover
plot(nmedian1~nav1,big.d)
summary(lm(nmedian1~nav1,big.d))
plot(nmax1~nav1,big.d)
summary(lm(nmax1~nav1,big.d))
plot(nmax1~nmedian1,big.d)
summary(lm(nmax1~nmedian1,big.d))

plot(nav_s~nav1,big.d)
summary(lm(nav_s~nav1,big.d))
plot(nav_s~nmedian1,big.d)
summary(lm(nav_s~nmedian1,big.d))
plot(nmed_s~nav1,big.d)
summary(lm(nmed_s~nav1,big.d))
plot(nmed_s~nmedian1,big.d)
summary(lm(nmed_s~nmedian1,big.d))
plot(nav_s~nmed_s,big.d)
summary(lm(nav_s~nmed_s,big.d))

plot(nmax_s~nav_s,big.d)
summary(lm(nmax_s~nav_s,big.d))
plot(nmax_s~nmed_s,big.d)
summary(lm(nmax_s~nmed_s,big.d))

plot(nmin_s~nav1,big.d)
summary(lm(nmin_s~nav1,big.d))
plot(nmin_s~nmedian1,big.d)
summary(lm(nmin_s~nmedian1,big.d))

plot(nmin_s~nav_s,big.d)
summary(lm(nmin_s~nav_s,big.d))
plot(nmin_s~nmax_s,big.d)
summary(lm(nmin_s~nmax_s,big.d))

# [1] "nav1"     "nmedian1" "nmax1"    "nmin1"    "nav_s"    "nmed_s"   "nmax_s"   "nmin_s"   "tav"      "tmedian"  "tmax"     "tmin"    
#[13] "FCover
plot(tav~nav1,big.d)
summary(lm(tav~nav1,big.d))
plot(tav~nmedian1,big.d)
summary(lm(tav~nmedian1,big.d))
plot(tav~nmax_s,big.d)
summary(lm(tav~nmax_s,big.d))
plot(tav~nmin_s,big.d)
summary(lm(tav~nmin_s,big.d))

plot(tmedian~nav1,big.d)
summary(lm(tmedian~nav1,big.d))
plot(tmedian~nmedian1,big.d)
summary(lm(tmedian~nmedian1,big.d))
plot(tmedian~nmax_s,big.d)
summary(lm(tmedian~nmax_s,big.d))
plot(tmedian~nmin_s,big.d)
summary(lm(tmedian~nmin_s,big.d))
plot(tmedian~tav,big.d)
summary(lm(tmedian~tav,big.d))

plot(tmax~nav1,big.d)
summary(lm(tmax~nav1,big.d))
plot(tmax~nmedian1,big.d)
summary(lm(tmax~nmedian1,big.d))
plot(tmax~nmax_s,big.d)
summary(lm(tmax~nmax_s,big.d))
plot(tmax~nmin_s,big.d)
summary(lm(tmax~nmin_s,big.d))
plot(tmax~tav,big.d)
summary(lm(tmax~tav,big.d))
plot(tmax~tmedian,big.d)
summary(lm(tmax~tmedian,big.d))

plot(tmin~nav1,big.d)
summary(lm(tmin~nav1,big.d))
plot(tmin~nmedian1,big.d)
summary(lm(tmin~nmedian1,big.d))
plot(tmin~nmax_s,big.d)
summary(lm(tmin~nmax_s,big.d))
plot(tmin~nmin_s,big.d)
summary(lm(tmin~nmin_s,big.d))
plot(tmin~tav,big.d)
summary(lm(tmin~tav,big.d))
plot(tmin~tmedian,big.d)
summary(lm(tmin~tmedian,big.d))
plot(tmin~tmax,big.d)
summary(lm(tmin~tmax,big.d))

# [1] "nav1"     "nmedian1" "nmax1"    "nmin1"    "nav_s"    "nmed_s"   "nmax_s"   "nmin_s"   "tav"      "tmedian"  "tmax"     "tmin"    
#[13] "FCover

plot(FCover~nav1,big.d)
summary(lm(FCover~nav1,big.d))
plot(FCover~nmedian1,big.d)
summary(lm(FCover~nmedian1,big.d))
plot(FCover~nmed_s,big.d)
summary(lm(FCover~nmed_s,big.d))
plot(FCover~nav_s,big.d)
summary(lm(FCover~nav_s,big.d))
plot(FCover~nmax_s,big.d)
summary(lm(FCover~nmax_s,big.d))
plot(FCover~nmin_s,big.d)
summary(lm(FCover~nmin_s,big.d))
plot(FCover~tav,big.d)
summary(lm(FCover~tav,big.d))
plot(FCover~tmedian,big.d)
summary(lm(FCover~tmedian,big.d))
plot(FCover~tmax,big.d)
summary(lm(FCover~tmax,big.d))
plot(FCover~tmin,big.d)
summary(lm(FCover~tmin,big.d))


temp<-as.data.frame(aggregate(dat_all$nav1,list(dat_all$plot),min,na.rm=T) )    
tgc8$ndvi_g_min<-temp$x[match(tgc8$plot,temp$Group.1)]     
temp<-as.data.frame(aggregate(dat_all$nav1,list(dat_all$plot),max,na.rm=T) )  
tgc8$ndvi_g_max<-temp$x[match(tgc8$plot,temp$Group.1)] 
cor.test(tgc8$lai,tgc8$ndvi_g_max,method="spearman")    
cor.test(tgc8$lai,tgc8$ndvi_g_min,method="spearman")    


temp<-as.data.frame(aggregate(dat_all$nav_s,list(dat_all$plot),min,na.rm=T) )    
tgc7$ndvi_s_min<-temp$x[match(tgc7$plot,temp$Group.1)]     
temp<-as.data.frame(aggregate(dat_all$nav_s,list(dat_all$plot),max,na.rm=T) )  
tgc7$ndvi_s_max<-temp$x[match(tgc7$plot,temp$Group.1)] 
cor.test(tgc7$lai,tgc7$ndvi_s_max,method="spearman")    
cor.test(tgc7$lai,tgc7$ndvi_s_min,method="spearman")    

temp<-as.data.frame(aggregate(dat_all$tav,list(dat_all$plot),min,na.rm=T) )    
tgc6$tav_min<-temp$x[match(tgc6$plot,temp$Group.1)]     
temp<-as.data.frame(aggregate(dat_all$tav,list(dat_all$plot),max,na.rm=T) )  
tgc6$tav_max<-temp$x[match(tgc6$plot,temp$Group.1)] 
cor.test(tgc6$lai,tgc6$tav_max,method="spearman")    
cor.test(tgc6$lai,tgc6$tav_min,method="spearman")    

length(tgc6$lai)
length(tgc7$lai)
length(tgc8$lai)
# ###(6) Habitat dependencies of variables --------------------------------


#for the ground ndvi
windowsFonts(Times=windowsFont("Times New Roman"))
dat_ng<-droplevels(na.omit(dat_all[,c(1,2,3,7,36)]))
Colourfigng<- c("red","#E69F00", "#56B4E9","#009E73", "#CC79A7")

dat_ng$Habitat<-factor(dat_ng$Habitat, levels = c("Grassland","Bush", "Edge","Forest","Plantation"))


fig1<-ggplot(dat_ng, aes(x=Habitat, y=nav1,fill=Habitat)) + 
         geom_boxplot(notch=F)+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fig1<-fig1+scale_fill_manual(values=Colourfigng)
fig1<-fig1+ylab("NDVIground") # Set axis labels
fig1<-fig1+ theme(axis.text.y=element_text(family="Times",size=27,color="black"), 
                         axis.text.x=element_text(family="Times",size=27,face="bold",color="black"),
                         axis.title.y=element_text(family="Times",size=27,face="bold",margin=margin(0,20,0,0))) +
         theme( axis.title.x=element_blank())+
         theme(legend.title=element_blank(),legend.text = element_text(family="Times",colour="black", size = 27))
fig1<-fig1+theme(legend.position="none")
fig1<-fig1+theme(text=element_text(family="Times"))
nam.p<-paste("NDVIdown","Habitats_Boxes","tiff",sep=".")
ggsave(fig1,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)

##for the sky ndvi
Colourfigfc<- c("#E69F00", "#56B4E9","darkblue","#009E73", "#CC79A7")
dat_ns<-droplevels(na.omit(dat_all[,c(1,2,3,4,14,36)]))
dat_ns$Habitat<-factor(dat_ns$Habitat, levels = c("Bush", "Edge","Forest","Plantation"))
fig1<-ggplot(dat_ns, aes(x=Habitat, y=nav_s,fill=Habitat)) + 
  geom_boxplot(notch=F)+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fig1<-fig1+scale_fill_manual(values=Colourfigfc)
fig1<-fig1+ylab("NDVIup") # Set axis labels
fig1<-fig1+ theme(axis.text.y=element_text(family="Times",size=27,color="black"), 
                  axis.text.x=element_text(family="Times",size=27,face="bold",color="black"),
                  axis.title.y=element_text(family="Times",size=27,face="bold",margin=margin(0,20,0,0))) +
  theme( axis.title.x=element_blank())+
  theme(legend.title=element_blank(),legend.text = element_text(family="Times",colour="black", size = 27))
fig1<-fig1+theme(legend.position="none")
fig1<-fig1+theme(text=element_text(family="Times"))
nam.p<-paste("NDVIup","Habitats_Boxes","tiff",sep=".")
ggsave(fig1,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)

#for thermal ground
dat_th<-droplevels(na.omit(dat_all[,c(1,2,3,21,36)]))
dat_th$Habitat<-factor(dat_th$Habitat, levels = c("Grassland","Bush", "Edge","Forest","Plantation"))
fig1<-ggplot(dat_th, aes(x=Habitat, y=tav,fill=Habitat)) + 
  geom_boxplot(notch=F)+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fig1<-fig1+scale_fill_manual(values=Colourfigng)
fig1<-fig1+ylab("Thermalground") # Set axis labels
fig1<-fig1+ theme(axis.text.y=element_text(family="Times",size=27,color="black"), 
                  axis.text.x=element_text(family="Times",size=27,face="bold",color="black"),
                  axis.title.y=element_text(family="Times",size=27,face="bold",margin=margin(0,20,0,0))) +
  theme( axis.title.x=element_blank())+
  theme(legend.title=element_blank(),legend.text = element_text(family="Times",colour="black", size = 27))
fig1<-fig1+theme(legend.position="none")
fig1<-fig1+theme(text=element_text(family="Times"))
nam.p<-paste("Thermalground","Habitats_Boxes","tiff",sep=".")
ggsave(fig1,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)

#for FCover
dat_fc<-droplevels(na.omit(dat_all[,c(1,2,3,34,36)]))
dat_fc$Habitat<-factor(dat_fc$Habitat, levels = c("Bush", "Edge","Woodland","Forest","Plantation"))
Colourfigfc<- c("#E69F00", "#56B4E9","darkblue","#009E73", "#CC79A7")
fig1<-ggplot(dat_fc, aes(x=Habitat, y=FCover,fill=Habitat)) + 
  geom_boxplot(notch=F)+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fig1<-fig1+scale_fill_manual(values=Colourfigfc)
fig1<-fig1+ylab("FCover (%)") # Set axis labels
fig1<-fig1+ theme(axis.text.y=element_text(family="Times",size=27,color="black"), 
                  axis.text.x=element_text(family="Times",size=27,face="bold",color="black"),
                  axis.title.y=element_text(family="Times",size=27,face="bold",margin=margin(0,20,0,0))) +
  theme( axis.title.x=element_blank())+
  theme(legend.title=element_blank(),legend.text = element_text(family="Times",colour="black", size = 27))
fig1<-fig1+theme(legend.position="none")
fig1<-fig1+theme(text=element_text(family="Times"))
nam.p<-paste("FCover","Habitats_Boxes","tiff",sep=".")
ggsave(fig1,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)

##for LAI (data aggregated at plot level)

nhem$Biome<-ov$Biome[match(nhem$folder,ov$PlotID)]
describe.by(nhem$LAIt, nhem$Biome)
dat_lai<-droplevels(na.omit(nhem[,c(2,3,6,14)]))


#dat_fc<-droplevels(dat_fc[-which(dat_fc$Habitat=="Crop"),])
dat_lai$Habitat<-factor(dat_lai$Biome, levels = c("Bush", "Edge","Woodland","Forest","Plantation"))
Colourfigfc<- c("#E69F00", "#56B4E9","darkblue","#009E73", "#CC79A7")
fig1<-ggplot(dat_lai, aes(x=Habitat, y=LAIt,fill=Habitat)) + 
  geom_boxplot(notch=F)+theme_bw()+theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fig1<-fig1+scale_fill_manual(values=Colourfigfc)
fig1<-fig1+ylab("LAI") # Set axis labels
fig1<-fig1+ theme(axis.text.y=element_text(family="Times",size=27,color="black"), 
                  axis.text.x=element_text(family="Times",size=27,face="bold",color="black"),
                  axis.title.y=element_text(family="Times",size=27,face="bold",margin=margin(0,20,0,0))) +
  theme( axis.title.x=element_blank())+
  theme(legend.title=element_blank(),legend.text = element_text(family="Times",colour="black", size = 27))
fig1<-fig1+theme(legend.position="none")
fig1<-fig1+theme(text=element_text(family="Times"))
nam.p<-paste("LAI","Habitats_Boxes","tiff",sep=".")
ggsave(fig1,filename=paste(dirGG,nam.p,sep="/"),dpi=300,width=10,height=8)




###test for differences in pairs:
shapiro.test(rnorm(n=1000))#p = 0.2662 eaning p values above 0.05 indicate normal distribtuion
do.call("rbind", with(dat_fc, tapply(FCover, Habitat,
                                     function(x) unlist(shapiro.test(x)[c("statistic", "p.value")]))))
# bartlett.test(EI.sensitivity ~ labeln3, data = dfvdat) #parametric
# fligner.test(EI.sensitivity ~ labeln3, data = dfvdat) #non parametric
y<-dat_fc$FCover
g<-dat_fc$Habitat
boxplot(y~g)
pairwise.wilcox.test(y, g, p.adj = "bonf") 
summary(aov(y~g))
anova(lm(y~ g, dat_fc))
TukeyHSD(aov(y~g))

library(car)
leveneTest(FCover~Habitat,dat_fc)


y<-dat_lai$LAIt
g<-dat_lai$Habitat
boxplot(y~g)
pairwise.wilcox.test(y, g, p.adj = "bonf") 
summary(aov(y~g))
anova(lm(y~ g, dat_lai))
TukeyHSD(aov(y~g))
leveneTest(LAIt~Habitat,dat_lai)

y<-dat_th$tav
g<-dat_th$Habitat
boxplot(y~g)
pairwise.wilcox.test(y, g, p.adj = "bonf") 
summary(aov(y~g))
anova(lm(y~ g, dat_th))
TukeyHSD(aov(y~g))

y<-dat_ns$nav_s
g<-dat_ns$Habitat
boxplot(y~g)
pairwise.wilcox.test(y, g, p.adj = "bonf") 
summary(aov(y~g))
anova(lm(y~ g, dat_ns))
TukeyHSD(aov(y~g))

y<-dat_ng$nav1
g<-dat_ng$Habitat
boxplot(y~g)
pairwise.wilcox.test(y, g, p.adj = "bonf") 
summary(aov(y~g))
anova(lm(y~ g, dat_ng))
TukeyHSD(aov(y~g))

with(dat_fc, tapply(FCover, list("hab#"=Habitat, "plot#"=plot), length1))
with(dat_th, tapply(tav, list("hab#"=Habitat, "plot#"=plot), length1))
describeBy(dat_th, dat_th$Habitat)

library(doBy)
res_fc<-summaryBy(FCover ~ Habitat + plot, data = dat_fc, 
          FUN = function(x) { c(m = mean(x), s = sd(x),n=length1(x)) } )
length(which(res_fc$Habitat=="Bush"))
sum(res_fc$FCover.n[which(res_fc$Habitat=="Bush")])
length(which(res_fc$Habitat=="Edge"))
sum(res_fc$FCover.n[which(res_fc$Habitat=="Edge")])
length(which(res_fc$Habitat=="Forest"))
sum(res_fc$FCover.n[which(res_fc$Habitat=="Forest")])
length(which(res_fc$Habitat=="Plantation"))
sum(res_fc$FCover.n[which(res_fc$Habitat=="Plantation")])
length(which(res_fc$Habitat=="Woodland"))
sum(res_fc$FCover.n[which(res_fc$Habitat=="Woodland")])
sum(res_fc$FCover.n)
length(res_fc$plot)

res_th<-summaryBy(tav ~ Habitat + plot, data = dat_th, 
                  FUN = function(x) { c(m = mean(x), s = sd(x),n=length1(x)) } )
length(which(res_th$Habitat=="Bush"))
sum(res_th$tav.n[which(res_th$Habitat=="Bush")])
length(which(res_th$Habitat=="Edge"))
sum(res_th$tav.n[which(res_th$Habitat=="Edge")])
length(which(res_th$Habitat=="Forest"))
sum(res_th$tav.n[which(res_th$Habitat=="Forest")])
length(which(res_th$Habitat=="Plantation"))
sum(res_th$tav.n[which(res_th$Habitat=="Plantation")])
length(which(res_th$Habitat=="Grassland"))
sum(res_th$tav.n[which(res_th$Habitat=="Grassland")])
sum(res_th$tav.n)
length(res_th$plot)

res_ng<-summaryBy(nav1 ~ Habitat + plot, data = dat_ng, 
                  FUN = function(x) { c(m = mean(x), s = sd(x),n=length1(x)) } )
length(which(res_ng$Habitat=="Bush"))
sum(res_ng$nav1.n[which(res_ng$Habitat=="Bush")])
length(which(res_ng$Habitat=="Edge"))
sum(res_ng$nav1.n[which(res_ng$Habitat=="Edge")])
length(which(res_ng$Habitat=="Forest"))
sum(res_ng$nav1.n[which(res_ng$Habitat=="Forest")])
length(which(res_ng$Habitat=="Plantation"))
sum(res_ng$nav1.n[which(res_ng$Habitat=="Plantation")])
length(which(res_ng$Habitat=="Grassland"))
sum(res_ng$nav1.n[which(res_ng$Habitat=="Grassland")])
sum(res_ng$nav1.n)
length(res_ng$plot)

res_ns<-summaryBy(nav_s ~ Habitat + plot, data = dat_ns, 
                  FUN = function(x) { c(m = mean(x), s = sd(x),n=length1(x)) } )
length(which(res_ns$Habitat=="Bush"))
sum(res_ns$nav_s.n[which(res_ns$Habitat=="Bush")])
length(which(res_ns$Habitat=="Edge"))
sum(res_ns$nav_s.n[which(res_ns$Habitat=="Edge")])
length(which(res_ns$Habitat=="Forest"))
sum(res_ns$nav_s.n[which(res_ns$Habitat=="Forest")])
length(which(res_ns$Habitat=="Plantation"))
sum(res_ns$nav_s.n[which(res_ns$Habitat=="Plantation")])

sum(res_ns$nav_s.n)
length(res_ns$plot)