#=============================================================================== # 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]