# establish dirs
wd<-dirname(sys.frame(1)$ofile) # get dir of sourced R script and other dir
rdir<-"R Objects" # dir name for R objects within dir of R script
sdir<-"Data Tables" # dir name to store data tables within dir of R script
pdir<-"R Plots" # dir name to store R-plots within dir of R script

# load functions and data if needed
if (!exists("getHUC")) source(file.path(wd,"AqEECsFunctionsB.R")) # load functions
if (!exists("getHabitatHUC")) source(file.path(wd,"AqEECsFunctionsB.R")) # load functions
# if (!exists("TerrEECs.df")) load(file.path(wd,rdir,"TerrEECs.df")) # load data
if (!exists("overlaps.ar")) load(file.path(wd,rdir,"overlaps.ar")) # load data
if (!exists("cdl.names")) load(file.path(wd,rdir,"cdl.names")) # load data
if (!exists("huc_convert.df")) load(file.path(wd,rdir,"huc_convert.df")) # load data
if (!exists("huc_acres.df")) load(file.path(wd,rdir,"huc_acres.df")) # load data
if (!exists("speciesHUC.df")) load(file.path(wd,rdir,"speciesHUC.df")) # load data
if (!exists("habitatHUC.df")) load(file.path(wd,rdir,"habitatHUC.df")) # load data
if (!exists("species_info.df")) load(file.path(wd,rdir,"species_info.df")) # load data

TerrEECs.df<-read.csv(file.path(wd,"TerrestrialEECs.csv"))
names(TerrEECs.df)[1]<-"Pesticide"

chemlist<-c("chlorpyrifos","diazinon","malathion","prometryn","bromoxynil octanoate")
huclist<-c("1","2","3","4","5","6","7","8","9","10a","10b","11a","11b","12a","12b","13","14","15a","15b","16a","16b","17a","17b","18a","18b","19a","19b","20a","20b","21")
binlist<-c("0","2","3","4","5","6","7")
readoptions<-c("Species range","Read from file","Choose HUC2s")

# select data to summarize ***********  USER selection section of inputs
usedialog<-FALSE # ********************* use dialog boxes to enter data, leave FALSE for now
chemtext<-chemlist[5] # *************** select chemical from chemlist above e.g. 1 for cyp, 2 for dia, 3 for ma
toxfilename<-"bro.terr.tox.txt" # file with info for toxicity data section of the plot
readhucs<-readoptions[1] # ******************* how to get HUC12 and HUC2 lists, from readoptions list e.g. 1 for "Species range"
speciesID<-5180 # ************* Select species Entity ID
plotHabitat<-TRUE # set to true to plot critical habitat rather than sepcies range
useHUC12s<-T # ******************* use HUC12 info otherwise rely on HUC2 and don't plot overlap
useSpeciesBins<-F # **************** plot all the bins assigned to a species or specific bins
OmitBin3<-T  # **************** don't plot bin 3 values for individual uses
OmitBin4<-T  # **************** don't plot bin 4 values for individual uses
bintext<-c("2","5") # *************** bins to plot if not based on species
huctext<-c("18a") # *************** list of HUC2s for EEC data if not specified by the species HUC12 list
avgper<-c(1,2,3) # ************** EEC averaging periods to plot, 1-4 correspond to 1-d, 4-d, 21-d, 60-d 
multiHUC2adj<-TRUE # if multiple HUC2s adjust EECs based on proportion of total area
# *************************************

toxdatafile<-file.path(wd,toxfilename)
ESUfile<-file.path(wd,"SalmonESUs.csv")
ESUnames<-read.csv(ESUfile)

