setwd("H:/BigLake/Paper/DataandCodeUpdatedSummer2021/Figure4")

data <- read.csv("epi_waterquality19952013.csv",header=TRUE)
data <- data[which(data$summer==1),]


######### 1995 to 2013 data #############

summersecchimean9513 <- aggregate(data[c("secchi")],
                                  list(data$lagoslakeid), FUN=mean, na.rm=TRUE)

summersecchisd9513 <- aggregate(data[c("secchi")],
                                  list(data$lagoslakeid), FUN=sd, na.rm=TRUE)

summersecchimean9513 <- as.data.frame(summersecchimean9513)
summersecchisd9513 <- as.data.frame(summersecchisd9513)

colnames(summersecchimean9513) <- c("id","mean")
colnames(summersecchisd9513) <- c("id","sd")

summersecchi9513 <- merge(summersecchimean9513,summersecchisd9513,by="id")

summersecchi9513<-summersecchi9513[!(summersecchi9513$mean=="NaN"),]

summersecchi9513$ngreater1 <- 0
summersecchi9513$ngreater1[summersecchi9513$sd>=0] <- 1

summersecchi9513$sd[is.na(summersecchi9513$sd)] <- 0

write.csv(summersecchi9513,"summersecchi9513.csv")

IDsinDataset <- read.csv("IDsinDataset.csv",header=TRUE)
colnames(IDsinDataset) <- c("id")
summersecchi9513dataset <- merge(summersecchi9513,IDsinDataset,by="id")


regions <- read.csv("RegionData.csv",header=TRUE)
summersecchi9513wregions <- merge(summersecchi9513dataset,regions,by="id")

mean.secchi = mean(summersecchi9513dataset$mean, na.rm = T)
median.secchi = median(summersecchi9513dataset$mean, na.rm = T)
var.secchi = var(summersecchi9513dataset$mean, na.rm = T)
sd.secchi = sd(summersecchi9513dataset$mean, na.rm = T)
min.secchi = min(summersecchi9513dataset$mean, na.rm = T)
max.secchi = max(summersecchi9513dataset$mean, na.rm = T)
num.obs = nrow(summersecchi9513dataset)
count.eutrophic = length(which(summersecchi9513dataset$mean<=1.83))
perceutrophic = count.eutrophic/num.obs
list(mean.secchi,median.secchi,sd.secchi,min.secchi,max.secchi,num.obs,
     count.eutrophic,perceutrophic)


HistSecchi <- summersecchi9513dataset$mean
HistSecchi[is.na(HistSecchi)] <- -1
library(plyr)
count(HistSecchi)

hist(HistSecchi, breaks=c(-1,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16),
     border="black", col="grey",xlab="Mean summer Secchi depth, 1995 and later",
     main=NA)

HistSecchi2 <- HistSecchi[which(HistSecchi>-1)]

HistSecchi3 <- summersecchi9513wregions[,c("mean","lowag","highag","other")]
HistSecchi3lowag <- HistSecchi3[which(HistSecchi3$lowag==1),1]
HistSecchi3highag <- HistSecchi3[which(HistSecchi3$highag==1),1]
HistSecchi3other <- HistSecchi3[which(HistSecchi3$other==1),1]

num.obs.lowag = length(HistSecchi3lowag)
count.eutrophiclowag = length(which(HistSecchi3lowag<=1.83))
perceutrophiclowag = count.eutrophiclowag/num.obs.lowag
list(num.obs.lowag,count.eutrophiclowag,perceutrophiclowag)

num.obs.highag = length(HistSecchi3highag)
count.eutrophichighag = length(which(HistSecchi3highag<=1.83))
perceutrophichighag = count.eutrophichighag/num.obs.highag
list(num.obs.highag,count.eutrophichighag,perceutrophichighag)

num.obs.other = length(HistSecchi3other)
count.eutrophicother = length(which(HistSecchi3other<=1.83))
perceutrophicother = count.eutrophicother/num.obs.other
list(num.obs.other,count.eutrophicother,perceutrophicother)

