##==== HEADER ========================================================================
## Ttools version v.1.6
##
## Purpose:
## This is an R toolbox to process temperatur time series data from British Columbia
## Ministry of Forest, Land and Natural Resources Operation
##
## Modification History:
## beta version v.1.0 code generation and adaption -- January/February 2013, hasi
## v.1.1 integration and mod. of preexisting tools from PERMOS & PermaSense , jn, hasi
## v.1.2 simplification in ts.df.runmean -- 7 March 2013, hasi
## v.1.3 tools 'elements', 'ts.normal' and 'ts.diffnorm' -- March 12 2013, hasi
## v.1.4 version numbers, 'count.sim', corrected 'ts.df.runmean'  -- March 19 2013, hasi
## v.1.4.1 correction in ts.normal  -- March 25 2013, hasi
## v.1.4.2 introduce ts.vars.mean and nice.abline -- April 3 2013, hasi
## v.1.5 revision of ts.agg, introduced ts.amp -- April 15 2013, hasi
## v.1.5.1 introduced ts.wdt.runmean, added quantiles to ts.stat -- April 20 2013, hasi
## v.1.5.2 correction in ts.reg for R-bug with posix-vectors -- May 11 2013, hasi
## v.1.5.3 update ts.stat to run with R 3.0.0
## v.1.5.4 addaption of ts.reg and ts.agg
## v.1.6	conversion to R time series class and first integration of related functions
## 
## Autors: 
## hasi (Andreas Hasler), andreas_hasler@gmx.ch, BC Forest Dept. / hiahaz
## jn (Jeannette Noetzli), jeannette.noetzli@geo.uzh.ch, PERMOS / University of Zurich
## 
## Description:
## makes diverse tools available of T timeseries data treatment (see list after header)
##
## Notations:
## 't' refers to time and 'T' to temperature
## "%Y-%m-%d" is the R notation for date formate "yyyy-mm-dd"
## 'd' or 'df' refers to data frame 
## 'ts' refers to time series (usually in a dataframe: e.g. ts.cut = cuts a data frame 
## 		which contains t in first columne) 
## 		note that the standart data formate used in this toolbox are "data frame objects" and 
##		not R "time series objects", however there are tools to transform between these objects
## 'Rts' refers to ts class objects, 
##
## Note for Rts: the use of ts objects is only recomended for long regular
## time series (e.g. monthly data), otherwise the acf function on  is 
##								
## Note!	t is always in POSIXt! This uses more memory than "Date" but allows common 
##		plotting and use of the same tools for all time series data.
##
## Compatibility:
## recomended to use Ttools with R version: R 3.0.0 
## veryfied function on R 3.0.0 on mac OS 10.7.5 and R 2.8.0 on PC Windows Vista
## Dependencies:
## library "base" and "stat" (usually loaded with installation)
## package "caTools" required (for runmean and all funcions using runmean)
##
##======================================================================================

## List of tools (functions):
## lower listed functions may depend on functions above. No separate copying recomended!
##--------------------------- # 

## &	concatenate strings with ampersand

## count(x,na.rm=TRUE); counts elements not NA (or all elements) in vector or df

## elements(x); gives all different elements of vector (no duplicates)

## count.sim(x); counts similar elements of vector and returns vector with count for each 

## df.maketime(d,format="auto"); make first columne posix from UTC in format specified
##								use format = "auto" if time given in "yyyy-mm-dd" or "yyyy-mm-dd HH:MM:SS"

## chk.class.t(t) / chk.class.d(d); checks class of vector(t) or df (d) and return t

## df.sort(d, col) sort data frame by columne

## quant01, quant09 (x,na.rm=TRUE); calculate 10% and 90% quantiles 

## posix(x) ;converts x = "yyyy-mm-dd" to POSIXt(GMT) if not allready POSIXt
## t2c(x) ;converts POSIXt to character "yyyy-mm-dd"; reverse of posix()

## pldev(x,y,pld="scr",plfn=plfn);x,y(device units); 
## 	pld = plot device: "scr" for screen, "png" for png file or "pdf"; 
##	plfn = filename(without extension))

## ts.cut(d,beg,end,format="%Y-%m-%d")	;d is input data frame with t in first col

## ts.reg(d,tbeg=NA,tend=NA,int=NA) ;makes regular timeseries: fills with NA
##	if tbeg/tend = NA start and end taken from d, otherwise timerange as defined
##	tbeg, tend in "%Y-%m-%d" or posix (exact)
##  if int is NA the timestep is taken from first two rows, otherwise given in "int" (minutes)

## ts.agg(indf,byform="%Y",fun="mean", qnum=1)	
## byform = "%Y-%m-%d %H" for hourly (daily, etc.) means, 
##	'fun' is aggreg. function; 'qnum' is number of values needed for valid output
##	returns POSIX if possible, otherwise numeric or character (for %m-%d)

## seas.mean <- function(d,cm=1); d=dataframe with POSIX in col 1, 
##	centre month: cm=1 seasons are DJF,MAM,JJA,SON; cm=2: JFM, AMJ,..; cm=3: FMA,... 

## ts.wdt.runmean(d,wdt,end=FALSE, qlim = 1) 
## own function to make weighted temporal running mean of data frames, SLOW!!!
## d = dataframe with t(POSIX) in first row, win = average window in days, 
## end (mean at end or centre), qlim for quality limit