if (usedialog) {
	chemtext<-select.list(chemlist,preselect=chemtext,title="Select pesticide")
	toxdatafile<-file.choose()
	bintext<-select.list(binlist,preselect=bintext,title="Select desired bins",multiple=TRUE)
	readhucs<-select.list(readoptions,preselect=readhucs,title="Select HUCs")
	toxfilename<-basename(toxdatafile)
}
if (readhucs=="Read from file") {
	foo<-read.delim(file.choose()) # use file.choose() for windows tk_choose.files() for Mac
	listname<-select.list(names(foo))
	speciesID<-paste0(listname," (List)")
	HUCs<-sort(na.omit(foo[,grep(listname,names(foo))])) # make sure list is sorted in HUC12 order
	huc2list<-subset(huc_convert.df, HUC12 %in% HUCs)
	huc2list<-huc2list[order(huc2list$HUC12),]
	huctext<-unique(as.character(huc2list$HUC02))
	useHUC12s<-TRUE
	if (length(HUCs)==0) {
		stop("No HUC12s in list")
	}
}
if (readhucs=="Species range") {
	if (useSpeciesBins) {
		bintext<-as.character(unique(species_info.df[species_info.df$EntityID==speciesID,]$Bin))	
	}
	if (plotHabitat==FALSE) {
		HUCs<-sort(getHUC(speciesID)) # species entityID to get list of HUCs in order
		HUCtype<-"Range"
	} else {
		HUCs<-sort(getHabitatHUC(speciesID)) # species entityID to get list of HUCs in order
		HUCtype<-"Habitat"
	}	
	if (length(ESUnames[ESUnames$ID==speciesID,6])>0) {
		speciesLabel<-ESUnames[ESUnames$ID==speciesID,6]
	} else {
		speciesLabel<-species_info.df[species_info.df$EntityID==speciesID,][1,5]
	}
	speciesID<-paste0(speciesLabel," (",speciesID,") ",HUCtype)
	huc2list<-subset(huc_convert.df, HUC12 %in% HUCs)
	huc2list<-huc2list[order(huc2list$HUC12),]
	huctext<-unique(as.character(huc2list$HUC02))
	useHUC12s<-TRUE
	if (length(HUCs)==0) {
		stop("No HUC12s in list")
	}
}
if (readhucs=="Choose HUC2s") {
	if (usedialog) {
		huctext<-select.list(huclist,preselect=huctext,title="Select desired HUC2s",multiple=TRUE)		
	}
	if (useSpeciesBins) {
		bintext<-as.character(unique(species_info.df[species_info.df$EntityID==speciesID,]$Bin))	
	}
	speciesID<-paste0(huctext,collapse="_")
	useHUC12s<-FALSE
}

toxdata.df<-read.delim(toxdatafile) # read tox data from tab-delimited file, columns in pairs with conc col and label col
numToxRows<-dim(toxdata.df)[2]/2
minTox<-min(toxdata.df[,c(1,3)],na.rm=TRUE)
maxTox<-max(toxdata.df[,c(1,3)],na.rm=TRUE)

# *************************************
# get use overlap data based on HUC12s or set to dummy fill data
GroupedHUCdata<-data.frame(cdl.names[3:27])
if (useHUC12s) {
	tempHUCs<-subset(huc_acres.df, HUC12 %in% HUCs) # get list of HUC12s with acres
	tempHUCs<-tempHUCs[order(tempHUCs$HUC12),] # order list by HUC12 to make sure and match order of HUC12s in other data objects
	HUCsTotalAcres<-tempHUCs[,2]
	HUCsRows<-as.numeric(row.names(tempHUCs))
	HUCsPercents<-overlaps.ar[,HUCsRows,3:27] # overlaps omit HUC12 and acres
	
	GroupedUseAcres<-array(0,dim=c(6,25)) # array to collect acres for each HUC12 based on use for each of 6 years
	if (length(HUCsRows)==1) {
		GroupedUseAcres<-HUCsPercents*HUCsTotalAcres/100 # acres foe each of 6 years for each use
	} else {
		for (i in 1:length(HUCsRows)) {
			GroupedUseAcres<-GroupedUseAcres + HUCsPercents[,i,]*HUCsTotalAcres[i]/100 # add acres for each of 6 years by use for the ith HUC in the list
		}
	}
	
	GroupedUsePercents<-100*GroupedUseAcres/sum(HUCsTotalAcres)
	GroupedHUCdata<-cbind(GroupedHUCdata,apply(GroupedUsePercents,2,median),apply(GroupedUsePercents,2,min),apply(GroupedUsePercents,2,max))	
} else {
	GroupedHUCdata<-cbind(GroupedHUCdata,array(-1,dim=25),array(-1,dim=25),array(-1,dim=25)) # dummy data to fill object	
	HUCsRows<-0
	HUCsTotalAcres<-0
}
names(GroupedHUCdata)<-c("scenario","PercentUse","min","max")

