{
# read in the data file as a matrix of 13 columns
# first column is the year and the next 12 are the monthly
# rainfall values

#commands to read data, calculate basic monthly statistics, plot
#them for each month and also boxplots and histograms

#commands to perform robust stats - median, IQR, etc.

#commands to do a bivariate histogram


#read in data..

rain=matrix(scan("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

#calculate monthly mean/median, variance/iqr and skews

zmean=1:12

zvars=1:12

zskew=1:12

for(i in 1:12){
zmean[i]=mean(rain1[,i])

zmean[i]=quantile(rain1[,i],c(0.5))

zvars[i]=var(rain1[,i])
zvars[i]=diff(quantile(rain1[,i],c(0.25,0.75)))

zskew[i]=skew(rain1[,i])
}

#or
zmean=colSums(rain1)/N

#or
zmean=apply(rain1,2,mean)
zvars=apply(rain1,2,var)
zskew=apply(rain1,2,skew)

par(mfrow=c(3,1))

months=c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec")

plot(1:12,zmean,xlab="Months",ylab="Mean Rainfall", type="l", axes=F)
axis(1,at=1:12,labels=months)
axis(2)

plot(1:12,zvars,xlab="Months",ylab="Var Rainfall", type="l", axes=F)
xaxis(1,at=1:12,labels=months)
axis(2)

plot(1:12,zskew,xlab="Months",ylab="Skew of Rainfall", type="l", axes=F)
axis(1,at=1:12,labels=months)
axis(2)

#boxplot the raw data
par(mfrow=c(1,1))

rain2=as.data.frame(rain1)
xs=1:12
zz=myboxplot(split(t(rain2),xs),plot=F,cex=1.0)
zz$names=rep(" ",length(zz$names))
z1=bxp(zz,ylim=range(rain2,zmean),xlab="Months",ylab="Monthly rainfall",axes=F)
box(lty="solid",lwd=2)

axis(1,at=z1,labels=months)
axis(2)
lines(z1,zmean,lty=1,lwd=2)

#Plot histogram of each month and comment on the features..
par(mfrow=c(3,4))
for(i in 1:12){
hist(rain1[,i],xlab="Monthly rainfall", ylab="", probability=T, main="")
title(main=c(months[i],' rainfall'))
}

#to get wateryear from a matrix of values that corresponds to Jan ~ Dec..
zmat=cbind(rain1[1:N1,10:12],rain1[2:N,1:9])
wateryearsum=apply(zmat,1,sum)
plot(years[2:N],wateryearsum,xlab="Years",ylab="Water year rainfall", type="l")
abline(h=mean(wateryearsum))

#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)

nino3years=unique(test[,1])		#1900 to 2004

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

X=sesrain[nino3ses >= 0.75]
Y=sesrain[nino3ses <= -0.75]

myboxplot(X,Y,axes=F)
groups=c("El Nino years", "La Nina years")
axis(1,at=1:2,labels=groups)
axis(2)
box()

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

#bivariate histogram of NINO3 and AllIndia rainfall

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

 h2d <- hist2d(sesrain,nino3ses,show=FALSE, same.scale=TRUE, nbins=c(20,30))

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



###### histogram of the data..
nmons=12		#number of months..
xs=1:nmons
xsimpdf1=as.data.frame(t(rain1))
zz=myboxplot(split(xsimpdf1,xs),plot=F,cex=1.0)
zz$names=rep("",length(zz$names))
z1=bxp(zz,ylim=range(rain1,zmean),xlab="",ylab="",cex=1.25)
axis(1,at=xs,labels=months,cex=1.0)
lines(z1,zmean,lwd=2)


}