##PLOT OF TRENDS
   #REQUIRES PACKAGES ncdf

    ##LOAD DATA SET

	library(ncdf)
	url<-"d:/climate/data/jones"
 	loc<-file.path(url,"HadCRUT3.nc" ) #[1] "c"
 	v<-open.ncdf(loc)
		#test <- v$var[[1]]
		#test$size # 72   36    1 1883
		#test$dim[[1]]$vals #[1] -177.5 -172.5 -167.5 -162.5 -157.5 -152.5 -147
		#test$dim[[2]]$vals # [1] -87.5 -82.5 -77.5 -72.5 -67.5 -62.5 -57.5

	instr <- get.var.ncdf( v, v$var[[1]]) # 1856 2005
	dim(instr)# [1]   72   36 1883
		#this is organized in 72 longitudes from -177.5 to 177.5 and 36 latitudes from -87.5 to 87.5


##MAKE ARRAY 
     #I did a couple of other calculations and have a larger array here 
	stat<-array(NA,dim=c(72,36,4));	##LS trend
	year<-(1:1883)/12
	index<- (12*(1978-1849)-1):(12*(2006-1849)-4)

    for (k in 1:2592){
	i<- k%%72 +1;  j<-1+ (k-i)/72 
	x<-instr[i,j,];
	count<-sum(!is.na(x[index]))
	if ( count>50) {fm<-lm(instr[i,j,][index] ~ year[index])
		stat[i,j,4]<-coef(fm)[2]
		}
		}

##PLOT
	range(stat[,,4],na.rm=TRUE) #[1] 0.1734019  0.2377499

	quantile(stat[,,4],c(0,.05,.25,.5,.75,.95,1),na.rm=TRUE) #[1] -0.1734019  0.2377499
 	#         0%           5%          25%          50%          75%          95%         100% 
	#-0.173401891 -0.019404259  0.003856024  0.016665475  0.030756937  0.057195758  0.237749851 

	layout(1);par(mar=c(3,3,1,1)) 
 	image.plot(seq(-177.5,177.5,5),seq(-87.5,87.5,5),stat[,,4],xlab="LONGITUDE", ylab="LATITUDE",col=tim.colors(136)[c(seq(2,94,2),96:136)])#,
       	world(xlim=c(-87.5,87.5),ylim=c(-180,180),add=TRUE)


		#saved as solar/ccsp.h70.gif