##==== HEADER ========================================================================
## 
## Purpose:
## This is a script to process temperatur data from British Columbia
## Ministry of Forest, Land and Natural Resources Operation
##
## Modification History:
## code generation -- January / Feburary 2013, Hasi
## introduced quantiles -- April 2013, hasi
## correct rounding -- May 2013, hasi
##
## Autors: 
## Hasi (Andreas Hasler), andreas_hasler@gmx.ch, BC Forest Dept.
## 
## Description:		Dataframes are loaded(files TTSws.RData; run "load_TTS.R" 
## 				if not available) and statistics is generated. 
##
## Dependencies:		Ttools.R (own toolbox for temperature timeseries treatment
## 				caTools  (library)
##
## I/O, variables and units:
## input:			list of GST and WX sensorIDs (given by csv file specified 
##				in 'gstlist' and 'wxlist' )
##				list of reference variables if corr = TRUE; see*
##				path (given below under control variable "path") 		
##				temperature time series (TTS) data from RData files
##				 
## variables:		sensorID = ID for selection on the files to be loaded
##				path = absolute path of data files
##				TTS-xxx = dataframe with all temperature data per site
##		
## output:			dataframes and lists of dataframes (and corr. csv files)  
##					per field site (XXX): 
##				DMT_XXX: daily mean temperatures
##				YRM_XXX$d/q: running mean annual temperatures and quality 
##					(q=1: all data, q=0.5: half of data, ...)	
##				MAT_XXX$mean/min/max/stdv = mean, min ect of running mean 
##					annual temperature 
##				DeD_XXX$tdd(thawdegday), fdd(freezedegday), tdc(thawdaycount), 
## 					fdc(freezedaycount), dc(daycount)
##
## units:			t is in inported as "mm-dd-yyyy HH:MM:SS AM/PM"
##				t is internaly in R POSIXct
##				T, T1, T2 is temperature in ?C
##
## ** note regarding correction of total values for h
##			
##
##						
##======================================================================================

## clear current workspace to avoid env. dependency
#rm(list=ls())

##======================================================================================
##   control variables
##======================================================================================

## set file name for meta data file
wxlist <- "meta/PermafrostWxID.csv"
gstlist <- "meta/PermafrostProbeID.csv"
reflist <- "meta/MAT_ref.csv"

## set correct true if 
correct <- TRUE

## set qmat as min precent of data available for RMAT calc
qmat <- 99

##======================================================================================
##   preparation
##======================================================================================


## load data to workspace
if (!exists("TTSWS")) load("RData/TTSws.RData")

## load toolbox "Ttools" (requires R-library "caTools")
source("R_source/Ttools.R")

## saving csv file with MATs
save <- FALSE

## get sensor lists and list length
wxlist <- read.csv(wxlist, as.is=TRUE,sep=";")
if(dim(wxlist)[2]==1) wxlist <- read.csv(wxlist, as.is=TRUE)
gstlist <- read.csv(gstlist, as.is=TRUE,sep=";")
if(dim(gstlist)[2]==1) gstlist <- read.csv(gstlist, as.is=TRUE)
lwx<-length(wxlist$site)
lgst<-length(gstlist$site)

## make one large list with all sensors:
slist <- rbind(wxlist,gstlist)
slist <- slist[order(slist$site),]
## get list of sites:
sites <- levels(factor(slist$site))
lsl <- length(sites) ## length of short list

## reflist for overall MAT value and spreads
if (correct){
	x <- read.csv(reflist, as.is=TRUE,sep=";")
	if(dim(x)[2]==1) reflist <- read.csv(reflist, as.is=TRUE) else reflist<-x
}

## dummy for workspace indication
STAWS <- TRUE

##======================================================================================
##======================================================================================
##======================================================================================
##   main loop through sites