# ************************************
# get EEC info using TerrEECs instead
pks.df<-TerrEECs.df[TerrEECs.df$Pesticide==chemtext,-1]
scentext<- levels(pks.df$Scenario)
croptext<- levels(pks.df$Crop) 

temp.df<-subset(pks.df,HUC2 %in% huctext & Bin %in% bintext)
if (dim(temp.df[temp.df$Scenario=="Other Rowcrops",])[1]>0) {
	temp.df[temp.df$Scenario=="Other Rowcrops",]$Scenario<-"Other RowCrops" # make Scenario names consistent
}
rm(pks.df)

temp.df<-temp.df[temp.df$X1.day>0,] # omit rows with 0 EEC

# check for multiple huc2 EECs
# if desired pad EECs from each HUC2 to reflect contribution to overall distribution proportional to area
if (multiHUC2adj & length(huctext)>1 & useHUC12s) {
	hucsWithAcres<-merge(huc2list,tempHUCs)
	huc2Acres<-aggregate(hucsWithAcres[,3],by=list(hucsWithAcres$HUC02),FUN="sum")
	huc2Acres<-cbind(huc2Acres,round(100*huc2Acres$x/sum(huc2Acres$x)))
	names(huc2Acres)<-c("HUC2","Acres","Percent")
	for (i in 1:dim(huc2Acres)[1]) {
		t.df<-temp.df[as.character(temp.df$HUC2)==as.character(huc2Acres$HUC2[i]),]
		for (j in 1:(huc2Acres$Percent[i]-1)) {
			temp.df<-rbind(temp.df,t.df)
		}
	}
}

temp_medians<-aggregate(temp.df[,8:11],by=list(temp.df$Bin,temp.df$Scenario,temp.df$Crop),FUN="median")
temp_quants<-aggregate(temp.df[,8:11],by=list(temp.df$Bin,temp.df$Scenario,temp.df$Crop),FUN="quant95")

# **********************************
# combined plot of use and EECs
temp_y<-array(-1,dim=length(temp_medians$Group.2)) # vector for the median percent use data over 6 years 
temp_min<-array(-1,dim=length(temp_medians$Group.2)) # vector for the min percent use data over 6 years
temp_max<-array(-1,dim=length(temp_medians$Group.2)) # vector for the max percent use data over 6 years
if (useHUC12s) {
	for (i in 1:length(temp_y)) {
		temp_y[i]<-sum(na.omit(GroupedHUCdata[GroupedHUCdata$scenario==as.character(temp_medians$Group.2[i]),])[,2])
		temp_min[i]<-sum(na.omit(GroupedHUCdata[GroupedHUCdata$scenario==as.character(temp_medians$Group.2[i]),])[,3])
		temp_max[i]<-sum(na.omit(GroupedHUCdata[GroupedHUCdata$scenario==as.character(temp_medians$Group.2[i]),])[,4])
	}
}

temp_data<-cbind(temp_medians,temp_quants[4:7]) # combine median and range EECs
temp_data<-cbind(temp_data,temp_y,temp_min,temp_max) # combine percent areas

if (useHUC12s & max(temp_y)>0) {
	temp_data<-temp_data[temp_data$temp_y>0,] # remove uses with zero percent area
	temp_order<-as.numeric(factor(temp_data$temp_y)) # establish order based on percent area
} else {
	temp_order<-as.numeric(factor(temp_data$Group.2)) # establish order based on Scenario
}
temp_data<-cbind(temp_data,temp_order) # add order for y plotting of EEC data

