####################################
############Function to examine trends in GW level data
Wells_trends=function(nrecords=20,file_input_path="C:\\dave\\scratch\\GWUseCase\\default.dbf")
{
###function requires package foreign
library(foreign)
########################function to plot trends in water level using the equation
########################Yk(t)=Xk(t)-Xk(tk)+1/(k-1)SUM{(from i to k-1) Yi(tk)}
########################where ti  equals first time for series i
########################default.dbf contains sites Feature ID's (HydroID's), time and values 
###########################nreocrds is the threshold needed for selecting wells
##########################reading the data base file
raw_data<-read.dbf(file=file_input_path)
##########################getting the number of wells 
wells_feature<-unique(raw_data$FeatureID)
###########################getting total number of wells
nwells<-length(wells_feature)
################################fixing the data
##################################extracting the data with GW level records >= nrecords and sorting them based on time 
TSValue<-list()
Time<-list()
well_no<-list()
i=0
for(j in 1:nwells)
	{
		if(length(raw_data$TSDateTime[raw_data$FeatureID==wells_feature[j]])>=nrecords)
			{
			well_no[[i+1]]<-wells_feature[j]
			itt<-sort(unclass(raw_data$TSDateTime[raw_data$FeatureID==wells_feature[j]]),method="sh",index.return=TRUE)$ix#sorting the records based on time
			TSValue[[i+1]]<-raw_data$TSValue[raw_data$FeatureID==wells_feature[j]][itt]
			Time[[i+1]]<-raw_data$TSDateTime[raw_data$FeatureID==wells_feature[j]][itt]
			i=i+1
			}
	}

#################################The number of wells which satisfies our criteria
nwells<-dim(summary(well_no))[1]

#################################getting the date of each well in conversion format
first_date_wells<-rep(NA,nwells)
	for(i in 1:nwells)
		{
			first_date_wells[i]<-unclass(Time[[i]][1])
		}
################################sorting these dates
itt2<-sort(first_date_wells,method="sh",index.return=TRUE)$ix

#################################getting the order of wells based upon their start time 
wells_ID_ordered<-well_no[itt2]

#################################the list order based upon the start date of all wells is found as itt2
##################################let's get the scaled levels that is X(t)-X(t1) where t1 represents the first time record in each well records
scaled_levels<-list()
	for(i in 1:nwells)
		{
			list<-itt2[i]
			scaled_levels[[i]]<-TSValue[[list]]-TSValue[[list]][1]	
		}
#################################let's get the time associated with scaled_levels 
Time_for_scaled_levels<-list()
	for(i in 1:nwells)
		{
			list<-itt2[i]		
			Time_for_scaled_levels[[i]]<-Time[[list]]
		}

################################let's do the correction
	for(i in 2:nwells)
		{
#################################first well remains the same
it_index<-rep(NA,i-1)
			for(j in 1:(length(it_index)))
				{
					it_index[j]<-findInterval(unclass(Time_for_scaled_levels[[i]]),unclass(Time_for_scaled_levels[[i-j]]))[1]
				} 
			uu<-rep(NA,2*length(it_index))
			temp<-rep(NA,2*length(it_index))
			temp_wells<-rep(NA,length(temp)/2)	
				for(l in 1:(length(it_index)))
					{
						uu[(1+(l-1)*2):(2*l)]<-c(it_index[l],it_index[l]+1)
					}
				for(k in 1:(length(it_index)))
					{
						temp[(1+(k-1)*2):(2*k)]<-c(scaled_levels[[(i-k)]][uu[k+(k-1)]],scaled_levels[[(i-k)]][uu[2*k]])
						temp_wells[k]<-mean(temp[(1+(k-1)*2):(2*k)],na.rm=TRUE)
					}
			scaled_levels[[i]]<-scaled_levels[[i]]+mean(temp_wells,na.rm=TRUE)
		}


###############################getting the min and max level record for plotting as well as start and end date 
maxy<-rep(NA,nwells)
miny<-rep(NA,nwells)
	for(i in 1:nwells)
		{
			maxy[i]<-max(-scaled_levels[[i]])
			miny[i]<-min(-scaled_levels[[i]])
		}
maxy<-ceiling(max(maxy))
miny<-ceiling(min(miny))
###############################making start date as 30 days before the beginning of first series (graphing feature)
start_date<-Time_for_scaled_levels[[1]][1]-30
################################making end date as 30 days after the end of last series (graphing feature)
pp<-rep(NA,nwells)
	for(i in 1:nwells)
		{
			pp[i]<-unclass(Time_for_scaled_levels[[i]][length(Time_for_scaled_levels[[i]])])
		}
end_date<-as.Date(max(pp), origin="1970-01-01")+30
list(nwells=nwells,levels=scaled_levels,time=Time_for_scaled_levels,start_date=start_date,end_date=end_date,maxy=maxy,miny=miny)
###################################end
}