## runmean(x(vector), k(window in rows), alg="C",endrule="mean",align)
##   	;function from library "caTools"

## ts.df.runmean(d,win(days), qlim = 1 (quality level; 1=all data needed), qout=FALSE, end=FALSE) 
##	;makes running means on data frames, uses "caTools"
## 	;returns list of "d" and "q" if qout = TRUE, otherwise dataframe "d"
##	;end= TRUE returns means at end of averaging period

## int.sh.gap(d,gapsize=1): linear interpolation of short gaps (gapsize in datapoints)

## day.mean(d); calculates daily means with gap filling dependent on measurement interval:
##			gapsize = 4 if hourly meas., = 2 if 2-hour intv., = 1 if less

## ts.stat(d,agglevel="%Y"); calulates list with vectors or data frames for 
##					mean, min, max, stdv and count

## deg.day(dm,daycount=TRUE); dm = df with daily means, returns list with df,s for pos 
##					and neg degree days: "thawdeg" and "freezdeg" 
##					if daycount=TRUE "thawdaycount", "freezdaycount" and daycount 
## 				output: list with dataframes: tdd(thawdeg), fdd(freezdeg),  
## 					tdc(thawdaycount),fdc(freezdaycount), dc(daycount)

## ts.MAT.plot(d,plv="T",tbeg=NA,tend=NA,title="");dataframe and variable name
##	timerange (tbeg, tend) in %Y (yyyy)(entire years considered)or posix (exact cuting)	

## ts.normal(d,ybeg=1961,yend=1990,yout=NA,spread=NA); 
##					function to calc daily normals over period ybeg:yend
## 					d must contain daily values, d[,1] should contain POSIX
## 					yout is output year (yyyy) or range(yyyy:yyyy), 
##					for NA output is only in %m-%d including "02-29"
##					spread="quant" gives 0.1 and 0.9 quantiles, 
##					spread="maxmin" gives max and min	
##					spread="stdv" gives norm +- sd (only for normallly dist. residuals)	

## ts.diffnorm(d,tbeg=1961,tend=1990,quant=FALSE)
## 					function to calc difference from daily normals over period tbeg:tend
## 					d must contain daily values, d[,1] should contain POSIX
## 					tbeg, tend are start and end year o period

## ts.aa(d,type="annual",simple=TRUE)
## 					function to calc AA (annual amplitude) from timeseries
## 					calc mean amplitude(max - min) for all years 
## 					type =  "daily" or "annual"
##					simple=TRUE returns only max-min value and uses Jan and Jul average for calc
##					simple=FALSE returns upper and lower deviation from mean and is based on
##							running mean stat. is about 5% higher values

## ts.vars.mean(d)
## 					function to calc a vector with means for each point in time  
##					averaging all variables in a time series data frame d
##					returns a vector as long as rows in d

## nice.abline(ix = "y",iy = 1); function to make ablines, ix is only on "y"=1 year
##					april 2013, hasi


     
						
##======================================================================================
## start functions:
##======================================================================================
##   &
## construct method to concatenate strings with ampersand; Henrik Bengtsson, 2005
# ------------------------------------------------------- # 

"&" <- function(...) UseMethod("&") 
"&.default" <- .Primitive("&") 
"&.character" <- function(...) paste(...,sep="")

##======================================================================================
##   count(x) counts all elements not NA in vector or df
# ------------------------------------------------------- # 
count <- function(x,na.rm=TRUE){
 	if (na.rm) x <- sum(!is.na(x)) else {
 		if (is.vector(x)) x <- length(x) else x <- (dim(x)[1])*(dim(x)[2])
 	}	
	return(x)
}

##======================================================================================
##   gives different elements of vector (no duplicates)
# ------------------------------------------------------- # 
elements <- function(x) x <- levels(as.factor(x))

##======================================================================================
##   counts similar elements of vector 
# ------------------------------------------------------- # 
count.sim <- function(x) {
y <- elements(x)
for (i in 1:length(y)) y[i] <- sum(y[i]==x,na.rm=TRUE)
return(as.numeric(y))
}

##======================================================================================
##   sort dataframe by columne
# ------------------------------------------------------- # 
df.sort <- function(d, col) d[order(d[,col]),]

##======================================================================================
## functions to make 10% and 90% quantile
# ------------------------------------------------------- # 
quant01 <- function(x,na.rm=TRUE) quantile(x,probs=0.1,na.rm=na.rm,names=FALSE) 
quant09 <- function(x,na.rm=TRUE) quantile(x,probs=0.9,na.rm=na.rm,names=FALSE)

##======================================================================================
##   make first columne posix from UTC
# ------------------------------------------------------- # 
df.maketime <- function(d,format="auto"){
t <- d[,1] 
t1 <- t[1] 	
## make posix 
if (format=="auto"){
	if (nchar(t1)==10) t <- posix(t)
	if (nchar(t1)==19) t <- strptime(t, format="%Y-%m-%d %H:%M:%S", tz="GMT")
} else t <- strptime(t, format=format, tz="GMT")
d[,1] <- as.POSIXct(t)
return(d)
}

##======================================================================================
##   chk.class.t(t) ; chk.class.d()
## checks t or first col for posix
## makes system sensitive: unix requires other plot window than windows
## return: chk.class.t(t) -> TRUE/FALSE, chk.class.d(d) -> t
# ------------------------------------------------------- # 
chk.class.t <- function(t)cl <- max(class(t)=="POSIXt")

