{
# Plot a 2-D histogram
# Get the All Indian Seasonal Rainfall
# Get the seasonal NINO3

#commands to do a bivariate histogram


#read in data..

rain=matrix(scan("http://civil.colorado.edu/~balajir/r-session-files/aismr.txt"),ncol=13,byrow=T)

#pick the monthly rainfall values excluding the year.
rain1=rain[,2:13]
years=rain[,1]		#1871 - 1999		
N=length(years)		#of years
N1=N-1

#read the NINO3 index..
test=matrix(scan("nino3-index.txt"),ncol=3,byrow=T)
nino3=matrix(test[,3],ncol=12,byrow=T)

#get the nino3 of the summer (June - Sep) season
nino3ses=apply(nino3[,6:9],1,mean)

#get the seasonal rainfall for the common years - i.e. 1900 - 1999
N=1999-1900+1

N2=1900-1871+1
N3=N2+N-1

sesrain=apply(rain1[,6:9],1,mean)
sesrain=sesrain[N2:N3]

nino3ses=nino3ses[1:N]

#Both seasonal rainfall and seasonal NINO3 are now of the same length..

#############################

#bivariate histogram of NINO3 and AllIndia rainfall

library(gplots)
library(akima)          #perhaps not required wth gplots..

ncls = round(log(length(sesrain), 2) + 1)
h2d <- hist2d(sesrain,nino3ses,show=FALSE, same.scale=FALSE, nbins=c(ncls,ncls))
#h2d <- hist2d(sesrain,nino3ses,show=FALSE, same.scale=FALSE, nbins=c(15,15))

# You can change the values in nbins and also other arguements in the above
# command - Please see the help and play around.

#Draw a perspective plot
persp( h2d$x, h2d$y, h2d$counts, ticktype="detailed", theta=30, phi=30, expand=0.5, 
shade=0.5, col="cyan", ltheta=-30)

# the arugments in the above command can be changed to change the viewing
# angle of the plot and also other aspects - color, mesh etc. 

}