##SCRIPT TO CHECK ADJUSTMENTS FOR ALICE SPRINGS

##GET ALICE SPRINGS ID
#load("d:/climate/data/observations/station/ghcn/ghcn.info.tab")
#info=ghcn.info
#info[grep("ALICE",info$name),]
#id0="50194326000"

#ADJUSTED VERSIONS
idm=c("17jan2012","11mar2012","9May2012","21jun2012","9oct2012","27mar2013","22jul2014")

##UTILITY FUNCTIONS

#make timeseries from data frame extracted from GHCN data
make_timeseries=function(A) {
  names(A)[1]="year"
  year=min(A$year):max(A$year)
  N=length(year)
  x=array(NA,dim=c(N,12))
  x[match(A$year,year),]=as.matrix(A[,2:13])
  y=ts(c(t(x)),start=min(A$year),freq=12)
  y[y<  -900]=NA
  y=y/100
  return(y)
}

#extract station data from bulk GHCN dset by id
get_station=function(idn, workx=work, out="d:/temp/temp1.txt"){
# workx=work; out="d:/temp/temp1.txt"
  temp=substr(workx,1,11)==idn;
  # sum(temp)
   work1=workx[temp]
   work1=work1[grep("TAVG",work1)]	
   writeLines(work1,out)
  fred=read.fwf(out,widths=c(11,4,4,rep(c(5,3),12)) )
  fredx= fred[,seq( 2,26,2)]
  y=make_timeseries(fredx)
   return(y)
}
######################

### GUNZIP NOT NEEDED 
# file.copy("d:/climate/scripts/compression/gzip/gunzip.exe", "d:/temp/gunzip.exe")
###


##COLLATE ADJ VERSIONS
K=length(idm)
station=rep(list(NA),K)
names(station)=idm



setwd("d:/temp")
rm("temp.gz")


for(i in 1:K) {
loc=paste("http://www.climateaudit.info/data/observations/gridcell/ghcn/ghcnm.tavg.latest.qca.tar.",idm[i],".gz",sep="")
download.file(loc,"d:/temp/temp.gz",mode="wb")
setwd("d:/temp")

#compression method 1
handle=gzfile("temp.gz")
work=readLines(handle,n=-1)
close(handle)
rm("temp.gz")

##method 2 
#system(paste("gunzip -f","temp.gz", sep=" "))
#work=readLines("temp")
#rm("temp")
#rm("temp.gz")

test=get_station(idn=id0)
station[[i]]=test
rm(test)
}


##UNADJUSTED VERSIONS
unadj=rep(list(NA),2)
idu=c("9oct2012", "22jul2014")

for(i in 1:2) {
loc=paste("http://www.climateaudit.info/data/observations/gridcell/ghcn/ghcnm.tavg.latest.qcu.tar.",idu[i],".gz",sep="")
download.file(loc,"d:/temp/temp.gz",mode="wb")

setwd("d:/temp")

#compression method 1
handle=gzfile("temp.gz")
work=readLines(handle,n=-1)
close(handle)
rm("temp.gz")

#compressoin method 2
#system(paste("gunzip -f","temp.gz", sep=" "))
#work=readLines("temp")

test=get_station(idn=id0)
unadj[[i]]=test
rm(test)
}


##COMPARE UNADJ
work=ts.union(unadj[[1]],unadj[[2]],NA)
work[,3]=work[,2]-work[,1]
plot.ts(work[,3])
  #no difference
unadj=unadj[[2]]


####COLLECT ADJUSTED AND ADJ DIFFEERENCES
A=unadj
for(i in 1:K) A=ts.union(A,station[[i]])
dimnames(A)[[2]]=c("unadj",idm)

B=array(NA,dim=c(nrow(A),7))
for(j in 1:7) B[,j]=A[,j+1]-A[,1]
dimnames(B)[[2]]=idm
work=data.frame(year=c(time(A)),B)
names(work)[2:8]=idm

##PLOT
library(ggplot2)
library(reshape2)
x=melt(work,id.vars="year")
names(x)[2]="version"
ggplot(data=x,aes(year,value,group=version,col=version))+geom_line()+labs(x="",y="Adjustment (deg C)")+
   ggtitle ("GHCN Adjustments: Alice Springs")