####################
## FUNCTIONS TO COLLATE NEW USHCN INFORMATION
#
# get.info.ushcn2 - returns informaiton from ftp://ftp.ncdc.noaa.gov/pub/data/ushcn/v2/monthly/ushcn-stations.txt
#
# get.data.ushcn2 - returns monthly time series start 1895 of 1218 USHCN stations 
#
# extractf - utility to return time series from datra frame where not all years necessarily present

get.info.ushcn2=function() { 
	url="ftp://ftp.ncdc.noaa.gov/pub/data/ushcn/v2/monthly/ushcn-stations.txt"
           # [3] "013160  32.8347  -88.1342   38.0 AL GAINESVILLE/LOCK               011694 ------ ------ +6"
 	download.file(url,"temp.dat")
	info=read.fwf("temp.dat",widths=c(6,9,10,8,2,1,31,6,1,6,1,6,3),na.strings="------",fill=TRUE)
	info=info[,c(1:5,7,8,9,11,13)]
	names(info)=c("id","lat","long","alt","state","name","code1","code2","code3","code4")
	info$name=gsub(" +$","",info$name)
	return(info)
}

	#info=get.info.ushcn2() 

get.data.ushcn2=function() {
	url="ftp://ftp.ncdc.noaa.gov/pub/data/ushcn/v2/monthly/9641C_200904_F52.avg.gz"
	download.file(url,"temp.gz",mode="wb")
 	handle=gzfile("temp.gz")
	Data=readLines(handle,n=-1)
 	   #[1] "01108431895   487X   445X   599X   662X   726X   797X   833X   826X   822X   635X   561E   491E   657E"
	   # [2] "01108431896   484X   535X   568X   702E   778X   792X   820X   835X   786X   665X   599X   481X   671E"
	   # all values in column 7 are 3
 	N=length(Data);N #[1] 139973
	   #was #555080 in old USHCN #increase to 562586
	widths0=cumsum( c(6,1,4,rep(c(6,1),13)) );
	K=length(widths0)
	widths0=cbind(c(1,widths0[1:(K-1)]+1),widths0) 
	hcn=array(NA,dim=c(length(Data),K));
	hcn=data.frame(hcn)
	for(i in 1:K) hcn[,i]=substr(Data,widths0[i,1],widths0[i,2])
	for (i in c(2,3,seq(4,28,2)) ) hcn[,i]=as.numeric(hcn[,i])
	names(hcn)[c(1,3)]=c("id","year")
	temp= hcn== -9999; hcn[temp]=NA
	   #x=tapply(!is.na(hcn$id),hcn$id,sum)
	   #tapply(!is.na(x),x,sum)
		# 106  107  108  110  111  112  113  114  115 
		#   1    2    1    5    2    4    8    4 1191 
	id=levels(factor( hcn$id));
	M=length(id);M #1218

	X=NULL
	for( i in 1:M) {
	  temp=hcn$id==id[i]; sum (temp)
	  working=hcn[temp,]
	  name0=paste(id[i])
	  X=ts.union(X,extractf(working))   #this is in Fahrenheit tenths
	}
	dimnames(X)[[2]]=id
	X=ts(X,start=1895,freq=12)
	return(X)
}


##utility function for USHCN
extractf=function(A){
	  target=min(A$year):max(A$year);N=length(target)
	  dummy= array(NA,dim=c(N,12))
	  test=match(target,A$year);temp=!is.na(test)
	  dummy[temp,]=as.matrix(working[test[temp],seq(4,26,2)])
	  extractf= ts( c(t(dummy)),start=target[1],freq=12)/10
	  extractf
	}


	# hcn=get.data.ushcn2()


#utility function for matching id form
numbertochar=function(x) if (nchar(x)==5) paste("0",x,sep="") else paste(x)


Hcodes=read.table("http://www.climateaudit.org/data/station/ushcn/history_codes.dat",sep="\t",skip=1,header=TRUE,fill=TRUE)
dim(Hcodes) # 199 6
USHCN.history=readLines("http://cdiac.ornl.gov/ftp/ushcn_monthly/station_history")

get.history.ushcn=function(id0,ushcn.history=USHCN.history,hcodes=Hcodes,Details=details) {
	k=(1:nrow(Details))[Details$id==id0]
	title0=paste(Details$state[k],Details$name[k]);title0 #  "046506"
	histid=as.numeric(substr(ushcn.history,1,6))
	temp1=!is.na(hcodes[,1]);hcodes=hcodes[temp1,]
	widths0=hcodes[,1:2]
	what0=hcodes[,3];
	M=nrow(widths0)

	temp=!is.na(match(histid,id0));sum(temp)
	x=ushcn.history[temp][2:sum(temp)]

	X=array(NA,dim=c(sum(temp)-1,M) );
	X=data.frame(X)
	f=function(i){ f=substr(x,widths0[i,1],widths0[i,2]);f}
	for(i in 1:M){X[,i]=f(i)}
	names(X)=hcodes[,4]
		#X$name=gsub(" +$","",X$name)
		#data.frame(t(X))
	K=nrow(X)

  ##MOVE REPORT
	print(paste(id0,title0,sep="   "))
	A=data.frame(X[,c(2,21,6,23,65,31,32,36,37,51)]) #22 is dist flag
	A$elev_chg=X$elev;	A$hgt_chg=X$hgt_temp;
	A$dir_move[A$dir_move=="999"]=NA
	A$hgt_chg[A$hgt_chg=="99"]=NA

	temp=(substr(A$dist_prev_loc,1,1)=="8");A=A[!temp,]
	temp=(substr(A$dist_prev_loc,1,1)=="9");
	A$dist_prev_loc=as.numeric(A$dist_prev_loc)
	temp=(A$dist_prev_loc==999);A$dist_prev_loc[temp]=NA
	temp=(A$dist_prev_loc>=900)&!is.na(A$dist_prev_loc);A$dist_prev_loc[temp]=A$dist_prev_loc[temp]%%900
	#A$dist_prev_loc
	#v=c("dist_prev_loc");for (i in v){A[,i]=as.numeric(A[,i])}
	#v=c("dist_prev_loc");temp=(A[,v]==999);A[,v][temp]=NA
	v=c("dist_prev_loc");A[,v]=A[,v]/10

	#A$elev_chg=c(NA,diff(as.numeric(A$elev_chg)) )
	#A$hgt_chg[A$hgt_chg==99]=NA;A$hgt_chg=c(NA,diff(A$hgt_chg))

	obs=paste(X[,64],X[,65],sep="");
	temp=(substr(obs,1,1)=="9")&(substr(obs,4,4)=="9")& !(substr(obs,1,2)=="99")& !(substr(obs,3,4)=="99")
	A[temp,"obstime_temp"]=substr(obs[temp],2,3)
	temp=(A$obstime_temp=="99");A$obstime_temp[temp]=NA
	move_report=A[,c(1:4,11,12,5:10)];move_report

##also see http://mi3.ncdc.noaa.gov/mi3qry/login.cfm
	
  ##instrumentation report  #SHELTER THERMOMETER
	print(paste(usid,title0,sep="   "))
	instrument_report= X[,c(2:3,31,40,46, 32,36,37,47,48,52,53,51,54,55,61)]

	get.history.ushcn=list(id=id0,name=title0,move=move_report,instrument=instrument_report)
	get.history.ushcn
}
	
	#k=grep("ORLAND",details$name)[1]
	#id0=details$id[k];id0 #  "046506"
	#get.history.ushcn(id0)