##SCRIPT TO CHECK ADJUSTMENTS FOR ALICE SPRINGS
#https://notalotofpeopleknowthat.wordpress.com/2012/03/15/an-adjustment-like-alice/
#https://cliscep.com/2017/02/06/instability-of-ghcn-adjustment-algorithm/


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

##UTILITY FUNCTIONS
require(ggplot2)
require(reshape2)

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

##download Matthews versions
 #returns list of unzipped files 
download_versions=function(dset="unadj") {
idmx=list(adj=c("17jan2012","11mar2012","9May2012","21jun2012","9oct2012","27mar2013","22jul2014"),
  unadj=c("22jul2014"))
idm=idmx[[dset]]
K=length(idm)
alias=c(adj="qca",unadj="qcu")

version=rep(list(NA),K)
names(version)=idm

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

#compression method 1
handle=gzfile("temp.gz")
work=readLines(handle,n=-1)
close(handle)
version[[j]]=work
rm(work)
}
return(version)
}
##method 2 
#system(paste("gunzip -f","temp.gz", sep=" "))
#work=readLines("temp")
#rm("temp")
#rm("temp.gz")

##
#function to collate by station id
collate_station=function(idn,dset="adj") {
  count=c(adj=7,unadj=1)
  K=count[dset]
  if(dset=="adj") version=versiona else version=versionu
  out=NULL
  for(i in 1:K) {
   station= get_station(idn,work=version[[i]])
   out=ts.union(out,station)
  }
  if(dset=="adj") dimnames(out)[[2]]= names(version)
  return(out)
} 	

##collect plotting dataframe
collect_plot=function() {
 B=array(NA,dim=c(nrow(adj),7))
 for(j in 1:7) B[,j]=adj[,j]-unadj
 dimnames(B)[[2]]=names(versiona)
 work=data.frame(year=c(time(adj)),B)
 names(work)[2:8]=names(versiona)
 x=melt(work,id.vars="year")
 names(x)[2]="version"
 return(x)
}
 
get_info=function(dset=versionu[[1]]) {
 work=dset
 y=substr(work,16,19)
 x=work[!y=="TAVG"]
 x=x[3:length(x)]
 writeLines(x,"d:/temp/temp1.txt")
 info=read.fwf("d:/temp/temp1.txt",widths=c(11,9,10,7,1,22),fill=T,
    colClasses=c("character",rep("numeric",3),rep("character",2)))
 info=info[,c(1:4,6)]
 names(info)=c("id","lat","long","alt","name")
 info$name=gsub(" +$","",info$name)
 return(info)
}
   
versionu=download_versions(dset="unadj")
  length(versionu)
versiona=download_versions(dset="adj")
  length(versiona)
info=get_info()


adj=collate_station(idn=id0,dset="adj")
unadj=collate_station(idn=id0,dset="unadj")

# year= floor(range(c(time(unadj))))
#this writes collated data
#out=data.frame(year=rep(year[1]:year[2],each=12),month=rep(1:12,136),unadj,adj)
#write.csv(out,file="d:/2017/alice_springs_versions.csv",row.names=F)

B=collect_plot()
p=ggplot(data=B,aes(year,value,group=version,col=version))+geom_line()+labs(x="",y="Adjustment (deg C)")+
   ggtitle ("GHCN Adjustments: Alice Springs")
p
ggsave(p,file="d:/climate/images/observations/alice-springs.png",w=5,h=5)

##WELLINGTON
info[grep("WELLINGTON",info$name),]
id0="50793436001"
adj=collate_station(idn=id0,dset="adj")
unadj=collate_station(idn=id0,dset="unadj")

B=collect_plot()
p=ggplot(data=B,aes(year,value,group=version,col=version))+geom_line()+labs(x="",y="Adjustment (deg C)")+
   ggtitle ("GHCN Adjustments: Wellington NZ")
p
ggsave(p,file="d:/climate/images/observations/ghcn_adjustment_wellingtonNZ.png",w=5,h=5)



##GRAND CANYON
info[grep("DETROIT LAKES",info$name),]
id0=info$id[grep("DETROIT LAKES",info$name)] # "42500212142"
adj=collate_station(idn=id0,dset="adj")
unadj=collate_station(idn=id0,dset="unadj")

B=collect_plot()
p=ggplot(data=B,aes(year,value,group=version,col=version))+geom_line()+labs(x="",y="Adjustment (deg C)")+
   ggtitle ("GHCN Adjustments: Detroit Lakes MN")
p

##ORLAND
info[grep("ORLAND",info$name),]
id0="42500046506"
adj=collate_station(idn=id0,dset="adj")
unadj=collate_station(idn=id0,dset="unadj")

B=collect_plot()
p=ggplot(data=B,aes(year,value,group=version,col=version))+geom_line()+labs(x="",y="Adjustment (deg C)")+
   ggtitle ("GHCN Adjustments: Orland")
p