chk.class.d <- function(d){
## check d for time in first row
t <- d[,1]
if (!max(class(t)=="POSIXt")) stop("no time format (POSIXt) in first columne!")	
return (t)	
}

##======================================================================================
##   posix() ; t2c()
## time conversion between posix and "yyyy-mm-dd"
# ------------------------------------------------------- # 

posix <- function(...)if(class(...)[1]!="POSIXt")strptime(...,format="%Y-%m-%d",tz="GMT")else(...)
t2c <- function(...) strftime(as.POSIXlt(...),format="%Y-%m-%d",tz="GMT")

##======================================================================================
##   pldev() 
## define grafic device type size 
## makes system sensitive: unix requires other plot window than windows
# ------------------------------------------------------- # 
pldev <- function(x,y,pld="scr",plfn="") {

syst <- .Platform$OS.type
fn <- ""& plfn &"."&pld
if (pld == "scr") if (syst == "unix")quartz (,x,y,pointsize=10) else windows(x,y,pointsize=10)
if (pld == "png") png(filename=fn,width=160*x,height=160*y,pointsize=23)
if (pld == "pdf") pdf(file=fn,width=x,height=y,pointsize=10)
if (pld == "stu") if (getOption("device")!="RStudioGD") dev.off()

} ## end pldev()
##======================================================================================
## Cut temperature series to a given interval  
## Jan 2013/hasi
## beg/end in format "%Y-%m-%d" or as POSIX 
## d is the input data frame with first col POSIXct
# ------------------------------------------------------- # 
ts.cut <- function(d,tbeg,tend,format="%Y-%m-%d"){

   	# convert to posix
	if (!chk.class.t(tbeg)) tbeg <-strptime(tbeg,format=format,tz="GMT") 
	if (!chk.class.t(tend)) tend <-strptime(tend,format=format,tz="GMT")
	
	tbeg <- as.POSIXct(tbeg)
	tend <- as.POSIXct(tend)

	#cut dates and temperature series to time range:
	lines   <- ((d[,1])>=tbeg) & ((d[,1])<=tend) & (!is.na(d[,1]))
	dCut    <- d[lines,]
	return(dCut)

}# end ts.cut()
##======================================================================================	
## function to fills data frame d with NAs where no data with regular 
## time steps for whole time range (first to last existing data in ts)
## with t-interval given by first two rows or with "int" in minutes
## date-time must be in first col and tbeg/tend is in posix or "yyyy-mm-dd"
## 2013 hasi
# ------------------------------------------------------- #
ts.reg <- function(d,tbeg=NA,tend=NA, int=NA) {

t <- chk.class.d (d)

## get time range
tr <- range(t, na.rm = TRUE)

## replace start and end
if (is.na(tbeg)) tbeg <- tr[1]
if (is.na(tend)) tend <- tr[2]

## if tbeg/end set cuts dataframe
if (!is.na(tbeg)| !is.na(tend)) d <- ts.cut(d, tbeg, tend)
	
## get interval
if (is.na(int)) intv <- as.numeric(difftime(t[2],t[1],unit="mins")) else intv <- int

#create regular time series ([intv] = in minutes)
rtime <- seq(from=tbeg,to=tend,by=intv*60,tz="GMT")
#create new data frame
reg<-data.frame(t=rtime)
d<-merge(reg,d,by="t",all.x=TRUE)

return(d)

} # end ts.reg

##======================================================================================
# Aggregate time series (daily, monthly, yearly ...)  
# indf is the input data frame with first col POSIXct 
# byform = "%Y-%m-%d %H" for hourly (daily, etc.) means
# aggregation over same month or day of all years is possible (used in "ts.normal")
# qnum is number of values needed
# Jan 2011, PERMOS/jn - modeified April, June 2013, hasi
# note!!!: timestamp for daily data is at beginning of day but data for the following 23:59 hours
# ------------------------------------------------------- #
ts.agg <- function(d,byform="%Y",fun="mean",qnum=1) {

	t <- chk.class.d (d)

	#make time format for aggregating and sorting
	dims    <- dim(d)
	t <- as.POSIXlt(t)
	aggChar <- strftime(t,byform)
	#print(aggChar[1])

	#variable names for later output
	vn <- names(d)

	#aggregate
	aggdf   <- aggregate(d[,2:dims[2]],by=list(aggChar),FUN=fun,na.rm=TRUE) 
	## note that time is added as character from  
	tmps    <- round(aggdf[,2:dims[2]],digits=3)

	#set aggregated values calculated with less than qnum values to NA  
	nval    <- aggregate(d[,2:dims[2]],by=list(aggChar),FUN=count)    
	notok   <- (nval[,2:dims[2]]<qnum)
	tmps[notok] <- NA

   	# make nice time labels for output
   	if(byform =="%Y-%m-%d %H") aggDate <- strptime(aggdf[,1],format=byform,tz="UTC")+3600
   	## add 3600 seconds for next hour (timestamp after averaging periode) 
   	if(byform =="%Y-%m-%d") aggDate <- posix(aggdf[,1]) 
   	if(byform =="%Y-%m") aggDate <- posix(paste(aggdf[,1],"-15",sep=""))
   	if(byform =="%Y") aggDate <- as.numeric(aggdf[,1])
	if(byform =="%m-%d") aggDate <- (aggdf[,1])
	if(byform =="%m") aggDate <- as.numeric(aggdf[,1])
	if(byform =="%d") aggDate <- as.numeric(aggdf[,1])
	if(byform =="%Y-%d") stop("aggragation non-sense!")

	out <- data.frame(t=aggDate,tmps)
  	names(out)<-vn

  	return(out)
}#end ts.agg