for (k in 1:lsl){	
	## get site name	
	site <- sites[k]
 	print("processing "&site)
 	
	## geting dataframe of respective sensor assigned to "d"
	dfn <- paste("TTS_",site,sep="") 
	d <- get(dfn)
	l <- dim(d)[2]

	## reflist for site for MAT calcualtion
	if (correct) shreflist <- reflist[(reflist$site == site),]
		
##======================================================================================
##   calculating stats (DMT, RMAT, MAT, DeD, SOTO )
##======================================================================================
	
	## **********************************************************
	## make daily mean temperature (DMT):
  print("calc daily mean T")
	dm <- day.mean(d)
	## asign df as 
	dfn <- "DMT_" & site 
	assign(dfn,dm)
	if (save) write.csv (dm,"CSV_pro/"&dfn&".csv",row.names=FALSE)

	## **********************************************************
	## running annual means
	if (site == "PIN") q <- 95 else q <- qmat 
	print("calc RMAT where "&q&"% of data present")
	rmat <- ts.df.runmean(dm,365,q/100,qout=TRUE,end=TRUE)
	## asign df as list with 'd' for data and 'q' for quality
	dfn <- "RMAT_" & site 
	assign(dfn,rmat)
	if (save) write.csv (rmat,"CSV_pro/"&dfn&".csv",row.names=FALSE)

	## **********************************************************
	## make MAT statistics (mean, max, min, stdv, count)
	## call MAT statistic
  print("calc MAT statistics")
	stat <- ts.stat(rmat$d)
	
	if (correct) {
		
		## only if varaible to correct in shreflist
		l <- dim(shreflist)[1]
		if (l > 0){
		## get line number with totals	
		lin <- stat$mean$t == "tot"
		mmat <- as.numeric(stat$mean[lin,])
 		Tmax <- as.numeric(stat$max[lin,])
		Tmin <- as.numeric(stat$min[lin,])
 		q90 <- as.numeric(stat$q90[lin,])
		q10 <- as.numeric(stat$q10[lin,])		
		SD <- as.numeric(stat$stdv[lin,])
		nam <- names(stat$mean)

		## loops throug reference list

		for (i in 1:l){
	
			## get Tcor (to calc Tmin and Tmax) and a reference temperature time series Tref
			vcor <- nam==shreflist$Tcor[i]
			vref <- nam==shreflist$Tref[i]
			Tcor <- rmat$d[,vcor]
			Tref <- rmat$d[,vref]

			## dT is mean difference 
			dT <- mean(Tcor-Tref, na.rm=TRUE)
      dT <- round(dT,2)
			
			## Trefref is Tref with values only where Tcor has values as well
			Trefref <- Tref
			Trefref[is.na(Tcor)] <- NA 
      

			## calc uncertainty: U = sd(T-Tref) * (sd(Tref)/sd(Trefref)) 
			U <- sd(Tcor-Tref, na.rm=TRUE) * (sd(Tref, na.rm=TRUE)/sd(Trefref, na.rm=TRUE))
      U <- round(U,2)
      
			## calc total values and write back to dfs in list
			## matcor = mat + (matref - matrefref)
			## Tmin/max,SD  = Trefmax/min + dT +- U
		
			stat$mean[lin,vcor] <- round(mmat[vcor] + mmat[vref] - mean(Trefref, na.rm=TRUE),2)
			stat$max[lin,vcor] <- Tmax[vref]+ dT + U
			stat$min[lin,vcor] <- Tmin[vref]+ dT - U
			stat$q90[lin,vcor] <- q90[vref]+ dT + U
			stat$q10[lin,vcor] <- q10[vref]+ dT - U
			stat$stdv[lin,vcor] <- SD[vref]+ dT + U
		}
		}
		
	}

	## asign df as list with 'mean', 'max', 'min', 'stdv', 'count'
	dfn <- "MAT_" & site 
	assign(dfn,stat)
	if (save) write.csv (stat,"CSV_pro/"&dfn&".csv",row.names=FALSE)

	## **********************************************************
	## call freezing / thawing degree days (deg.day)
	print("calc degree days")
	degday <- deg.day(dm)
	## asign df as list with 'd' for data and 'q' for quality
	dfn <- "DeD_" & site 
	assign(dfn,degday)
	if (save) write.csv (degday,"CSV_pro/"&dfn&".csv",row.names=FALSE)



}#   END of main loop through the sites
##======================================================================================
##======================================================================================
## save workspace in file

## remove unused variables
rm(list=c("Tcor","Tmax","Tmin","Tref","SD","dfn","mmat","stat","d","site","correct","lsl"))
rm(list=c("lin","vcor","vref","U","dm","rmat","dT","nam","l"))

 ## save all dataframes for later use (for faster data access)
if (save) save(list=ls(), file = "RData/STATws.RData") 

##======================================================================================
##   end
##======================================================================================