# xlab="Mean summer Secchi depth (1995 to 2013 measures)"

postscript(file="Figure4.eps",width=6, height=6)  

par(mfrow=c(2,2), las=2)
par(mai=c(0.5,0.4,0.1,0.1))

hist(HistSecchi2, breaks=seq(0,16,by=0.05),
     border="black", col="grey",main=NA)
abline(v = 1.83, col="red", lwd=3, lty=2)

hist(HistSecchi3lowag, breaks=seq(0,16,by=0.05), ylim=c(0,25),
     border="black", col="grey",main=NA)
abline(v = 1.83, col="red", lwd=3, lty=2)

hist(HistSecchi3highag, breaks=seq(0,16,by=0.05), ylim=c(0,40),
     border="black", col="grey",main=NA)
abline(v = 1.83, col="red", lwd=3, lty=2)

hist(HistSecchi3other, breaks=seq(0,16,by=0.05), ylim=c(0,150),
     border="black", col="grey",main=NA)
abline(v = 1.83, col="red", lwd=3, lty=2)

dev.off()



######### 2005 to 2013 data #############

data0513 <- data[which(data$sampleyear>=2005),]

summersecchimean0513 <- aggregate(data0513[c("secchi")],
                                  list(data0513$lagoslakeid), FUN=mean, na.rm=TRUE)

summersecchisd0513 <- aggregate(data0513[c("secchi")],
                                list(data0513$lagoslakeid), FUN=sd, na.rm=TRUE)

summersecchimean0513 <- as.data.frame(summersecchimean0513)
summersecchisd0513 <- as.data.frame(summersecchisd0513)

colnames(summersecchimean0513) <- c("id","mean")
colnames(summersecchisd0513) <- c("id","sd")

summersecchi0513 <- merge(summersecchimean0513,summersecchisd0513,by="id")

summersecchi0513<-summersecchi0513[!(summersecchi0513$mean=="NaN"),]

summersecchi0513$ngreater1 <- 0
summersecchi0513$ngreater1[summersecchi0513$sd>=0] <- 1

summersecchi0513$sd[is.na(summersecchi0513$sd)] <- 0

write.csv(summersecchi0513,"summersecchi0513.csv")

summersecchi0513wregions <- merge(summersecchi0513,regions,by="id")

mean.secchi = mean(summersecchi0513$mean, na.rm = T)
median.secchi = median(summersecchi0513$mean, na.rm = T)
var.secchi = var(summersecchi0513$mean, na.rm = T)
sd.secchi = sd(summersecchi0513$mean, na.rm = T)
min.secchi = min(summersecchi0513$mean, na.rm = T)
max.secchi = max(summersecchi0513$mean, na.rm = T)
num.obs = nrow(summersecchi0513)
list(mean.secchi,median.secchi,sd.secchi,min.secchi,max.secchi,num.obs)


HistSecchi <- summersecchi0513$mean
HistSecchi[is.na(HistSecchi)] <- -1
library(plyr)
count(HistSecchi)

hist(HistSecchi, breaks=c(-1,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16),
     border="black", col="grey",xlab="Mean summer Secchi depth, 1995 and later",
     main=NA)

HistSecchi2 <- HistSecchi[which(HistSecchi>-1)]

hist(HistSecchi2, breaks=seq(0,16,by=0.05),
     border="black", col="grey",xlab="Mean summer Secchi depth (based on measures from 2005 and later)",
     main=NA)
abline(v = 1.83, col="red", lwd=3, lty=2)

########### Compare two time periods ################

summersecchicompare <- merge(summersecchi9513,summersecchi0513,by="id")

library(Hmisc)
rcorr(summersecchicompare$mean.x, summersecchicompare$mean.y, type="spearman")
summary(lm(mean.x~mean.y, data=summersecchicompare))