##======================================================================================	
## calculates seasonal means from dataframe with POSIX in first col
## means x month of quartal: (DJF,MAM,JJA,SON): cm=1, (JFM,...):cm=2, (FMA,...):cm=3
## returns posix with centered data
## uses ts.agg()
## 2013, hasi
# ------------------------------------------------------- # 
seas.mean <- function(d,cm=1){

## check for time
t <- chk.class.d (d)
##get interval of dataset	
int <- as.numeric(difftime(t[2],t[1],unit="days"),origin=0)

## make regular ts from beg to end and get length in days
d <- ts.reg(d)

## daily mean
dm <- ts.agg(d,byform="%Y-%m-%d",qnum=(1/int))

## montly mean
mm <- ts.agg(dm,byform="%Y-%m",qnum=25)

## 3-mothly mean
l <- dim(mm)[1]
c <- dim(mm)[2]
sm <- mm[2:(l-1),]
for (i in 2:(l-1)) {
	for (k in 2:c) sm[i-1,k] <- mean(c(mm[i-1,k],mm[i,k],mm[i+1,k]))}

## select 
t <- as.POSIXlt(sm[,1])
mt <- as.numeric(strftime(t,format="%m"))
qm <- mt - (trunc((mt-1)/3)*3)
lin <- (qm==cm)
sm <- sm[lin,]

return(sm)

} ## end of seas.mean

##======================================================================================	
## own function to make weighted temporal running mean of data frames
## note: this function is very slow due to the for loops included.
## a faster version for weighted running means may be developed based on runmean
## d = dataframe with t(POSIX) in first row, win = average window in days, 
## end (mean at end or centre), qlim for quality limit
## 2013, hasi
# ------------------------------------------------------- # 
ts.wdt.runmean <- function(d,wdt,end=FALSE, qlim = 1) {

## check for time
t <- chk.class.d (d)

##get interval of dataset	
int <- as.numeric(difftime(t[2],t[1],unit="days"),origin=0)

## avg window in number of rows
binwidth <- length(wdt)

## short time
tsh <- t[((binwidth/2)):(length(t)-binwidth/2)]
l<-length(tsh)
rmean <- numeric(l)
valid <- numeric(l)

## loop through variabbles
cn <- dim(d)[2]
for (k in 2:cn){
var <- as.numeric(d[,k])
isnna <- as.numeric(!is.na(var))
if (length(var) < binwidth) {
        stop("'var' too short for this 'win'.")
    }

##loop with running mean
for (i in 1:l){
	
	rmean[i] <- weighted.mean(var[i:(i+binwidth-1)],wdt,na.rm=TRUE)
	valid[i] <- weighted.mean(isnna[i:(i+binwidth-1)],wdt,na.rm=TRUE)
	
}

rmean[valid < qlim] <- NA
temp<- data.frame(time=tsh,var=rmean)
reg<-data.frame(time=t)
temp<-merge(reg,temp,by="time",all.x=TRUE)
if(end){
	var<-temp$var
	tsh<-temp$time
	tsh <- tsh[(binwidth/2+1):(length(tsh))]
	var <- var[1:(length(var)-(binwidth/2))]
	temp<- data.frame(time=tsh,var=var)
	temp<-merge(reg,temp,by="time",all.x=TRUE)
	}

d[,k] <- temp$var
}
	
return(d)

} # end ts.wdt.runmean()

##======================================================================================
library("caTools")
## for function "runmean()"
##Usage
##  runmean(x, k, alg=c("C", "R", "fast", "exact"),
##         endrule=c("mean", "NA", "trim", "keep", "constant", "func"),
##         align = c("center", "left", "right"))
## with x vector or matrix
## with k windows of n elements