maxed_uses<-unique(temp_data[temp_data$temp_order==max(temp_data$temp_order),]$Group.2) # identify max uses (multiple could be at 100%)
for (i in 1:length(maxed_uses)) {
	temp_data[temp_data$Group.2==maxed_uses[i],]$temp_order<-temp_data[temp_data$Group.2==maxed_uses[i],]$temp_order + (i-1) # increment order to separate uses
}

miny<-min(temp_data$temp_y)
min1<-min(temp_data[,8:10]) # min of quant data omitting 60-d
max1<-max(temp_data[,8:10]) # max of quant data omitting 60-d

# if (maxTox>min1) {min1<-minTox}
# if (minTox>max1) {max1<-maxTox}

min1<-min(c(minTox,min1))
max1<-max(c(maxTox,max1))

# min1<-0.1 # overide auto scale if desired

if (useHUC12s & max(temp_y)>0) {
	# omit bin 3 data except for 100% Uses unless no data remains afterward
	if (OmitBin3 & dim(temp_data[!(temp_data$Group.1==3 & temp_data$temp_y<99),])[1]>0) {
		temp_data<-temp_data[!(temp_data$Group.1==3 & temp_data$temp_y<99),]	
	}
	# omit bin 4 data except for 100% Uses unless no data remains afterward
	if (OmitBin4 & dim(temp_data[!(temp_data$Group.1==4 & temp_data$temp_y<99),])[1]>0) {
		temp_data<-temp_data[!(temp_data$Group.1==4 & temp_data$temp_y<99),]	
	}
}

# **********************************************************
# plot data with uniform spacing along Y axis
# plot median, 5%, 95% on continuous Y axis

savePlot<-menu(c("Save","View"),graphics=TRUE,title="Save or View Plot?")

if (savePlot==1) {
	plotName<-file.path(wd,pdir,paste0(speciesID,"_",chemtext,"_Terrestrial.pdf"))
	pdf(file=plotName,width=6.5,height=5.5,pointsize=9)
	par(mar=c(8.1,16.5,0.5,9),cex.axis=1, bty="n")
} else {
	par(mar=c(8.1,16.5,0.5,10),cex.axis=1, bty="n")
}

maxy<-max(temp_data$temp_order) # number of uses to plot need to add more for bins 3 and 4 and tox data (2 plus 1 gap plus numToxRows)
ypad<-1 # use 3 to make room for 2 rows of bins 3 and 4	
plot(temp_data$X1.day, temp_data$temp_order, type="n", ylim = c(0, maxy+ypad+0.3+numToxRows), xlim = c(min1,max1), log = "x", xlab ="", yaxt="n", ylab="")

for (i in seq(1,maxy+ypad-1,by=2)) {
#	abline(h=i-0.35,lty="dotted")
	rect(min1,i-0.35,max1,1+i-0.35,border="transparent",col="grey90")
}

binmark<-c(21,21,21,21,21,23,23) # list of symbols for bins 0,2-7
bincol<-c("black","white","grey","black","black","grey","black") # list of colors for bin 0,2-7 symbols

avgperlabels<-c("AgDrift","TerrPlant-Dry","TerrPlant-Wet","60-d")
avgperinfo<-paste0("", paste0(avgperlabels[avgper],collapse=", "), " (from bottom up)")
if (useHUC12s) {
	HUClabel<-paste0(speciesID,": ",length(HUCsRows), " HUC12s")
} else {
	HUClabel<-"HUC2 info only"
}

titletext<-c(chemtext,avgperinfo,paste0("huc:",paste0(huctext,collapse=",")),HUClabel)

# add averaging periods 4:7 for 1, 4, 21, 60 d
for (i in avgper+3) {
	arrows(unlist(temp_data[,i+4])[,1], temp_data$temp_order + (i-4)/10, unlist(temp_data[,i+4])[,2], temp_data$temp_order + (i-4)/10, length = 0, angle = 90, code = 3)
	for (j in c(7:2,0)) {
		points(temp_data[temp_data$Group.1==j,i], temp_data[temp_data$Group.1==j,]$temp_order + (i-4)/10, pch=binmark[grep(j,c(0,2:7))], bg=bincol[grep(j,c(0,2:7))])
	}
}

