##==== 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
## change critaria for line style to abs(offset) > variation spread -- Jan 2014, Hasi
##
## Autors: 
## Hasi (Andreas Hasler), andreas_hasler@gmx.ch, BC Forest Dept.
## 
## Description:		Plots 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:			plot profiles
##
## 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
##======================================================================================

## set file name for meta data file
proflist <- "meta/SO_TO_proflist.csv"

## set plot mode: "png","pdf" or "scr" (for screen)
pl <- "scr"
allin1 <- T

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

##======================================================================================
##   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")

## 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=";")

## prepare plot
plfn <- "plots/PROallsites"
if (allin1) pldev(x=8,y=15,pld=pl,plfn=plfn)
## make rows and cols on sheet c(rows,cols):
if (allin1) par(mfrow=c(4,2),xpd=FALSE)


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

for (k in 1:lsl){	
	## get site name	
	site <- sites[k]
	sitn <- (slist$sitename[site==slist$site])[1]
	print("ploting "&site)
	
	## getting dataframe MAT of respective sensor assigned to "d"
	dfn <- "MAT_" & site 
	stat <- get(dfn)
	mmat <- stat$mean
	
	## get line number with totals	
	nam <- names(mmat)
	lin <- mmat$t == "tot"
	mmat <- as.numeric(mmat[lin,])

	## proflist for site
	lin <- proflist$site == site
	prlist <- proflist[lin,]
	nprof <- dim(prlist)[1]
	if (nprof > 0){
   ## only if profile defined... 
	  
	 ## getting dataframe spread of respective sensor assigned to "d"
	 dfn <- "SOTO_" & site 
	 soto <- get(dfn)  

##======================================================================================
##   ploting
##======================================================================================
	## define grafic device size
 	plfn <- paste("plots/PROpsite_",site,sep='')
	if (!allin1) pldev(x=6,y=6,pld=pl,plfn=plfn)

	## title 
	tlab <- "T-profil"
	title <- paste(sitn,tlab,sep=' ')
	## ploting variables
	xlim <- c(-4,4)
	if(site=="HUD") xlim <- c(-5,3)
	if(site=="POP") xlim <- c(-2,6)
	ylim <- c(-1,2)
	xlab <- "MAT (°C)"
	ylab <- "height / depth (m)"
	
	## point type
	pty <- c(1,19,2,6,17,0,5,15,18,3,4,8,7,9,10,11,12)

	##plot settings:
	col <- "black"
	
	## prep vector with location (or profil) names
	loc <- 1:nprof

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

		## prepare plotting vectors
		var <- match(prlist[i,2:4],nam)
		yval <- nam[var]
		if (!is.na(yval[2])) vn <- yval[2] else vn <- yval[3]
		vn <- unlist(strsplit(vn,"_"))[1]
		loc[i] <- "location "& vn
		y <- c(0,0,0)
		y[is.na(yval)] <- NA
		y <- y + c(1.25, 0,-1*prlist[i,5])
		x <- mmat[var] 
    vx <- as.numeric(soto[1:6,vn])
		
		## set lapse rate correction:		
		offs <- prlist$lapsrc[i]
		if (is.na(offs)) dT <- 0 else dT <- offs
		x[1] <- x[1] + dT
		
		## prepare line plot
		## split to different parts for separate ploting of line style
		xso <- x[1:2]
		xto <- x[2:3]
		xdT <- x[c(1,3)]
		yso <- y[1:2]
		yto <- y[2:3]
		ydT <- y[c(1,3)] 
		
		## make linetype vector c(SO, dT, TO)
    ## dependent on offset vs. variation
    ## if offset larger than variation in temp. difference make line solid
		lt <- as.numeric(abs(vx[1:3]) < (vx[4:6]/2))	+ 1
    
		## **************************
		## plot 
		## plot lines and points of profile
		if(i==1)plot(x,y,col=col,pch=pty[i],xlim=xlim,ylim=ylim,xlab=xlab, ylab=ylab, main=sitn)
		if(!is.na(vx[1]))lines(xso,yso,lty=lt[1])
		if(is.na(x[2])&!is.na(vx[2])) lines(xdT,ydT,lty=lt[2])  # plot dTline only if Tsurf missing
		if(!is.na(vx[3]))lines(xto,yto,lty=lt[3])
		points(x,y,col=col,pch=pty[i],cex=1.2)
 	}
 	abline(h=0,v=0,lty=1)
 	legend("topright",loc,pch=pty[1:nprof],pt.cex=1.3,bty="n")
 	## close ploting dev for png or pdf
	if ((pl!="scr") &!allin1) dev.off()
 	} ## end of "if nprof>0"


}#   END of main loop through the sites

if ((pl!="scr") & allin1) dev.off()

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