##======================================================================================	
## function to make temporal running mean (package caTools required)
## d = dataframe with t(POSIX) in first row, win = average window in days, 
## end (mean at end or centre(or centre-1/2), mask for masking values if NA in win
## returns list of d and q if qout = TRUE
## Jan / Feb 2013, hasi
## March 2013, correct b2 to integer, introduce option end, hasi
# ------------------------------------------------------- # 
ts.df.runmean <- function(d,win, qlim = 1, qout=FALSE, end=FALSE) {

## check for time
t <- chk.class.d (d)
    
##get interval of dataset	
int <- as.numeric(difftime(t[2],t[1],unit="days"),origin=0)

## number of rows for averaging
binwidth <- win/int # x rows mean
b2 <- floor(binwidth / 2)

## make matix for values
cols <- dim(d)[2]
rows <- dim(d)[1]

## make df for quality parameter:
isnna <- d
isnna[,2:cols] <- as.numeric(!is.na(d[,2:cols]))
q <- isnna 						##preserve isnna for debuging

if (rows < (binwidth*qlim)) {
        stop(" timseries in 'd' too short for this averaging window 'win'.")
    }

## loop through cols
for (c in 2:(cols)){
	## calc runmean for each col
	d[,c] <- runmean(d[,c],k=binwidth, alg="C", endrule="mean")
	q[,c] <- runmean(q[,c],k=binwidth, alg="C", endrule="mean")
}

## define ends of quality indicator:
q[1:(1+ b2),2:cols] <- q[1:(1+ b2),2:cols]*((0:b2)/(2*b2)+0.5)
q[(rows - b2):rows,2:cols] <- q[(rows - b2):rows,2:cols]*((b2:0)/(2*b2)+0.5)

## round 
d[,2:cols] <- round(d[,2:cols],3)
q[,2:cols] <- round(q[,2:cols],2)

## filter for low quality
mkna <- matrix(FALSE,rows,cols)
mkna[,2:cols] <- (q[,2:cols] < qlim)
d[mkna] <- NA

if(end){
	d[(1+b2):rows,2:cols]  <- d[1:(rows-b2),2:cols]
	q[(1+b2):rows,2:cols]  <- q[1:(rows-b2),2:cols]
	d[1:b2,2:cols]  <- NA
	q[1:b2,2:cols]  <- NA
	}

## return with or without q
if (qout) return(list(d=d,q=q)) else return (d)

} # end ts.df.runmean()

##======================================================================================	
## linear temporal interpolation short value gaps 
## input: 	d[,1]=t in POSIX; d[,2:...]=data, gap size to fill: "gapsize"
## output: d
## implemented with runmean for performance reasons:
## watch out for boundry effects if gapsize large than data between gaps!!!
## Feb 2013, hasi
# ------------------------------------------------------- #
int.sh.gap <- function(d,gapsize=1){

## check for time
t <- chk.class.d (d)

k <- (gapsize*2)+1

vn <- dim(d)[2]
for (i in 2:vn){
	## make vectors without NAs
	var <- d[,i]
	lin <- is.na(var)
	a<-t[!lin]
	var <- var[!lin]
	#interpolate and add to data frame
	T <- approx(a,var,t,rule=1,method="linear")
	var <- T$y	
	lin2 <-  (runmean(lin,k=k, alg="C", endrule="mean")) > 0.5
	var[lin2] <- NA
	d[,i] <- var 
}

return(d)

} ## end int.sh.gap()

##======================================================================================	
## function to make daily means
## input: 	d[,1]=t in POSIX; d[,2:...]=data, 
##		gaps are filled up to 'gapsize'; gapsize depends in meas. rate
##			 gapsize = 2 if hourly m., 1 if 2-hour intv., 0 if less
## output: dm with "date format" 
## Feb 2013, hasi
# ------------------------------------------------------- # 
day.mean <- function(d){

## check for time
t <- chk.class.d (d)

##get interval of dataset	
int <- as.numeric(difftime(t[2],t[1],unit="days"),origin=0)

## fill gaps dependent on interval (gapsize = meas./day / 12)
gapsize <- floor((1/int)/12)
x <- int.sh.gap(d,gapsize)

## calculate aggregates:
##daily:
dm <- ts.agg(x,byform="%Y-%m-%d",qnum=(1/int))

return(dm)

} # end dm()

##======================================================================================	
## function to make time series statistics
## input: 	d[,1]=t in POSIX; d[,2:...]=data, 
##		of any aggreg level: e.g. "ym", "mm", "dm", "raw", "%m"(for moths of all years) 
##		aggregation level of statistics: agglevel= "%Y-%m-%d","%m-%d","%d",%m
## output: list with vectors or data frames for mean, min, max, stdv, count
## Feb 2013, hasi
# ------------------------------------------------------- # 
ts.stat <- function(d,agglevel="%Y"){

##get length
l<-dim(d)[2]

xmean <- ts.agg(d,byform=agglevel,fun="mean")
lin <- (dim(xmean)[1])+1
for (i in 2:l) xmean[lin,i] <- round(mean(d[,i], na.rm=TRUE),2)
xmean[lin,1] <- "tot"

xmin <- ts.agg(d,byform=agglevel,fun="min")
for (i in 2:l) xmin[lin,i] <- min(xmin[,i], na.rm=TRUE)
xmin[lin,1] <- "tot"

xmax <- ts.agg(d,byform=agglevel,fun="max")
for (i in 2:l) xmax[lin,i] <- max(xmax[,i], na.rm=TRUE)
xmax[lin,1] <- "tot"

q10 <- ts.agg(d,byform=agglevel,fun="quant01")
for (i in 2:l) q10[lin,i] <- round(quant01(d[,i], na.rm=TRUE),2)
q10[lin,1] <- "tot"

q90 <- ts.agg(d,byform=agglevel,fun="quant09")
for (i in 2:l) q90[lin,i] <- round(quant09(d[,i], na.rm=TRUE),2)
q90[lin,1] <- "tot"

xstdv <- ts.agg(d,byform=agglevel,fun="sd")
for (i in 2:l) xstdv[lin,i] <- round(sd(d[,i], na.rm=TRUE),2)
xstdv[lin,1] <- "tot"

xcount <- ts.agg(d,byform=agglevel,fun="count")
for (i in 2:l) xcount[lin,i] <- sum(xcount[,i], na.rm=TRUE)
xcount[lin,1] <- "tot"


stat <- list(mean=xmean,min=xmin,max=xmax,q10=q10,q90=q90,stdv=xstdv,count=xcount)

return(stat)

} # end ts.stat()

