#===============================================================================
# Name : polyplot
# Original author: Adam Wilson, with inspiration from Paul Murrell's book "R Graphics" (Especially chapter 3 and figure 3.25)
# Date (dd/mm/yy): 3/5/08
# Aim : Produce a line plot with filled regions (polygons) between the line and some horizontal line, such as zero.
#===============================================================================
## This plot will probably be most useful for timeseries where the user hopes to emphasize the value with respect to some constant (such as zero or a mean)
##Define Function
polyplot=function(x,y,hline=median(y),col.high="red",col.low="blue",NAValue=-9999,...){
## hline enter a number or function (such as "mean(y)")
## col.high color for area above hline
## col.low color for area below hline
## NAValue used as a placeholder, select another value if your data contains a value of -9999 or there will be problems with the polygons
## ... other terms passed to plot()

##define vectors to be used
n=length(x)
int=vector(mode="numeric", length=n-1)
int2=vector(mode="numeric", length=n-1)
int3=vector(mode="numeric", length=n-1)
##computes intercepts of line with hline
for(i in 1:n-1){
ii=c(i,i+1)
int[i]=((hline-y[i])/lm(y[ii]~x[ii])$coefficient[2])+x[i]
int2[i]=ifelse(int[i]>=ii[1],int[i],NA)
int3[i]=ifelse(int2[i]<=ii[2],int2[i],NA)
}
## rearrange the vectors to the correct form for the polygon function. There
## may be a better way to do this, but it works....
int3=c(int3,NA)
int3[is.na(int3)]=NAValue
yint=ifelse(int3==NAValue,NAValue,hline)
newdata=(cbind(x,int3,y,yint,NA))
newdata[int3==NAValue,5] = NAValue
xnew=as.vector(t(cbind(x,int3,newdata[,5],int3)))
ynew=as.vector(t(cbind(y,yint,newdata[,5],yint)))
xnew=xnew[xnew!=NAValue]
ynew=ynew[ynew!=NAValue]
## draw the plots
plot(ynew~xnew,type="n", axes=T, ann=T,...)
n=length(xnew)
col=ifelse(rep(ynew[1]<hline,2),c(col.low,col.high),c(col.high,col.low))
polygon(c(xnew[1],xnew,xnew[n]),c(hline,ynew,hline),col=col,border=NA)
abline(h=hline,col="grey")
lines(x,y)
}

polyplot(1:50,rnorm(50),las=1,cex.lab=2,cex.axis=2)