# add levels "Bin 3" and "Bin 4" to temp_data columns 2 and 3
levels(temp_data[,2])<-c(levels(temp_data[,2]),"Bin 3","Bin 4")
levels(temp_data[,3])<-c(levels(temp_data[,3]),"Bin 3","Bin 4")

# add aggregate data for bins 3 and 4 if in bintext and overlap info is available
if (useHUC12s & max(temp_y)>0) {
	t.bins<-as.numeric(bintext[bintext %in% c("3","4")]) # list bins 3 and/or 4 if in bintext
	if (length(t.bins>0)) {
		for (t.b in t.bins) {
			t.df<-aggregate(temp_medians[temp_medians$Group.1==t.b,4:7],by=list(temp_medians[temp_medians$Group.1==t.b,]$Group.2),FUN="max") # maximum run for each use and bin 3
			names(t.df)[1]<-"scenario"
			
			t.df<-merge(t.df,GroupedHUCdata)
			t.df<-t.df[t.df$PercentUse<99,]
			
			# get row names of temp_medians that contain the data in t.df
			t.rows<-as.numeric(row.names(temp_medians[temp_medians$Group.1==t.b & temp_medians$X60.day==t.df$X60.day[1] & temp_medians$Group.2==t.df$scenario[1],]))[1] # start with first row of t.df (pick only one row if dup X60day eecs)
			for (i in 2:dim(t.df)[1]) {
				t.rows<-c(t.rows,as.numeric(row.names(temp_medians[temp_medians$Group.1==t.b & temp_medians$X60.day==t.df$X60.day[i] & temp_medians$Group.2==t.df$scenario[i],]))[1]) # pick only one row if dup 60-d eecs
			}
			
			tq.df<-temp_quants[t.rows,] # get quantile data from same runs
			t.df<-cbind(t.df,tq.df[,4:7])
	
			temp_data<-rbind(temp_data,temp_data[1,]) # append row to temp_data
			newrow<-dim(temp_data)[1] # row number of row added
			# update new row in temp_data for bin 3 or 4 info
			temp_data[newrow,1]<-t.b
			temp_data[newrow,2:3]<-paste0("Bin ",t.b)
			temp_data[newrow,12:14]<-100
			temp_data[newrow,15]<-maxy + t.b - 2
			
			# calculate bin 3 and 4 averaging periods 4:7 for 1, 4, 21, 60 d
			for (i in avgper+3) {
				t.median<-sum(t.df[,i-2]*t.df[,6]/100)
				t.max<-sum(unlist(t.df[,5+i])[,2]*t.df[,6]/100)
				t.min<-sum(unlist(t.df[,5+i])[,1]*t.df[,6]/100)

				# update new row in temp_data 
				temp_data[newrow,i]<-t.median
				temp_data[newrow,i+4][1]<-t.min
				temp_data[newrow,i+4][2]<-t.max

				# plot handling off scale points with < or >
				if (t.max<min1) {
					points(min1,maxy + t.b - 2 + (i-4)/10,pch=60)
				} else {
					if (t.min>max1) {
						points(max1,maxy + t.b - 2 + (i-4)/10,pch=62)				
					} else {
						arrows(t.min, maxy + t.b - 2 + (i-4)/10, t.max, maxy + t.b - 2 + (i-4)/10, length = 0, angle = 90, code = 3)
						points(t.median, maxy + t.b - 2 + (i-4)/10, pch=binmark[t.b-1], bg=bincol[t.b-1])
					}
				}
			}
		}
	}	
}

# add various labels
axis(2,at=c(0,temp_data$temp_order,maxy+seq(1:(ypad+numToxRows))),labels=c(paste0("All HUCs (",format(sum(HUCsTotalAcres),digits=1,nsmall=0,trim=TRUE,scientific=FALSE)," acres)"),as.character(paste0(temp_data$Group.2," (",format(sum(HUCsTotalAcres)*temp_data$temp_y/100,digits=1,nsmall=0,trim=TRUE,scientific=FALSE),")")),rep("",ypad-1),"--------------",names(toxdata.df[seq(1,numToxRows*2-1,by=2)])),las=1,padj=0.4,lwd=0) # addition of Bin 3 and 4 and tox data