##======================================================================================	
## function to make time freezing thawind days
## input: 	dm: data frame with timeseries of daily mean (func. dm()), form=%Y-%m-%d 
## output: list with dataframes tdd(thawdeg), fdd(freezdeg), tdc(thawdaycount), 
## fdc(freezdaycount), dc(daycount)
## Feb 2013, hasi
# ------------------------------------------------------- # 
deg.day <- function(dm,daycount=TRUE){

## check for time
t <- chk.class.d (dm)

## reduce t to years
t <- as.POSIXlt(t)
aggChar <- strftime(t,"%Y")

#variable names for later output
vname <- names(dm)

## variable are of df
vn <- dim(dm)[2]
## make pos and neg dfs 
pos <- dm[,2:vn] * (dm[,2:vn] >= 0)
neg <- dm[,2:vn] * (dm[,2:vn] < 0)
## calc yearly sum
tdd <- aggregate(pos,by=list(aggChar),FUN=sum,na.rm=TRUE) 
names(tdd)<-vname
fdd <- aggregate(neg,by=list(aggChar),FUN=sum,na.rm=TRUE)
names(fdd)<-vname
degday <- list(tdd=tdd,fdd=fdd)

if (daycount){
#set aggregated values calculated with less than qnum values to NA  
	tdc    <- aggregate((pos>0),by=list(aggChar),FUN=sum,na.rm=TRUE) 
	names(tdc)<-vname
	fdc    <- aggregate((neg<0),by=list(aggChar),FUN=sum,na.rm=TRUE) 
	names(fdc)<-vname
	dc    <- aggregate(dm[,2:vn],by=list(aggChar),FUN=count)
	names(dc)<-vname
	degday <- list(tdd=tdd,fdd=fdd,tdc=tdc,fdc=fdc,dc=dc) 
}

return(degday)
}

##======================================================================================	
## function to plot annual and running mean 
## d = dataframe with t(POSIX) , timesteps 1 day or later 
## end (mean at end or centre), mask for masking values if NA in win
## Jan 2013, hasi
# ------------------------------------------------------- # 
ts.MAT.plot <- function(d,plv="T",tbeg=NA,tend=NA,title="",yrml=TRUE){

## check for time
t <- chk.class.d (d)

##get interval of dataset	
int <- as.numeric(difftime(t[2],t[1],unit="days"),origin=0)

## make internal dataframe with t and variable of interest
d <- data.frame(t=t,T=d[,plv])

## prepare to cut dataframe
if (!is.na(tbeg)){
	if(!(chk.class.t(tbeg))) tbeg <- posix(""&tbeg&"-01-01")
} else tbeg <- min(t,na.rm=TRUE)
if (!is.na(tend)){
 	if(!chk.class.t(tend))tend <- strptime(""&tend&"-12-31 23",format="%Y-%m-%d %H",tz="GMT")
} else tend <- max(t,na.rm=TRUE)
## make regular ts from tbeg to tend and get length in days
d <- ts.reg(d,tbeg,tend)
tr <- as.numeric(difftime(as.POSIXct(tend),as.POSIXct(tbeg),unit="days"),origin=0)

## calculate aggregates:
##daily:
if (int==1) dm <- d else dm <- day.mean(d)

##annuals
ym <- ts.agg(dm,byform="%Y",qnum=(345)) 

##running means
## monthly
mrm <- ts.df.runmean(dm,win=30,qlim=0.95)
## yearly
if (count(dm$t)>365) yrm <- ts.df.runmean(dm,win=365,qlim=0.95)

## set margins c(bottom, left, top, right)
par(mar = c(2, 3, 2, 2) + 0.2, mgp=c(2,0.7,0))   # Leave space for z axis

## general titles and labels 
title <- title
xlab <- "time"
ylab <- "T (°C)" ## (C)

## parameters for left plot (barplot)
l <- length(ym$T)
xlim <- c(0,l)+(1/6)
ylim1 <- c(-4.1,4.1)
if (max(!is.na(ym$T))==0) yaxt="n" else yaxt="s"
col <- ym$T >= 0
col[col] <- "orange"
col[ym$T < 0] <- "blue"

## parameters for right plot (lineplot)
ylim2 <- c(-38,38)
at<-posix(""&c(ym$t,ym$t[l]+1)&"-01-01")
barplot(ym$T,width=2/3,main=title,ylim=ylim1,space=0.5,col=col,xlim=xlim,yaxt="n")
##make ablines
abline(h=c(-4,-3,-2,-1,0,1,2,3,4),col="black",lty=2)
## make axis on right side 
axis (side=4,yaxt=yaxt)

## plot running mean


## right plot: line plots (T time series and running means of it)

par (new=TRUE)
plot(d$t,d$T, ylab=ylab,type="l",col="lightgrey",main=title,ylim=ylim2)
abline(v=as.POSIXct(at),col="black",lty=2)
lines(dm$t,dm$T,type="l",col="darkgrey")
lines(mrm$t,mrm$T,col="blue")
if (tr >= 365) lines(yrm$t,yrm$T,col="darkblue",lty=1,lwd=1)
if ((tr >= 365)&(max(!is.na(ym$T))>0)&yrml) lines(yrm$t,yrm$T*10,col="orange",lty=1,lwd=2)
if ((tr >= 365)&(max(!is.na(ym$T))>0)&yrml) lines(yrm$t,yrm$T*10,col="blue",lty=2,lwd=2)

abline(h=0)

} # end of "ts.MAT.plot()"

