##==== HEADER ========================================================================
## 
## Purpose:
## This is a script to plot temperatur data from British Columbia
## Ministry of Forest, Land and Natural Resources Operation
##
## Modification History:
## code generation -- April 2013, Hasi
##
## Autors: 
## Hasi (Andreas Hasler), andreas_hasler@gmx.ch, BC Forest Dept.
## 
## Description:		calc of SO and TO profiles for each field site.
##
## Dependencies:		Ttools.R (own toolbox for temperature timeseries treatment,
## 				(contained in STATws.R workspace)
##
## I/O, variables and units:
## input:		STATws.RData workspace (generated by procedure "stat_TTSMAT.R")
##				contains toolbox and different aggregates of time series
##				for details see header of "stat_TTSMAT.R"
##				SO_TO_proflist.csv = list of variables in each profile
##		
## output:		calculation of SO and TO and 
##
## units:		t is in inported as "mm-dd-yyyy HH:MM:SS AM/PM"
##				t is internaly in R POSIXct
##				temperature is in ?C
##						
##======================================================================================

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

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


## supress variation from MAT calc for sites:
supress <- c("","")

## use 10% / 90% quantiles (TRUE) or min max (FALSE)
quantiles <- FALSE


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

## load workspace with all data and tools 
if (!exists("STAWS")) load("RData/STATws.RData")

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

## save csv files?
save <- FALSE

## get list of sites:
sites <- levels(factor(slist$site))
lsl <- length(sites) ## length of short list

## getting list with profil definitions
proflist <- read.csv("meta/SO_TO_proflist.csv", as.is=TRUE)
if (dim(proflist)[2]==1) proflist <- read.csv("meta/SO_TO_proflist.csv", as.is=TRUE, sep=";")

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

for (k in 1:lsl){	
	## get site name	
	site <- sites[k]
	sitn <- (slist$sitename[site==slist$site])[1]
	print("processing SOTO for "&site)

 	## getting dataframe RMAT of respective sensor assigned to "d"
	dfn <- "RMAT_" & site 
	rmat <- get(dfn)
	rmat <- rmat$d
  
	## getting dataframe MAT of respective sensor assigned to "d"
	dfn <- "MAT_" & site 
	stat <- get(dfn)
	mmat <- stat$mean
	if (quantiles) 	xmax <- stat$q90 else xmax <- stat$max
	if (quantiles) xmin <- stat$q10 else xmin <- stat$min
	daycount <- stat$count
	
	## get totals	(mean/max/min) of MAT statistics
	nam <- names(mmat)
	lin <- mmat$t == "tot"
	mmat <- as.numeric(mmat[lin,])
	xmax <- as.numeric(xmax[lin,])
	xmin <- as.numeric(xmin[lin,])
	daycount <- as.numeric(daycount[lin,])

	## maximum of daycount
	dcmax <- max(daycount[2:length(daycount)],na.rm=TRUE)

	## proflist for site
	lin <- proflist$site == site
	prlist <- proflist[lin,]
	nprof <- dim(prlist)[1]
	if (nprof > 0){

##======================================================================================
##   calc
##======================================================================================
	
	## loop for profiles
	for(i in 1:nprof){

		## prepare vectors
		var <- match(prlist[i,2:4],nam)
		typ <- prlist$surftype[i]	
		yval <- nam[var]
		if (!is.na(yval[2])) vn <- yval[2] else vn <- yval[3]
		vn <- unlist(strsplit(vn,"_"))[1]
		y <- c(0,0,0)
		y[is.na(yval)] <- NA
		y <- y + c(1.25, 0,-1*prlist[i,5])
		x <- mmat[var] 
		xmaxi <- xmax[var]
		xmini <- xmin[var]
		
		## set lapse rate correction:		
		offs <- prlist$lapsrc[i]
		if (is.na(offs)) dT <- 0 else dT <- offs
		x[1] <- x[1] + dT
		
		## **************************
		## make variation ranges:
		if (!is.na(var[1])) T1 <- rmat[,var[1]] else T1 <- NA
		if (!is.na(var[2])) T2 <- rmat[,var[2]] else T2 <- NA
		if (!is.na(var[3])) T3 <- rmat[,var[3]] else T3 <- NA
		
		if (quantiles){
		  v1_2 <- (quant09(T1-T2,na.rm=T) - quant01(T1-T2,na.rm=T))
		  v1_3 <- (quant09(T1-T3,na.rm=T) - quant01(T1-T3,na.rm=T))
		  v2_3 <- (quant09(T2-T3,na.rm=T) - quant01(T2-T3,na.rm=T))
		} else {
		  v1_2 <- (max(T1-T2,na.rm=T) - min(T1-T2,na.rm=T))
		  v1_3 <- (max(T1-T3,na.rm=T) - min(T1-T3,na.rm=T))
		  v2_3 <- (max(T2-T3,na.rm=T) - min(T2-T3,na.rm=T))
		}
		
		## take spreads from MAT calc for short time series if not supressed (Gunnel):
		if(!max(supress==site,na.rm=TRUE)){
		  if (count(T1-T2) < (dcmax/2)) v1_2 <- ((xmaxi[1] -xmini[1])+ (xmaxi[2] -xmini[2]))
		  if (count(T1-T3) < (dcmax/2)) v1_3 <- ((xmaxi[1] -xmini[1])+ (xmaxi[3] -xmini[3]))  
		  if (count(T2-T3) < (dcmax/2)) v2_3 <- ((xmaxi[2] -xmini[2])+ (xmaxi[3] -xmini[3]))
		}
		
		v1_2 <- v1_2 + (dT/2)	## here we assume lapserate error beeing +-25% (-3.75 to -6.25 deg/km)  
		v1_3 <- v1_3 + (dT/2)  
    
		## vector with variations and their plot (means of plotvar to plot in between)
		v <- c(v1_2,v1_3,v2_3) 
    v[v==-Inf] <- NA
    print (paste("raw variation: ",v,sep=""))

    ## add uncertainty of measurements +- 0.3°C:
    v <- sqrt((v*v)+(0.6*0.6))
		print (paste("uncert.: ",v,sep=""))	
		
		## **************************

		## make vector 
		SO <- x[2]-x[1]
		dT <- x[3]-x[1]	
		TO <- x[3]-x[2]
		loc_i <- round(c(SO,dT,TO,v),2)
		## add type
		loc_i <- c(loc_i,typ) 
		
		if (i == 1) d <- data.frame(var = c("SO","dT","TO","vSO","vdT","vTO","typ")) 
		d[,i+1] <- loc_i
		names(d)[i+1] <- vn

 	}
 	
 	## here (nprof > 0) because no data for saxton (yet)
 	## estimate missing SO values
 	#for (i in 2:nprof+1){
 	#	dT<-as.numeric(d[3,i])
 	#	mTO <- mean(as.numeric(d[2,2:nprof+1]),na.rm=T)
 	#	if (site=="POP") mTO <-0
 	#	if (is.na(d[1,i])) d[1,i] <- dT - mTO 
 		
 	#}
 		
	## asign df as list with 'd' for data and 'q' for quality
	dfn <- "SOTO_" & site 
	assign(dfn,d)
	if (save) write.csv (d,"CSV_pro/"&dfn&".csv",row.names=FALSE)
 	
 	
 	}


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

## remove unused variables
rm(list=c("SO","TO","dT","x","y","i","k","loc_i","my","mx","v1_2","v2_3","v1_3","xmax","xmin","var"))
rm(list=c("xmaxi","xmini","nam","lin","nprof","dfn","rmat","stat","d","site","supress","vn","T1","T2","T3"))

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