if (useHUC12s & max(temp_y)>0) {
	l.med<-array(data="",dim=length(temp_data$temp_y))
	l.min<-array(data="",dim=length(temp_data$temp_y))
	l.max<-array(data="",dim=length(temp_data$temp_y))
	prec<-2
	minprec<-10^(-1*prec)
# round percent areas
	for (i in 1:length(temp_data$temp_y)) {
		if (temp_data$temp_y[i]<minprec) {
			l.med[i]<-paste0("<",minprec)
		} else {
			l.med[i]<-round(temp_data$temp_y[i],digits=prec)
		}	
		if (temp_data$temp_min[i]<minprec) {
			l.min[i]<-paste0("<",minprec)
		} else {
			l.min[i]<-round(temp_data$temp_min[i],digits=prec)
		}	
		if (temp_data$temp_max[i]<minprec) {
			l.max[i]<-paste0("<",minprec)
		} else {
			l.max[i]<-round(temp_data$temp_max[i],digits=prec)
		}	
	}
	axis(4,at=temp_data$temp_order,labels=paste0(l.med," % (",l.min,", ",l.max,")"),las=1,lwd=0) # percent area on right Y-axis
#	mtext("Median percent of total acres", side = 4, line=13, adj=0, at=1)
#	mtext("(over 6 years with min, max)", side = 4, line=14, adj=0, at=1)	
}

# mtext("Uses (median total acres over 6 yrs)", side = 2, line=15, adj=0, at=1)
mtext("Uses", side = 2, line=15, adj=0, at=1)
mtext("Toxicity Data", side = 2, line=15, adj=1)

mtext("Terrestrial EEC (lbs a.i./A)", side=1, line=2)
mtext(paste0(titletext[1],": ",toxfilename),side=1, line=4)
mtext(titletext[2],side=1, line=5)
mtext(titletext[3],side=1, line=6)
mtext(titletext[4],side=1, line=7)

# add additional data to show magnitude of effect info
for (i in seq(1,numToxRows*2-1,by=2)) {
	if (!is.na(toxdata.df[1,i])) {
		if (max(na.omit(toxdata.df[,i]))<min1) {
			points(min1,maxy+ypad+(i+1)/2,pch=60)
		} else {
			if (min(na.omit(toxdata.df[,i]))>max1) {
				points(max1,maxy+ypad+(i+1)/2,pch=62)				
			} else {
				points(toxdata.df[,i],rep(maxy+ypad+(i+1)/2,length(toxdata.df[,i])),pch=3)
				arrows(min(na.omit(toxdata.df[,i])),maxy+ypad+(i+1)/2,max(na.omit(toxdata.df[,i])),maxy+ypad+(i+1)/2,length=0,angle=90,code=3)
				text(toxdata.df[,i], rep(maxy+ypad+(i+1)/2+0.2,length(toxdata.df[,i])), toxdata.df[,i+1], adj=c(0.5,1), cex=0.6)
			}
		}
	}
}

# add legend
legend(x="bottomleft",pch=binmark[c(0,2:7) %in% bintext],pt.bg=bincol[c(0,2:7) %in% bintext],bty="n",horiz=TRUE,legend=c("G","A"),cex=0.7)

# add X limit line
# abline(v=2000, lty=2)

if (savePlot==1) {
	dev.off()
}

# save table of EECs to file if desired
if (menu(c("Yes","No"),graphics=TRUE,title="Save Table?")==1) {
	tableName<-file.path(wd,sdir,paste0(speciesID,"_",chemtext,"_Terrestrial.csv"))
	temp_table<-temp_data
	names(temp_table)[1:3]<-c("Bin","Use","Crop")
	write.csv(temp_table,file=tableName,row.names=FALSE)
}

