##==== 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:		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 plot mode: "png","pdf" or "scr" (for screen)
pl <- "png"

## save ANN_AMP file and dataframe?
save <- TRUE

##======================================================================================
##   preparation & import
##======================================================================================


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

## get list of sites:
sites <- c("HUD","GUN","PIN","NON","MID","TET","POP")
lsl <- length(sites) ## length of short list

## list with definitions of air temp sensor per site:
siteamp <- data.frame(	gun=1,
						mid=1)

## prepare ploting
## define grafic device size
plfn <- "plots/SOvsAA"
pldev(x=5,y=5,pld=pl,plfn=plfn)
## make rows and cols on sheet c(rows,cols):
par(mfrow=c(1,1),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 SOTO of respective sensor assigned to "soto"
	dfn <- "SOTO_" & site 
	soto <- get(dfn)
	cn <- dim(soto)[2] 
	typ <- soto[soto$var=="typ",2:cn]
	offset <- as.numeric(soto[1,2:cn])
  offset[is.na(offset)] <- as.numeric(soto[2,2:cn])[is.na(offset)]
	mso <- mean(offset,na.rm=TRUE)
	
	## getting dataframe DMT of respective sensor assigned to "dmt"
	dfn <- "DMT_" & site 
	dmt <- get(dfn)
	
	## calc annual amplitude
	aa <- ts.amp(dmt)
	siteaa <- as.numeric(aa[1])
	if (site=="TET") siteaa <- as.numeric(aa$N1_Tair)
	if (site=="POP") siteaa <- as.numeric(aa$N4_Tair)
	
	## get number of offset per site
	nprof <- length(offset)
	
	## make list with aa&k data-frames
	if (k==1) ANN_AMP <- list(0)
	ANN_AMP[site]<- siteaa 
	## make list with aa&k data-frames
	if (k==1) ALL_AMP <- list(0)
	ALL_AMP[[site]]<- aa 

##======================================================================================
##   ploting
##======================================================================================
	
	## ploting limits
	xlim <- c(7,16)
	ylim <- c(0,6)
  ylab <- "SO (°C)"
  xlab <- "annual amplitude (°C)"
	
	## colours for different surface types
	col <- as.character(typ)
	col[typ=="r"] <- "grey"
	col[typ=="s"] <- "brown"
	col[typ=="f"] <- "darkgreen"
	col[typ=="d"] <- "black"
	
	## point type
	pty <- c(1,19,2,6,17,0,5,15)
	
	## points for different surface types
	pt2 <- as.numeric(typ)
	pt2[typ=="r"] <- 18
	pt2[typ=="s"] <- 21
	pt2[typ=="f"] <- 17
	pt2[typ=="d"] <- 22
	

	## plot 
		
	if(k==1) plot(siteaa,mso,type="p",xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,cex=1.4)
	points(siteaa,mso,col=1,pch=pty[k],cex=1.4)
	for (i in 1:nprof) points(siteaa,offset[i],col="grey",pch=pt2[i],cex=0.8)
		

	
		
	## plot variation arrows:
	#arrows(mx-v, my, mx+v, my, length = 0.03, angle = 90, code = 3, lwd = 1)
	

		
}#   END of main loop through the sites
##======================================================================================

#abline(h=0,v=0,lty=1,lwd=0.5)	
legend("topleft",sites,pch=pty[1:lsl],pt.cex=1.3,bty="n")
legend("bottomleft",c("rock","debris","soil","forest"),pch=c(18,22,21,17),cex=1,pt.cex=0.8,bty="n",col="grey",text.col="darkgrey")

ANN_AMP <- ANN_AMP[2:(lsl+1)]

## close ploting dev for png or pdf
if (pl!="scr") dev.off()
if (save) {
	
	## remove unused variables
	rm(list=c("aa","dfn","i","k","pl","lsl","sitn","site","mso","soto","dmt","offset"))
	
 	## save all dataframes for later use (for faster data access)
	save(list=ls(), file = "RData/STATws.RData") 
	
	
	## save ANN_AMP csv file
	write.csv (ANN_AMP,"CSV_pro/ANN_AMP.csv",row.names=FALSE)
	write.csv (ALL_AMP,"CSV_pro/ALL_AMP.csv",row.names=FALSE)
	
}

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