##======================================================================================
##======================================================================================	
## function to calc daily normals over period tbeg:tend
## d must contain daily values, d[,1] should contain POSIX
## tbeg, tend are start and end year o period
## march 2013, hasi
# ------------------------------------------------------- # 
ts.normal <- function(d,ybeg=1961,yend=1990,yout=NA,spread=NA){

## check for time
t <- chk.class.d(d)
tn <- names(d)[1]

## cut data frame
tbeg <- posix(""&ybeg&"-01-01")
tend <- posix(""&yend&"-12-31")
d <- ts.cut(d,tbeg,tend)

## calcualte normals
dnorm <- ts.agg(d,"%m-%d")

## make month and day labels
yearday <- elements(strftime(as.POSIXlt(t),"%m-%d"))
dnorm[,1] <- yearday

if (!is.na(spread)) {
if (spread == "quant"){
	sp_low <- ts.agg(d,"%m-%d",fun=quant01)
	sp_low[,1] <- yearday
	sp_high <- ts.agg(d,"%m-%d",fun=quant09) 
	sp_high[,1] <- yearday
}

if (spread == "minmax"){
	sp_low <- ts.agg(d,"%m-%d",fun=min)
	sp_low[,1] <- yearday
	sp_high <- ts.agg(d,"%m-%d",fun=max) 
	sp_high[,1] <- yearday
}

if (spread == "stdv"){
	sp_low <- dnorm - ts.agg(d,"%m-%d",fun=sd)
	sp_low[,1] <- yearday
	sp_high <- dnorm + ts.agg(d,"%m-%d",fun=sd) 
	sp_high[,1] <- yearday
}
}

if (!is.na(yout[1])) {
	l <- length(yout)
	if (l==1) tb <- posix(""&yout&"-01-01") else tb <- posix(""&yout[1]&"-01-01")
	if (l==1) te <- posix(""&yout&"-12-31") else te <- posix(""&yout[l]&"-12-31")

	#create regular time series ([intv] = in minutes)
	rtime <- seq(from=tb,to=te,by=24*3600,tz="GMT")

	## make df for yout with normals
	x <- data.frame(t=rtime, x=strftime(as.POSIXlt(rtime),"%m-%d"))
	names(dnorm)[1] <- "x"
	dnorm <- merge(x,dnorm,by="x", all.x=TRUE) 
	dnorm <- dnorm[,(names(dnorm)!= "x")]
	## order data by date
	dnorm<-df.sort(dnorm,1)

	if (!is.na(spread)) {
		## make df for yout with spreads
		names(sp_low)[1] <- "x"
		sp_low <- merge(x,sp_low,by="x", all.x=TRUE) 
		sp_low <- sp_low[,(names(sp_low)!= "x")]
		sp_low <- df.sort(sp_low,1)

		## make df for entire period with normals
		names(sp_high)[1] <- "x"
		sp_high <- merge(x,sp_high,by="x", all.x=TRUE) 
		sp_high <- sp_high[,(names(sp_high)!= "x")]
		sp_high <- df.sort(sp_high,1)
	}
}

if (!is.na(spread)) dnorm <- list(norm=dnorm,sp_low=sp_low,sp_high=sp_high)

return(dnorm)

} # end of "ts.norm()"

##======================================================================================
##======================================================================================	
## function to calc difference from daily normals over period tbeg:tend
## d must contain daily values, d[,1] should contain POSIX
## tbeg, tend are start and end year o period
## march 2013, hasi
# ------------------------------------------------------- # 
ts.diffnorm <- function(d,ybeg=1961,yend=1990){

## check for time
t <- chk.class.d(d)

## calc normal
dnorm <- ts.normal(d,ybeg,yend)
dnorm$x <- dnorm$t
dnorm <- dnorm[,(names(dnorm)!= "t")]

## make df for entire period with normals
x <- data.frame(t=t, x=strftime(as.POSIXlt(t),"%m-%d"))
dn <- merge(x,dnorm,by="x", all.x=TRUE) 
dn <- dn[,(names(dn)!= "x")]

## order data by date
dn<-dn[order(dn$t),]

## subtract data frames
diff <- d - dn
diff$t <- t

return(diff)

} # end of "ts.diffnorm()"