##################################Plotting
Global_data_Upstream<-Wells_trends(20,"C:\\Watershed_Model\\GroundWater\\wells\\Exercise_Weber_Wells\\GWUseCase\\UpStreamWells_records.dbf")
win.graph()
y_scale<-c(pretty(Global_data_Upstream$maxy:Global_data_Upstream$miny))
plot(Global_data_Upstream$time[[1]],-Global_data_Upstream$levels[[1]],col="black",type="l",ylab="Adjusted Well Level (ft)",xlim=c(as.Date(Global_data_Upstream$start_date),as.Date(Global_data_Upstream$end_date)),yaxt="n",ylim=c(y_scale[1],y_scale[length(y_scale)]),lty=1) 
axis(2,-pretty(Global_data_Upstream$maxy:Global_data_Upstream$miny),pretty(Global_data_Upstream$maxy:Global_data_Upstream$miny))
for(i in 2:Global_data_Upstream$nwells)
	{
		lines(Global_data_Upstream$time[[i]],-Global_data_Upstream$levels[[i]],col="black",type="l",lty=i)
	}

title("Adjuested groundwater level for Upstream wells",font=2)
mtext("90 wells",side=3,line=0,adj=1,font=4)

Global_data_Midstream<-Wells_trends(20,"C:\\Watershed_Model\\GroundWater\\wells\\Exercise_Weber_Wells\\GWUseCase\\MidStreamWells_records.dbf")
win.graph()
y_scale<-c(pretty(Global_data_Midstream$maxy:Global_data_Midstream$miny))
plot(Global_data_Midstream$time[[1]],-Global_data_Midstream$levels[[1]],col="red",type="l",ylab="Adjusted Well Level (ft)",xlim=c(as.Date(Global_data_Midstream$start_date),as.Date(Global_data_Midstream$end_date)),yaxt="n",ylim=c(y_scale[1],y_scale[length(y_scale)]),lty=1) 
axis(2,-pretty(Global_data_Midstream$maxy:Global_data_Midstream$miny),pretty(Global_data_Midstream$maxy:Global_data_Midstream$miny))
for(i in 2:Global_data_Midstream$nwells)
	{
		lines(Global_data_Midstream$time[[i]],-Global_data_Midstream$levels[[i]],col="red",type="l",lty=i)
	}

title("Adjuested groundwater level for Midstream wells",font=2)
mtext("16 wells",side=3,line=0,adj=1,font=4)

Global_data_Downstream<-Wells_trends(20,"C:\\Watershed_Model\\GroundWater\\wells\\Exercise_Weber_Wells\\GWUseCase\\DownStreamWells_records.dbf")
win.graph()
y_scale<-c(pretty(Global_data_Downstream$maxy:Global_data_Downstream$miny))
plot(Global_data_Downstream$time[[1]],-Global_data_Downstream$levels[[1]],col="blue",type="l",ylab="Adjusted Well Level (ft)",xlim=c(as.Date(Global_data_Downstream$start_date),as.Date(Global_data_Downstream$end_date)),yaxt="n",ylim=c(y_scale[1],y_scale[length(y_scale)]),lty=1) 
axis(2,-pretty(Global_data_Downstream$maxy:Global_data_Downstream$miny),pretty(Global_data_Downstream$maxy:Global_data_Downstream$miny))
for(i in 2:Global_data_Downstream$nwells)
	{
		lines(Global_data_Downstream$time[[i]],-Global_data_Downstream$levels[[i]],col="blue",type="l",lty=i)
	}

title("Adjuested groundwater level for Downstream wells",font=2)
mtext("58 wells",side=3,line=0,adj=1,font=4)