##======================================================================================	
## function to calc amplitude from timeseries
## calc mean amplitude(max - min) for all years 
## type = (daily or annual amplitude): type = "daily" or "annual"
## simple = TRUE for annual amplitude = (mean (July) - mean (Jannuary))/2
## April 2013, hasi
# ------------------------------------------------------- # 
ts.amp <- function(d,type="annual",simple=TRUE){

## check for time
t <- chk.class.d (d)

## check intervall for daily amplitude	
int <- as.numeric(difftime(t[2],t[1],unit="days"),origin=0)
if ((type=="daily")&((1/int) < 4)) stop("less than four measurements per day do not allow option daily!")

## running mean or aggregation
if(simple){
	c <- length(d[1,])
	if (type == "annual") {
		n <- ts.agg(d,"%m","mean",qnum = 15)
		c <- length(n[1,])
		jan <- n[(n[,1]==1),]
		jul <- n[(n[,1]==7),]
		amp <- (jul[2:c] - jan[2:c])/2
	}	
	if (type == "daily") {
		n <- (ts.agg(d,"%Y-%m-%d",fun="max")-ts.agg(d,"%Y-%m-%d",fun="min"))/2
		amp <- n[1,2:c]
		amp[1:(c-1)]<-NA
		for (i in 2:c) amp[1,i-1] <- mean(n[,i],na.rm=TRUE) 
	}
}else {
	
	## data frame for output
	c <- length(d[1,])
	amp <- as.data.frame(matrix(NA,4,c-1))
	names(amp) <- names(d)[2:c]
	row.names(amp) <- c("atop","abot","a","q")

	## annual amplitude
	if (type == "annual") {
		## calc normals (two years because min around dec 31st)
		n <- ts.normal(d,1900,2100,2012:2013)
		## smoothing for montly means
		n <- ts.df.runmean(n,win=30,qlim=0.95)

		## calc amplitude: mean ( max-mean, mean-min)
		for (i in 2:c) {
			amp[3,i-1] <- mean(n[,i],na.rm=TRUE)
			amp[1,i-1] <- max(n[,i],na.rm=TRUE) - amp[3,i-1]
			amp[2,i-1] <- amp[3,i-1] - min(n[,i],na.rm=TRUE)
			amp[3,i-1] <- mean(c(amp[1,i-1],amp[2,i-1]),na.rm=TRUE)
			amp[4,i-1] <- count(d[,i],na.rm=TRUE)*int / (365 * 2)
		}
	}
	
	## daily amplitude
	if (type == "daily") {
		## substracting daily means from data
		dm <- day.mean(d)
		top <- ts.agg(d,"%Y-%m-%d","max")[,2:c] - dm[,2:c]
		bot <- dm[,2:c] - ts.agg(d,"%Y-%m-%d","min")[,2:c]
		## calc amplitude: mean ( max-mean, mean-min)
		for (i in 2:c) {
			amp[1,i-1] <- mean(top[,i-1],na.rm=TRUE)
			amp[2,i-1] <- mean(bot[,i-1],na.rm=TRUE)
			amp[3,i-1] <- mean(c(amp[1,i-1],amp[2,i-1]),na.rm=TRUE)
			amp[4,i-1] <- count(d[,i],na.rm=TRUE)*int / (365 * 2)
		}
	}	
}

return(amp)

} # end of "ts.amp()"

##======================================================================================	
## function to calc a TS with means of all variables in a time series for each point in time
## returns a vector as long as rows in d
## april 2013, hasi
# ------------------------------------------------------- # 
ts.vars.mean <- function(d){

## check for time
t <- chk.class.d(d)

## get rid of first col
x<- d[,2:dim(d)[2]]

## calc means
l <- dim(d)[1]
varmean <- vector(,l)
for (i in 1:l) varmean[i] <- mean(x[i,],na.rm =TRUE) 

return(varmean)

}
##======================================================================================	
## function to make ablines
## april 2013, hasi
# ------------------------------------------------------- # 
nice.abline <- function(ix = "y",iy = 1){

if (ix == "y") x <- as.POSIXct(posix(""&(1950:2050)&"-01-01"))
y <- ((0:20)-10)*iy
abline(h=y,v=x,lty=2)
abline(h=0,lty=1,lwd=1.1)

return(d)

}

##======================================================================================	
## function to convert dataframe class timeseries to R time series class
## makes time series object with 24 observations for int="hour"
## 
## june 2013, hasi
# ------------------------------------------------------- # 
df2Rts <- function(d, int="auto"){

t <- chk.class.d (d)
cols <- dim(d)[2]

## get time range
tr <- range(t, na.rm = TRUE)

## get start day
st <- t2c(tr[1])

## get interval
if (int=="auto") intv <- as.numeric(difftime(t[2],t[1],unit="days")) 
## make monthly if int="month"
if (int=="day") d <- day.mean(d)
if (int=="month") d <- ts.agg(d,"%Y-%m")
if (int=="season") d <- seas.mean(d,cm=2)

if (intv==1/24) int <- "hour"
if (intv==1) int <- "day"
if (intv >= 28) int <- "month"
if (intv > 90) int <- "season"

if (int=="hour"){
	st <- 1
	f <- 24
}
if (int=="day"){
	st <- 1
	f <- 1
}
if (int=="month"){
	st <- as.numeric(c(substr(st,0,4),substr(st,6,7)))
	f <- 12
}
if (int=="season"){
	st <- as.numeric(c(substr(st,0,4),substr(st,6,7)))
	st[2] <- ceiling(st[2] / 3) 
	f <- 4
}

Rts <- ts(d[,2:cols],start=st,frequency=f)

return(Rts)

}

##======================================================================================	
## function to convert a R time series class to a dataframe  
## june 2013, hasi
## not yet implemented for multivariat ts!!!!
# ------------------------------------------------------- # 
Rts2df <- function(Rts){
	
	x <- Rts
	d <- data.frame( t = c(time(x)),x = c(x))
	d$t <- d$t
	
	return(d)
}

##======================================================================================	
## plot decomposition of time series class object
## june 2013, hasi
# ------------------------------------------------------- # 
plot.ts.decomp <- function(Rts){

plot(decompose(Rts))

return(d)

}
